Changeset 3869 for palm/trunk


Ignore:
Timestamp:
Apr 8, 2019 11:54:20 AM (2 years ago)
Author:
knoop
Message:

moving the furniture in the BCM around ;-)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/bulk_cloud_model_mod.f90

    r3786 r3869  
    2525! -----------------
    2626! $Id$
     27! moving the furniture around ;-)
     28!
     29! 3786 2019-03-06 16:58:03Z raasch
    2730! unsed variables removed
    2831!
     
    325328           bcm_check_data_output, &
    326329           bcm_check_data_output_pr, &
    327            bcm_header, &
    328330           bcm_init_arrays, &
    329331           bcm_init, &
     332           bcm_header, &
     333           bcm_actions, &
    330334           bcm_3d_data_averaging, &
    331335           bcm_data_output_2d, &
     
    336340           bcm_wrd_global, &
    337341           bcm_wrd_local, &
    338            bcm_actions, &
    339342           calc_liquid_water_content
    340343
     
    370373    END INTERFACE bcm_check_data_output_pr
    371374
     375    INTERFACE bcm_init_arrays
     376       MODULE PROCEDURE bcm_init_arrays
     377    END INTERFACE bcm_init_arrays
     378
     379    INTERFACE bcm_init
     380       MODULE PROCEDURE bcm_init
     381    END INTERFACE bcm_init
     382
    372383    INTERFACE bcm_header
    373384       MODULE PROCEDURE bcm_header
    374385    END INTERFACE bcm_header
    375 
    376     INTERFACE bcm_init_arrays
    377        MODULE PROCEDURE bcm_init_arrays
    378     END INTERFACE bcm_init_arrays
    379 
    380     INTERFACE bcm_init
    381        MODULE PROCEDURE bcm_init
    382     END INTERFACE bcm_init
    383 
    384     INTERFACE bcm_3d_data_averaging
    385        MODULE PROCEDURE bcm_3d_data_averaging
    386     END INTERFACE bcm_3d_data_averaging
    387 
    388     INTERFACE bcm_data_output_2d
    389        MODULE PROCEDURE bcm_data_output_2d
    390     END INTERFACE bcm_data_output_2d
    391 
    392     INTERFACE bcm_data_output_3d
    393        MODULE PROCEDURE bcm_data_output_3d
    394     END INTERFACE bcm_data_output_3d
    395 
    396     INTERFACE bcm_swap_timelevel
    397        MODULE PROCEDURE bcm_swap_timelevel
    398     END INTERFACE bcm_swap_timelevel
    399 
    400     INTERFACE bcm_rrd_global
    401        MODULE PROCEDURE bcm_rrd_global
    402     END INTERFACE bcm_rrd_global
    403 
    404     INTERFACE bcm_rrd_local
    405        MODULE PROCEDURE bcm_rrd_local
    406     END INTERFACE bcm_rrd_local
    407 
    408     INTERFACE bcm_wrd_global
    409        MODULE PROCEDURE bcm_wrd_global
    410     END INTERFACE bcm_wrd_global
    411 
    412     INTERFACE bcm_wrd_local
    413        MODULE PROCEDURE bcm_wrd_local
    414     END INTERFACE bcm_wrd_local
    415386
    416387    INTERFACE bcm_actions
     
    419390    END INTERFACE bcm_actions
    420391
    421     INTERFACE adjust_cloud
    422        MODULE PROCEDURE adjust_cloud
    423        MODULE PROCEDURE adjust_cloud_ij
    424     END INTERFACE adjust_cloud
    425 
    426     INTERFACE activation
    427        MODULE PROCEDURE activation
    428        MODULE PROCEDURE activation_ij
    429     END INTERFACE activation
    430 
    431     INTERFACE condensation
    432        MODULE PROCEDURE condensation
    433        MODULE PROCEDURE condensation_ij
    434     END INTERFACE condensation
    435 
    436     INTERFACE autoconversion
    437        MODULE PROCEDURE autoconversion
    438        MODULE PROCEDURE autoconversion_ij
    439     END INTERFACE autoconversion
    440 
    441     INTERFACE autoconversion_kessler
    442        MODULE PROCEDURE autoconversion_kessler
    443        MODULE PROCEDURE autoconversion_kessler_ij
    444     END INTERFACE autoconversion_kessler
    445 
    446     INTERFACE accretion
    447        MODULE PROCEDURE accretion
    448        MODULE PROCEDURE accretion_ij
    449     END INTERFACE accretion
    450 
    451     INTERFACE selfcollection_breakup
    452        MODULE PROCEDURE selfcollection_breakup
    453        MODULE PROCEDURE selfcollection_breakup_ij
    454     END INTERFACE selfcollection_breakup
    455 
    456     INTERFACE evaporation_rain
    457        MODULE PROCEDURE evaporation_rain
    458        MODULE PROCEDURE evaporation_rain_ij
    459     END INTERFACE evaporation_rain
    460 
    461     INTERFACE sedimentation_cloud
    462        MODULE PROCEDURE sedimentation_cloud
    463        MODULE PROCEDURE sedimentation_cloud_ij
    464     END INTERFACE sedimentation_cloud
    465 
    466     INTERFACE sedimentation_rain
    467        MODULE PROCEDURE sedimentation_rain
    468        MODULE PROCEDURE sedimentation_rain_ij
    469     END INTERFACE sedimentation_rain
    470 
    471     INTERFACE calc_precipitation_amount
    472        MODULE PROCEDURE calc_precipitation_amount
    473        MODULE PROCEDURE calc_precipitation_amount_ij
    474     END INTERFACE calc_precipitation_amount
    475 
    476     INTERFACE supersaturation
    477        MODULE PROCEDURE supersaturation
    478     END INTERFACE supersaturation
     392    INTERFACE bcm_swap_timelevel
     393       MODULE PROCEDURE bcm_swap_timelevel
     394    END INTERFACE bcm_swap_timelevel
     395
     396    INTERFACE bcm_3d_data_averaging
     397       MODULE PROCEDURE bcm_3d_data_averaging
     398    END INTERFACE bcm_3d_data_averaging
     399
     400    INTERFACE bcm_data_output_2d
     401       MODULE PROCEDURE bcm_data_output_2d
     402    END INTERFACE bcm_data_output_2d
     403
     404    INTERFACE bcm_data_output_3d
     405       MODULE PROCEDURE bcm_data_output_3d
     406    END INTERFACE bcm_data_output_3d
     407
     408    INTERFACE bcm_rrd_global
     409       MODULE PROCEDURE bcm_rrd_global
     410    END INTERFACE bcm_rrd_global
     411
     412    INTERFACE bcm_rrd_local
     413       MODULE PROCEDURE bcm_rrd_local
     414    END INTERFACE bcm_rrd_local
     415
     416    INTERFACE bcm_wrd_global
     417       MODULE PROCEDURE bcm_wrd_global
     418    END INTERFACE bcm_wrd_global
     419
     420    INTERFACE bcm_wrd_local
     421       MODULE PROCEDURE bcm_wrd_local
     422    END INTERFACE bcm_wrd_local
    479423
    480424    INTERFACE calc_liquid_water_content
     
    989933! Description:
    990934! ------------
     935!> Header output for bulk cloud module
     936!------------------------------------------------------------------------------!
     937    SUBROUTINE bcm_header ( io )
     938
     939
     940       IMPLICIT NONE
     941
     942       INTEGER(iwp), INTENT(IN) ::  io  !< Unit of the output file
     943
     944!
     945!--    Write bulk cloud module header
     946       WRITE ( io, 1 )
     947
     948       WRITE ( io, 2 )
     949       WRITE ( io, 3 )
     950
     951       IF ( microphysics_kessler )  THEN
     952          WRITE ( io, 4 ) 'Kessler-Scheme'
     953       ENDIF
     954
     955       IF ( microphysics_seifert )  THEN
     956          WRITE ( io, 4 ) 'Seifert-Beheng-Scheme'
     957          IF ( cloud_water_sedimentation )  WRITE ( io, 5 )
     958          IF ( collision_turbulence )  WRITE ( io, 6 )
     959          IF ( ventilation_effect )  WRITE ( io, 7 )
     960          IF ( limiter_sedimentation )  WRITE ( io, 8 )
     961       ENDIF
     962
     963       WRITE ( io, 20 )
     964       WRITE ( io, 21 ) surface_pressure
     965       WRITE ( io, 22 ) r_d
     966       WRITE ( io, 23 ) rho_surface
     967       WRITE ( io, 24 ) c_p
     968       WRITE ( io, 25 ) l_v
     969
     970       IF ( microphysics_seifert )  THEN
     971          WRITE ( io, 26 ) 1.0E-6_wp * nc_const
     972          WRITE ( io, 27 ) c_sedimentation
     973       ENDIF
     974
     975
     9761   FORMAT ( //' Bulk cloud module information:'/ &
     977               ' ------------------------------------------'/ )
     9782   FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/ &
     979             '    total water content is used.' )
     9803   FORMAT ( '--> Condensation is parameterized via 0% - or 100% scheme.' )
     9814   FORMAT ( '--> Precipitation parameterization via ', A )
     982
     9835   FORMAT ( '--> Cloud water sedimentation parameterization via Stokes law' )
     9846   FORMAT ( '--> Turbulence effects on precipitation process' )
     9857   FORMAT ( '--> Ventilation effects on evaporation of rain drops' )
     9868   FORMAT ( '--> Slope limiter used for sedimentation process' )
     987
     98820  FORMAT ( '--> Essential parameters:' )
     98921  FORMAT ( '       Surface pressure             :   p_0   = ', F7.2, ' hPa')
     99022  FORMAT ( '       Gas constant                 :   R     = ', F5.1, ' J/(kg K)')
     99123  FORMAT ( '       Density of air               :   rho_0 = ', F6.3, ' kg/m**3')
     99224  FORMAT ( '       Specific heat cap.           :   c_p   = ', F6.1, ' J/(kg K)')
     99325  FORMAT ( '       Vapourization heat           :   L_v   = ', E9.2, ' J/kg')
     99426  FORMAT ( '       Droplet density              :   N_c   = ', F6.1, ' 1/cm**3' )
     99527  FORMAT ( '       Sedimentation Courant number :   C_s   = ', F4.1 )
     996
     997
     998    END SUBROUTINE bcm_header
     999
     1000
     1001!------------------------------------------------------------------------------!
     1002! Description:
     1003! ------------
     1004!> Control of microphysics for all grid points
     1005!------------------------------------------------------------------------------!
     1006    SUBROUTINE bcm_actions
     1007
     1008       IMPLICIT NONE
     1009
     1010       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
     1011!
     1012!--       Calculate vertical profile of the hydrostatic pressure (hyp)
     1013          hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
     1014          d_exner = exner_function_invers(hyp)
     1015          exner = 1.0_wp / exner_function_invers(hyp)
     1016          hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
     1017!
     1018!--       Compute reference density
     1019          rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
     1020       ENDIF
     1021
     1022!
     1023!--    Compute length of time step
     1024       IF ( call_microphysics_at_all_substeps )  THEN
     1025          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
     1026       ELSE
     1027          dt_micro = dt_3d
     1028       ENDIF
     1029
     1030!
     1031!--    Reset precipitation rate
     1032       IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
     1033
     1034!
     1035!--    Compute cloud physics
     1036       IF ( microphysics_kessler )  THEN
     1037
     1038          CALL autoconversion_kessler
     1039          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
     1040
     1041       ELSEIF ( microphysics_seifert )  THEN
     1042
     1043          CALL adjust_cloud
     1044          IF ( microphysics_morrison )  CALL activation
     1045          IF ( microphysics_morrison )  CALL condensation
     1046          CALL autoconversion
     1047          CALL accretion
     1048          CALL selfcollection_breakup
     1049          CALL evaporation_rain
     1050          CALL sedimentation_rain
     1051          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
     1052
     1053       ENDIF
     1054
     1055       CALL calc_precipitation_amount
     1056
     1057    END SUBROUTINE bcm_actions
     1058
     1059
     1060!------------------------------------------------------------------------------!
     1061! Description:
     1062! ------------
     1063!> Control of microphysics for grid points i,j
     1064!------------------------------------------------------------------------------!
     1065
     1066    SUBROUTINE bcm_actions_ij( i, j )
     1067
     1068       IMPLICIT NONE
     1069
     1070       INTEGER(iwp) ::  i                 !<
     1071       INTEGER(iwp) ::  j                 !<
     1072
     1073       IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
     1074!
     1075!--       Calculate vertical profile of the hydrostatic pressure (hyp)
     1076          hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
     1077          d_exner = exner_function_invers(hyp)
     1078          exner = 1.0_wp / exner_function_invers(hyp)
     1079          hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
     1080!
     1081!--       Compute reference density
     1082          rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
     1083       ENDIF
     1084
     1085!
     1086!--    Compute length of time step
     1087       IF ( call_microphysics_at_all_substeps )  THEN
     1088          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
     1089       ELSE
     1090          dt_micro = dt_3d
     1091       ENDIF
     1092!
     1093!--    Reset precipitation rate
     1094       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
     1095
     1096!
     1097!--    Compute cloud physics
     1098       IF( microphysics_kessler )  THEN
     1099
     1100          CALL autoconversion_kessler_ij( i,j )
     1101          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
     1102
     1103       ELSEIF ( microphysics_seifert )  THEN
     1104
     1105          CALL adjust_cloud_ij( i,j )
     1106          IF ( microphysics_morrison ) CALL activation_ij( i,j )
     1107          IF ( microphysics_morrison ) CALL condensation_ij( i,j )
     1108          CALL autoconversion_ij( i,j )
     1109          CALL accretion_ij( i,j )
     1110          CALL selfcollection_breakup_ij( i,j )
     1111          CALL evaporation_rain_ij( i,j )
     1112          CALL sedimentation_rain_ij( i,j )
     1113          IF ( cloud_water_sedimentation )  CALL sedimentation_cloud_ij( i,j )
     1114
     1115       ENDIF
     1116
     1117       CALL calc_precipitation_amount_ij( i,j )
     1118
     1119    END SUBROUTINE bcm_actions_ij
     1120
     1121
     1122!------------------------------------------------------------------------------!
     1123! Description:
     1124! ------------
    9911125!> Swapping of timelevels
    9921126!------------------------------------------------------------------------------!
     
    10281162
    10291163    END SUBROUTINE bcm_swap_timelevel
    1030 
    1031 
    1032 !------------------------------------------------------------------------------!
    1033 ! Description:
    1034 ! ------------
    1035 !> Header output for bulk cloud module
    1036 !------------------------------------------------------------------------------!
    1037     SUBROUTINE bcm_header ( io )
    1038 
    1039 
    1040        IMPLICIT NONE
    1041 
    1042        INTEGER(iwp), INTENT(IN) ::  io  !< Unit of the output file
    1043 
    1044 !
    1045 !--    Write bulk cloud module header
    1046        WRITE ( io, 1 )
    1047 
    1048        WRITE ( io, 2 )
    1049        WRITE ( io, 3 )
    1050 
    1051        IF ( microphysics_kessler )  THEN
    1052           WRITE ( io, 4 ) 'Kessler-Scheme'
    1053        ENDIF
    1054 
    1055        IF ( microphysics_seifert )  THEN
    1056           WRITE ( io, 4 ) 'Seifert-Beheng-Scheme'
    1057           IF ( cloud_water_sedimentation )  WRITE ( io, 5 )
    1058           IF ( collision_turbulence )  WRITE ( io, 6 )
    1059           IF ( ventilation_effect )  WRITE ( io, 7 )
    1060           IF ( limiter_sedimentation )  WRITE ( io, 8 )
    1061        ENDIF
    1062 
    1063        WRITE ( io, 20 )
    1064        WRITE ( io, 21 ) surface_pressure
    1065        WRITE ( io, 22 ) r_d
    1066        WRITE ( io, 23 ) rho_surface
    1067        WRITE ( io, 24 ) c_p
    1068        WRITE ( io, 25 ) l_v
    1069 
    1070        IF ( microphysics_seifert )  THEN
    1071           WRITE ( io, 26 ) 1.0E-6_wp * nc_const
    1072           WRITE ( io, 27 ) c_sedimentation
    1073        ENDIF
    1074 
    1075 
    1076 1   FORMAT ( //' Bulk cloud module information:'/ &
    1077                ' ------------------------------------------'/ )
    1078 2   FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/ &
    1079              '    total water content is used.' )
    1080 3   FORMAT ( '--> Condensation is parameterized via 0% - or 100% scheme.' )
    1081 4   FORMAT ( '--> Precipitation parameterization via ', A )
    1082 
    1083 5   FORMAT ( '--> Cloud water sedimentation parameterization via Stokes law' )
    1084 6   FORMAT ( '--> Turbulence effects on precipitation process' )
    1085 7   FORMAT ( '--> Ventilation effects on evaporation of rain drops' )
    1086 8   FORMAT ( '--> Slope limiter used for sedimentation process' )
    1087 
    1088 20  FORMAT ( '--> Essential parameters:' )
    1089 21  FORMAT ( '       Surface pressure             :   p_0   = ', F7.2, ' hPa')
    1090 22  FORMAT ( '       Gas constant                 :   R     = ', F5.1, ' J/(kg K)')
    1091 23  FORMAT ( '       Density of air               :   rho_0 = ', F6.3, ' kg/m**3')
    1092 24  FORMAT ( '       Specific heat cap.           :   c_p   = ', F6.1, ' J/(kg K)')
    1093 25  FORMAT ( '       Vapourization heat           :   L_v   = ', E9.2, ' J/kg')
    1094 26  FORMAT ( '       Droplet density              :   N_c   = ', F6.1, ' 1/cm**3' )
    1095 27  FORMAT ( '       Sedimentation Courant number :   C_s   = ', F4.1 )
    1096 
    1097 
    1098     END SUBROUTINE bcm_header
    10991164
    11001165
     
    19932058    END SUBROUTINE bcm_wrd_local
    19942059
    1995 
    1996 !------------------------------------------------------------------------------!
    1997 ! Description:
    1998 ! ------------
    1999 !> Control of microphysics for all grid points
    2000 !------------------------------------------------------------------------------!
    2001     SUBROUTINE bcm_actions
    2002 
    2003        IMPLICIT NONE
    2004 
    2005        IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
    2006 !
    2007 !--       Calculate vertical profile of the hydrostatic pressure (hyp)
    2008           hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
    2009           d_exner = exner_function_invers(hyp)
    2010           exner = 1.0_wp / exner_function_invers(hyp)
    2011           hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
    2012 !
    2013 !--       Compute reference density
    2014           rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
    2015        ENDIF
    2016 
    2017 !
    2018 !--    Compute length of time step
    2019        IF ( call_microphysics_at_all_substeps )  THEN
    2020           dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
    2021        ELSE
    2022           dt_micro = dt_3d
    2023        ENDIF
    2024 
    2025 !
    2026 !--    Reset precipitation rate
    2027        IF ( intermediate_timestep_count == 1 )  prr = 0.0_wp
    2028 
    2029 !
    2030 !--    Compute cloud physics
    2031        IF ( microphysics_kessler )  THEN
    2032 
    2033           CALL autoconversion_kessler
    2034           IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
    2035 
    2036        ELSEIF ( microphysics_seifert )  THEN
    2037 
    2038           CALL adjust_cloud
    2039           IF ( microphysics_morrison )  CALL activation
    2040           IF ( microphysics_morrison )  CALL condensation
    2041           CALL autoconversion
    2042           CALL accretion
    2043           CALL selfcollection_breakup
    2044           CALL evaporation_rain
    2045           CALL sedimentation_rain
    2046           IF ( cloud_water_sedimentation )  CALL sedimentation_cloud
    2047 
    2048        ENDIF
    2049 
    2050        CALL calc_precipitation_amount
    2051 
    2052     END SUBROUTINE bcm_actions
    2053 
    20542060!------------------------------------------------------------------------------!
    20552061! Description:
     
    21072113
    21082114    END SUBROUTINE adjust_cloud
     2115
     2116!------------------------------------------------------------------------------!
     2117! Description:
     2118! ------------
     2119!> Adjust number of raindrops to avoid nonlinear effects in
     2120!> sedimentation and evaporation of rain drops due to too small or
     2121!> too big weights of rain drops (Stevens and Seifert, 2008).
     2122!> The same procedure is applied to cloud droplets if they are determined
     2123!> prognostically. Call for grid point i,j
     2124!------------------------------------------------------------------------------!
     2125    SUBROUTINE adjust_cloud_ij( i, j )
     2126
     2127       IMPLICIT NONE
     2128
     2129       INTEGER(iwp) ::  i                 !<
     2130       INTEGER(iwp) ::  j                 !<
     2131       INTEGER(iwp) ::  k                 !<
     2132
     2133       REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
     2134
     2135       DO  k = nzb+1, nzt
     2136!
     2137!--       Predetermine flag to mask topography
     2138          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     2139
     2140          IF ( qr(k,j,i) <= eps_sb )  THEN
     2141             qr(k,j,i) = 0.0_wp
     2142             nr(k,j,i) = 0.0_wp
     2143          ELSE
     2144!
     2145!--          Adjust number of raindrops to avoid nonlinear effects in
     2146!--          sedimentation and evaporation of rain drops due to too small or
     2147!--          too big weights of rain drops (Stevens and Seifert, 2008).
     2148             IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
     2149                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
     2150             ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
     2151                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
     2152             ENDIF
     2153
     2154          ENDIF
     2155
     2156          IF ( microphysics_morrison ) THEN
     2157             IF ( qc(k,j,i) <= eps_sb )  THEN
     2158                qc(k,j,i) = 0.0_wp
     2159                nc(k,j,i) = 0.0_wp
     2160             ELSE
     2161                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
     2162                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xamin * flag
     2163                ENDIF
     2164             ENDIF
     2165          ENDIF
     2166
     2167       ENDDO
     2168
     2169    END SUBROUTINE adjust_cloud_ij
    21092170
    21102171!------------------------------------------------------------------------------!
     
    22142275    END SUBROUTINE activation
    22152276
     2277!------------------------------------------------------------------------------!
     2278! Description:
     2279! ------------
     2280!> Calculate number of activated condensation nucleii after simple activation
     2281!> scheme of Twomey, 1959.
     2282!------------------------------------------------------------------------------!
     2283    SUBROUTINE activation_ij( i, j )
     2284
     2285       IMPLICIT NONE
     2286
     2287       INTEGER(iwp) ::  i                 !<
     2288       INTEGER(iwp) ::  j                 !<
     2289       INTEGER(iwp) ::  k                 !<
     2290
     2291       REAL(wp)     ::  activ             !<
     2292       REAL(wp)     ::  afactor           !<
     2293       REAL(wp)     ::  beta_act          !<
     2294       REAL(wp)     ::  bfactor           !<
     2295       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
     2296       REAL(wp)     ::  k_act             !<
     2297       REAL(wp)     ::  n_act             !<
     2298       REAL(wp)     ::  n_ccn             !<
     2299       REAL(wp)     ::  s_0               !<
     2300       REAL(wp)     ::  sat_max           !<
     2301       REAL(wp)     ::  sigma             !<
     2302       REAL(wp)     ::  sigma_act         !<
     2303
     2304       DO  k = nzb+1, nzt
     2305!
     2306!--       Predetermine flag to mask topography
     2307          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     2308!
     2309!--       Call calculation of supersaturation
     2310          CALL supersaturation ( i, j, k )
     2311!
     2312!--       Prescribe parameters for activation
     2313!--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
     2314          k_act  = 0.7_wp
     2315          activ  = 0.0_wp
     2316
     2317          IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
     2318!
     2319!--          Compute the number of activated Aerosols
     2320!--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
     2321             n_act     = na_init * sat**k_act
     2322!
     2323!--          Compute the number of cloud droplets
     2324!--          (see: Morrison + Grabowski, 2007, JAS, 64)
     2325!            activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
     2326
     2327!
     2328!--          Compute activation rate after Khairoutdinov and Kogan
     2329!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
     2330             sat_max = 0.8_wp / 100.0_wp
     2331             activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN             &
     2332                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /         &
     2333                  dt_micro
     2334
     2335             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
     2336          ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk )  THEN
     2337!
     2338!--          Curvature effect (afactor) with surface tension
     2339!--          parameterization by Straka (2009)
     2340             sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
     2341             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
     2342!
     2343!--          Solute effect (bfactor)
     2344             bfactor = vanthoff * molecular_weight_of_water *                  &
     2345                  rho_s / ( molecular_weight_of_solute * rho_l )
     2346
     2347!
     2348!--          Prescribe power index that describes the soluble fraction
     2349!--          of an aerosol particle (beta).
     2350!--          (see: Morrison + Grabowski, 2007, JAS, 64)
     2351             beta_act  = 0.5_wp
     2352             sigma_act = sigma_bulk**( 1.0_wp + beta_act )
     2353!
     2354!--          Calculate mean geometric supersaturation (s_0) with
     2355!--          parameterization by Khvorostyanov and Curry (2006)
     2356             s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *         &
     2357               ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
     2358
     2359!
     2360!--          Calculate number of activated CCN as a function of
     2361!--          supersaturation and taking Koehler theory into account
     2362!--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
     2363             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                    &
     2364                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
     2365             activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
     2366
     2367             nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)
     2368          ENDIF
     2369
     2370       ENDDO
     2371
     2372    END SUBROUTINE activation_ij
     2373
    22162374
    22172375!------------------------------------------------------------------------------!
     
    22972455
    22982456    END SUBROUTINE condensation
     2457
     2458!------------------------------------------------------------------------------!
     2459! Description:
     2460! ------------
     2461!> Calculate condensation rate for cloud water content (after Khairoutdinov and
     2462!> Kogan, 2000).
     2463!------------------------------------------------------------------------------!
     2464    SUBROUTINE condensation_ij( i, j )
     2465
     2466       IMPLICIT NONE
     2467
     2468       INTEGER(iwp) ::  i                 !<
     2469       INTEGER(iwp) ::  j                 !<
     2470       INTEGER(iwp) ::  k                 !<
     2471
     2472       REAL(wp)     ::  cond              !<
     2473       REAL(wp)     ::  cond_max          !<
     2474       REAL(wp)     ::  dc                !<
     2475       REAL(wp)     ::  evap              !<
     2476       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
     2477       REAL(wp)     ::  g_fac             !<
     2478       REAL(wp)     ::  nc_0              !<
     2479       REAL(wp)     ::  temp              !<
     2480       REAL(wp)     ::  xc                !<
     2481
     2482
     2483       DO  k = nzb+1, nzt
     2484!
     2485!--       Predetermine flag to mask topography
     2486          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     2487!
     2488!--       Call calculation of supersaturation
     2489          CALL supersaturation ( i, j, k )
     2490!
     2491!--       Actual temperature:
     2492          temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
     2493
     2494          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
     2495                              l_v / ( thermal_conductivity_l * temp )    &
     2496                              + r_v * temp / ( diff_coeff_l * e_s )      &
     2497                            )
     2498!
     2499!--       Mean weight of cloud drops
     2500          IF ( nc(k,j,i) <= 0.0_wp) CYCLE
     2501          xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
     2502!
     2503!--       Weight averaged diameter of cloud drops:
     2504          dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
     2505!
     2506!--       Integral diameter of cloud drops
     2507          nc_0 = nc(k,j,i) * dc
     2508!
     2509!--       Condensation needs only to be calculated in supersaturated regions
     2510          IF ( sat > 0.0_wp )  THEN
     2511!
     2512!--          Condensation rate of cloud water content
     2513!--          after KK scheme.
     2514!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
     2515             cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
     2516             cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
     2517             cond      = MIN( cond, cond_max / dt_micro )
     2518
     2519             qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
     2520          ELSEIF ( sat < 0.0_wp ) THEN
     2521             evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
     2522             evap      = MAX( evap, -qc(k,j,i) / dt_micro )
     2523
     2524             qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag
     2525          ENDIF
     2526       ENDDO
     2527
     2528    END SUBROUTINE condensation_ij
    22992529
    23002530
     
    24322662! Description:
    24332663! ------------
     2664!> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
     2665!------------------------------------------------------------------------------!
     2666    SUBROUTINE autoconversion_ij( i, j )
     2667
     2668       IMPLICIT NONE
     2669
     2670       INTEGER(iwp) ::  i                 !<
     2671       INTEGER(iwp) ::  j                 !<
     2672       INTEGER(iwp) ::  k                 !<
     2673
     2674       REAL(wp)     ::  alpha_cc          !<
     2675       REAL(wp)     ::  autocon           !<
     2676       REAL(wp)     ::  dissipation       !<
     2677       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
     2678       REAL(wp)     ::  k_au              !<
     2679       REAL(wp)     ::  l_mix             !<
     2680       REAL(wp)     ::  nc_auto           !<
     2681       REAL(wp)     ::  nu_c              !<
     2682       REAL(wp)     ::  phi_au            !<
     2683       REAL(wp)     ::  r_cc              !<
     2684       REAL(wp)     ::  rc                !<
     2685       REAL(wp)     ::  re_lambda         !<
     2686       REAL(wp)     ::  sigma_cc          !<
     2687       REAL(wp)     ::  tau_cloud         !<
     2688       REAL(wp)     ::  xc                !<
     2689
     2690       DO  k = nzb+1, nzt
     2691!
     2692!--       Predetermine flag to mask topography
     2693          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     2694          IF ( microphysics_morrison ) THEN
     2695             nc_auto = nc(k,j,i)
     2696          ELSE
     2697             nc_auto = nc_const
     2698          ENDIF
     2699
     2700          IF ( qc(k,j,i) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
     2701
     2702             k_au = k_cc / ( 20.0_wp * x0 )
     2703!
     2704!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
     2705!--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
     2706             tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ),     &
     2707                              0.0_wp )
     2708!
     2709!--          Universal function for autoconversion process
     2710!--          (Seifert and Beheng, 2006):
     2711             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
     2712!
     2713!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
     2714!--          (Use constant nu_c = 1.0_wp instead?)
     2715             nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
     2716!
     2717!--          Mean weight of cloud droplets:
     2718             xc = hyrho(k) * qc(k,j,i) / nc_auto
     2719!
     2720!--          Parameterized turbulence effects on autoconversion (Seifert,
     2721!--          Nuijens and Stevens, 2010)
     2722             IF ( collision_turbulence )  THEN
     2723!
     2724!--             Weight averaged radius of cloud droplets:
     2725                rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
     2726
     2727                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
     2728                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
     2729                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
     2730!
     2731!--             Mixing length (neglecting distance to ground and stratification)
     2732                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
     2733!
     2734!--             Limit dissipation rate according to Seifert, Nuijens and
     2735!--             Stevens (2010)
     2736                dissipation = MIN( 0.06_wp, diss(k,j,i) )
     2737!
     2738!--             Compute Taylor-microscale Reynolds number:
     2739                re_lambda = 6.0_wp / 11.0_wp *                                 &
     2740                            ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
     2741                            SQRT( 15.0_wp / kin_vis_air ) *                    &
     2742                            dissipation**( 1.0_wp / 6.0_wp )
     2743!
     2744!--             The factor of 1.0E4 is needed to convert the dissipation rate
     2745!--             from m2 s-3 to cm2 s-3.
     2746                k_au = k_au * ( 1.0_wp +                                       &
     2747                                dissipation * 1.0E4_wp *                       &
     2748                                ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
     2749                                ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
     2750                                                  sigma_cc )**2                &
     2751                                                ) + beta_cc                    &
     2752                                )                                              &
     2753                              )
     2754             ENDIF
     2755!
     2756!--          Autoconversion rate (Seifert and Beheng, 2006):
     2757             autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
     2758                       ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *            &
     2759                       ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
     2760                       rho_surface
     2761             autocon = MIN( autocon, qc(k,j,i) / dt_micro )
     2762
     2763             qr(k,j,i) = qr(k,j,i) + autocon * dt_micro                 * flag
     2764             qc(k,j,i) = qc(k,j,i) - autocon * dt_micro                 * flag
     2765             nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag
     2766             IF ( microphysics_morrison ) THEN
     2767                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *                &
     2768                                    autocon / x0 * hyrho(k) * dt_micro * flag )
     2769             ENDIF
     2770
     2771          ENDIF
     2772
     2773       ENDDO
     2774
     2775    END SUBROUTINE autoconversion_ij
     2776
     2777
     2778!------------------------------------------------------------------------------!
     2779! Description:
     2780! ------------
    24342781!> Autoconversion process (Kessler, 1969).
    24352782!------------------------------------------------------------------------------!
     
    24772824
    24782825   END SUBROUTINE autoconversion_kessler
     2826
     2827!------------------------------------------------------------------------------!
     2828! Description:
     2829! ------------
     2830!> Autoconversion process (Kessler, 1969).
     2831!------------------------------------------------------------------------------!
     2832    SUBROUTINE autoconversion_kessler_ij( i, j )
     2833
     2834
     2835       IMPLICIT NONE
     2836
     2837       INTEGER(iwp) ::  i      !<
     2838       INTEGER(iwp) ::  j      !<
     2839       INTEGER(iwp) ::  k      !<
     2840       INTEGER(iwp) ::  k_wall !< topography top index
     2841
     2842       REAL(wp)    ::  dqdt_precip !<
     2843       REAL(wp)    ::  flag              !< flag to indicate first grid level above surface
     2844
     2845!
     2846!--    Determine vertical index of topography top
     2847       k_wall = get_topography_top_index_ji( j, i, 's' )
     2848       DO  k = nzb+1, nzt
     2849!
     2850!--       Predetermine flag to mask topography
     2851          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     2852
     2853          IF ( qc(k,j,i) > ql_crit )  THEN
     2854             dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
     2855          ELSE
     2856             dqdt_precip = 0.0_wp
     2857          ENDIF
     2858
     2859          qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
     2860          q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
     2861          pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k)  &
     2862                                                                 * flag
     2863
     2864!
     2865!--       Compute the rain rate (stored on surface grid point)
     2866          prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
     2867
     2868       ENDDO
     2869
     2870    END SUBROUTINE autoconversion_kessler_ij
    24792871
    24802872
     
    25642956    END SUBROUTINE accretion
    25652957
     2958!------------------------------------------------------------------------------!
     2959! Description:
     2960! ------------
     2961!> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
     2962!------------------------------------------------------------------------------!
     2963    SUBROUTINE accretion_ij( i, j )
     2964
     2965       IMPLICIT NONE
     2966
     2967       INTEGER(iwp) ::  i                 !<
     2968       INTEGER(iwp) ::  j                 !<
     2969       INTEGER(iwp) ::  k                 !<
     2970
     2971       REAL(wp)     ::  accr              !<
     2972       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
     2973       REAL(wp)     ::  k_cr              !<
     2974       REAL(wp)     ::  nc_accr           !<
     2975       REAL(wp)     ::  phi_ac            !<
     2976       REAL(wp)     ::  tau_cloud         !<
     2977       REAL(wp)     ::  xc                !<
     2978
     2979
     2980       DO  k = nzb+1, nzt
     2981!
     2982!--       Predetermine flag to mask topography
     2983          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     2984          IF ( microphysics_morrison ) THEN
     2985             nc_accr = nc(k,j,i)
     2986          ELSE
     2987             nc_accr = nc_const
     2988          ENDIF
     2989
     2990          IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )  .AND.  &
     2991               ( nc_accr > eps_mr ) )  THEN
     2992!
     2993!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
     2994             tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
     2995!
     2996!--          Universal function for accretion process
     2997!--          (Seifert and Beheng, 2001):
     2998             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
     2999
     3000!
     3001!--          Mean weight of cloud drops
     3002             xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
     3003!
     3004!--          Parameterized turbulence effects on autoconversion (Seifert,
     3005!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
     3006!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
     3007             IF ( collision_turbulence )  THEN
     3008                k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
     3009                                 MIN( 600.0_wp,                                &
     3010                                      diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
     3011                               )
     3012             ELSE
     3013                k_cr = k_cr0
     3014             ENDIF
     3015!
     3016!--          Accretion rate (Seifert and Beheng, 2006):
     3017             accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *                      &
     3018                    SQRT( rho_surface * hyrho(k) )
     3019             accr = MIN( accr, qc(k,j,i) / dt_micro )
     3020
     3021             qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
     3022             qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
     3023             IF ( microphysics_morrison )  THEN
     3024                nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc *               &
     3025                                             hyrho(k) * dt_micro * flag        &
     3026                                         )
     3027             ENDIF
     3028
     3029
     3030          ENDIF
     3031
     3032       ENDDO
     3033
     3034    END SUBROUTINE accretion_ij
     3035
    25663036
    25673037!------------------------------------------------------------------------------!
     
    26223092
    26233093    END SUBROUTINE selfcollection_breakup
     3094
     3095
     3096!------------------------------------------------------------------------------!
     3097! Description:
     3098! ------------
     3099!> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
     3100!------------------------------------------------------------------------------!
     3101    SUBROUTINE selfcollection_breakup_ij( i, j )
     3102
     3103       IMPLICIT NONE
     3104
     3105       INTEGER(iwp) ::  i                 !<
     3106       INTEGER(iwp) ::  j                 !<
     3107       INTEGER(iwp) ::  k                 !<
     3108
     3109       REAL(wp)     ::  breakup           !<
     3110       REAL(wp)     ::  dr                !<
     3111       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
     3112       REAL(wp)     ::  phi_br            !<
     3113       REAL(wp)     ::  selfcoll          !<
     3114
     3115       DO  k = nzb+1, nzt
     3116!
     3117!--       Predetermine flag to mask topography
     3118          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     3119
     3120          IF ( qr(k,j,i) > eps_sb )  THEN
     3121!
     3122!--          Selfcollection rate (Seifert and Beheng, 2001):
     3123             selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * SQRT( hyrho(k) * rho_surface )
     3124!
     3125!--          Weight averaged diameter of rain drops:
     3126             dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
     3127!
     3128!--          Collisional breakup rate (Seifert, 2008):
     3129             IF ( dr >= 0.3E-3_wp )  THEN
     3130                phi_br  = k_br * ( dr - 1.1E-3_wp )
     3131                breakup = selfcoll * ( phi_br + 1.0_wp )
     3132             ELSE
     3133                breakup = 0.0_wp
     3134             ENDIF
     3135
     3136             selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
     3137             nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
     3138
     3139          ENDIF
     3140       ENDDO
     3141
     3142    END SUBROUTINE selfcollection_breakup_ij
    26243143
    26253144
     
    27533272! Description:
    27543273! ------------
    2755 !> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
    2756 !------------------------------------------------------------------------------!
    2757     SUBROUTINE sedimentation_cloud
    2758 
    2759 
    2760        IMPLICIT NONE
    2761 
    2762        INTEGER(iwp) ::  i             !<
    2763        INTEGER(iwp) ::  j             !<
    2764        INTEGER(iwp) ::  k             !<
    2765 
    2766        REAL(wp) ::  flag              !< flag to mask topography grid points
    2767        REAL(wp) ::  nc_sedi           !<
    2768 
    2769        REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
    2770        REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc !<
    2771 
    2772 
    2773        CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
    2774 
    2775        sed_qc(nzt+1) = 0.0_wp
    2776        sed_nc(nzt+1) = 0.0_wp
    2777 
    2778        DO  i = nxlg, nxrg
    2779           DO  j = nysg, nyng
    2780              DO  k = nzt, nzb+1, -1
    2781 !
    2782 !--             Predetermine flag to mask topography
    2783                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    2784 
    2785                 IF ( microphysics_morrison ) THEN
    2786                    nc_sedi = nc(k,j,i)
    2787                 ELSE
    2788                    nc_sedi = nc_const
    2789                 ENDIF
    2790 
    2791 !
    2792 !--             Sedimentation fluxes for number concentration are only calculated
    2793 !--             for cloud_scheme = 'morrison'
    2794                 IF ( microphysics_morrison ) THEN
    2795                    IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
    2796                       sed_nc(k) = sed_qc_const *                               &
    2797                               ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
    2798                               ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
    2799                    ELSE
    2800                       sed_nc(k) = 0.0_wp
    2801                    ENDIF
    2802 
    2803                    sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *           &
    2804                                     nc(k,j,i) / dt_micro + sed_nc(k+1)         &
    2805                                   ) * flag
    2806 
    2807                    nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *       &
    2808                                          ddzu(k+1) / hyrho(k) * dt_micro * flag
    2809                 ENDIF
    2810 
    2811                 IF ( qc(k,j,i) > eps_sb .AND.  nc_sedi > eps_mr )  THEN
    2812                    sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *  &
    2813                                ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
    2814                                                                            flag
    2815                 ELSE
    2816                    sed_qc(k) = 0.0_wp
    2817                 ENDIF
    2818 
    2819                 sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
    2820                                             dt_micro + sed_qc(k+1)             &
    2821                                ) * flag
    2822 
    2823                 q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
    2824                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
    2825                 qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
    2826                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
    2827                 pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
    2828                                         ddzu(k+1) / hyrho(k) * lv_d_cp *       &
    2829                                         d_exner(k) * dt_micro            * flag
    2830 
    2831 !
    2832 !--             Compute the precipitation rate due to cloud (fog) droplets
    2833                 IF ( call_microphysics_at_all_substeps )  THEN
    2834                    prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
    2835                                 * weight_substep(intermediate_timestep_count)  &
    2836                                 * flag
    2837                 ELSE
    2838                    prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
    2839                 ENDIF
    2840 
    2841              ENDDO
    2842           ENDDO
    2843        ENDDO
    2844 
    2845        CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
    2846 
    2847     END SUBROUTINE sedimentation_cloud
    2848 
    2849 
    2850 !------------------------------------------------------------------------------!
    2851 ! Description:
    2852 ! ------------
    2853 !> Computation of sedimentation flux. Implementation according to Stevens
    2854 !> and Seifert (2008). Code is based on UCLA-LES.
    2855 !------------------------------------------------------------------------------!
    2856     SUBROUTINE sedimentation_rain
    2857 
    2858        IMPLICIT NONE
    2859 
    2860        INTEGER(iwp) ::  i             !< running index x direction
    2861        INTEGER(iwp) ::  j             !< running index y direction
    2862        INTEGER(iwp) ::  k             !< running index z direction
    2863        INTEGER(iwp) ::  k_run         !<
    2864        INTEGER(iwp) ::  m             !< running index surface elements
    2865        INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
    2866        INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
    2867 
    2868        REAL(wp)     ::  c_run                      !<
    2869        REAL(wp)     ::  d_max                      !<
    2870        REAL(wp)     ::  d_mean                     !<
    2871        REAL(wp)     ::  d_min                      !<
    2872        REAL(wp)     ::  dr                         !<
    2873        REAL(wp)     ::  flux                       !<
    2874        REAL(wp)     ::  flag                       !< flag to mask topography grid points
    2875        REAL(wp)     ::  lambda_r                   !<
    2876        REAL(wp)     ::  mu_r                       !<
    2877        REAL(wp)     ::  z_run                      !<
    2878 
    2879        REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
    2880        REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
    2881        REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
    2882        REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
    2883        REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
    2884        REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
    2885        REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
    2886        REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
    2887 
    2888        CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
    2889 
    2890 !
    2891 !--    Compute velocities
    2892        DO  i = nxlg, nxrg
    2893           DO  j = nysg, nyng
    2894              DO  k = nzb+1, nzt
    2895 !
    2896 !--             Predetermine flag to mask topography
    2897                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    2898 
    2899                 IF ( qr(k,j,i) > eps_sb )  THEN
    2900 !
    2901 !--                Weight averaged diameter of rain drops:
    2902                    dr = ( hyrho(k) * qr(k,j,i) /                               &
    2903                           nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
    2904 !
    2905 !--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
    2906 !--                Stevens and Seifert, 2008):
    2907                    mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
    2908                                                      ( dr - 1.4E-3_wp ) ) )
    2909 !
    2910 !--                Slope parameter of gamma distribution (Seifert, 2008):
    2911                    lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
    2912                                 ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
    2913 
    2914                    w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
    2915                                                a_term - b_term * ( 1.0_wp +    &
    2916                                                   c_term /                     &
    2917                                                   lambda_r )**( -1.0_wp *      &
    2918                                                   ( mu_r + 1.0_wp ) )          &
    2919                                               )                                &
    2920                                 ) * flag
    2921 
    2922                    w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
    2923                                                a_term - b_term * ( 1.0_wp +    &
    2924                                                   c_term /                     &
    2925                                                   lambda_r )**( -1.0_wp *      &
    2926                                                   ( mu_r + 4.0_wp ) )          &
    2927                                              )                                 &
    2928                                 ) * flag
    2929                 ELSE
    2930                    w_nr(k) = 0.0_wp
    2931                    w_qr(k) = 0.0_wp
    2932                 ENDIF
    2933              ENDDO
    2934 !
    2935 !--          Adjust boundary values using surface data type.
    2936 !--          Upward-facing
    2937              surf_s = bc_h(0)%start_index(j,i)
    2938              surf_e = bc_h(0)%end_index(j,i)
    2939              DO  m = surf_s, surf_e
    2940                 k         = bc_h(0)%k(m)
    2941                 w_nr(k-1) = w_nr(k)
    2942                 w_qr(k-1) = w_qr(k)
    2943              ENDDO
    2944 !
    2945 !--          Downward-facing
    2946              surf_s = bc_h(1)%start_index(j,i)
    2947              surf_e = bc_h(1)%end_index(j,i)
    2948              DO  m = surf_s, surf_e
    2949                 k         = bc_h(1)%k(m)
    2950                 w_nr(k+1) = w_nr(k)
    2951                 w_qr(k+1) = w_qr(k)
    2952              ENDDO
    2953 !
    2954 !--          Model top boundary value
    2955              w_nr(nzt+1) = 0.0_wp
    2956              w_qr(nzt+1) = 0.0_wp
    2957 !
    2958 !--          Compute Courant number
    2959              DO  k = nzb+1, nzt
    2960 !
    2961 !--             Predetermine flag to mask topography
    2962                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    2963 
    2964                 c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
    2965                                       2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
    2966                           dt_micro * ddzu(k) * flag
    2967                 c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
    2968                                       2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
    2969                           dt_micro * ddzu(k) * flag
    2970              ENDDO
    2971 !
    2972 !--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
    2973              IF ( limiter_sedimentation )  THEN
    2974 
    2975                 DO k = nzb+1, nzt
    2976 !
    2977 !--                Predetermine flag to mask topography
    2978                    flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    2979 
    2980                    d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
    2981                    d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
    2982                    d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
    2983 
    2984                    qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
    2985                                                               2.0_wp * d_max,  &
    2986                                                               ABS( d_mean ) )  &
    2987                                                       * flag
    2988 
    2989                    d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
    2990                    d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
    2991                    d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
    2992 
    2993                    nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
    2994                                                               2.0_wp * d_max,  &
    2995                                                               ABS( d_mean ) )
    2996                 ENDDO
    2997 
    2998              ELSE
    2999 
    3000                 nr_slope = 0.0_wp
    3001                 qr_slope = 0.0_wp
    3002 
    3003              ENDIF
    3004 
    3005              sed_nr(nzt+1) = 0.0_wp
    3006              sed_qr(nzt+1) = 0.0_wp
    3007 !
    3008 !--          Compute sedimentation flux
    3009              DO  k = nzt, nzb+1, -1
    3010 !
    3011 !--             Predetermine flag to mask topography
    3012                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    3013 !
    3014 !--             Sum up all rain drop number densities which contribute to the flux
    3015 !--             through k-1/2
    3016                 flux  = 0.0_wp
    3017                 z_run = 0.0_wp ! height above z(k)
    3018                 k_run = k
    3019                 c_run = MIN( 1.0_wp, c_nr(k) )
    3020                 DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
    3021                    flux  = flux + hyrho(k_run) *                               &
    3022                            ( nr(k_run,j,i) + nr_slope(k_run) *                 &
    3023                            ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)  &
    3024                                               * flag
    3025                    z_run = z_run + dzu(k_run) * flag
    3026                    k_run = k_run + 1          * flag
    3027                    c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )    &
    3028                                               * flag
    3029                 ENDDO
    3030 !
    3031 !--             It is not allowed to sediment more rain drop number density than
    3032 !--             available
    3033                 flux = MIN( flux,                                              &
    3034                             hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
    3035                             dt_micro                                           &
    3036                           )
    3037 
    3038                 sed_nr(k) = flux / dt_micro * flag
    3039                 nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
    3040                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
    3041 !
    3042 !--             Sum up all rain water content which contributes to the flux
    3043 !--             through k-1/2
    3044                 flux  = 0.0_wp
    3045                 z_run = 0.0_wp ! height above z(k)
    3046                 k_run = k
    3047                 c_run = MIN( 1.0_wp, c_qr(k) )
    3048 
    3049                 DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
    3050 
    3051                    flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
    3052                                   qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
    3053                                   0.5_wp ) * c_run * dzu(k_run) * flag
    3054                    z_run = z_run + dzu(k_run)                   * flag
    3055                    k_run = k_run + 1                            * flag
    3056                    c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )    &
    3057                                                                 * flag
    3058 
    3059                 ENDDO
    3060 !
    3061 !--             It is not allowed to sediment more rain water content than
    3062 !--             available
    3063                 flux = MIN( flux,                                              &
    3064                             hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
    3065                             dt_micro                                           &
    3066                           )
    3067 
    3068                 sed_qr(k) = flux / dt_micro * flag
    3069 
    3070                 qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
    3071                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
    3072                 q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
    3073                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
    3074                 pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
    3075                                         ddzu(k+1) / hyrho(k) * lv_d_cp *       &
    3076                                         d_exner(k) * dt_micro            * flag
    3077 !
    3078 !--             Compute the rain rate
    3079                 IF ( call_microphysics_at_all_substeps )  THEN
    3080                    prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
    3081                                 * weight_substep(intermediate_timestep_count) &
    3082                                 * flag
    3083                 ELSE
    3084                    prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
    3085                 ENDIF
    3086 
    3087              ENDDO
    3088           ENDDO
    3089        ENDDO
    3090 
    3091        CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
    3092 
    3093     END SUBROUTINE sedimentation_rain
    3094 
    3095 
    3096 !------------------------------------------------------------------------------!
    3097 ! Description:
    3098 ! ------------
    3099 !> Computation of the precipitation amount due to gravitational settling of
    3100 !> rain and cloud (fog) droplets
    3101 !------------------------------------------------------------------------------!
    3102     SUBROUTINE calc_precipitation_amount
    3103 
    3104        IMPLICIT NONE
    3105 
    3106        INTEGER(iwp) ::  i             !< running index x direction
    3107        INTEGER(iwp) ::  j             !< running index y direction
    3108        INTEGER(iwp) ::  k             !< running index y direction
    3109        INTEGER(iwp) ::  m             !< running index surface elements
    3110 
    3111        IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
    3112             ( .NOT. call_microphysics_at_all_substeps .OR.                     &
    3113             intermediate_timestep_count == intermediate_timestep_count_max ) ) &
    3114        THEN
    3115 !
    3116 !--       Run over all upward-facing surface elements, i.e. non-natural,
    3117 !--       natural and urban
    3118           DO  m = 1, bc_h(0)%ns
    3119              i = bc_h(0)%i(m)
    3120              j = bc_h(0)%j(m)
    3121              k = bc_h(0)%k(m)
    3122              precipitation_amount(j,i) = precipitation_amount(j,i) +           &
    3123                                                prr(k,j,i) * hyrho(k) * dt_3d
    3124           ENDDO
    3125 
    3126        ENDIF
    3127 
    3128     END SUBROUTINE calc_precipitation_amount
    3129 
    3130 
    3131 !------------------------------------------------------------------------------!
    3132 ! Description:
    3133 ! ------------
    3134 !> Control of microphysics for grid points i,j
    3135 !------------------------------------------------------------------------------!
    3136 
    3137     SUBROUTINE bcm_actions_ij( i, j )
    3138 
    3139        IMPLICIT NONE
    3140 
    3141        INTEGER(iwp) ::  i                 !<
    3142        INTEGER(iwp) ::  j                 !<
    3143 
    3144        IF ( large_scale_forcing  .AND.  lsf_surf ) THEN
    3145 !
    3146 !--       Calculate vertical profile of the hydrostatic pressure (hyp)
    3147           hyp    = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)
    3148           d_exner = exner_function_invers(hyp)
    3149           exner = 1.0_wp / exner_function_invers(hyp)
    3150           hyrho  = ideal_gas_law_rho_pt(hyp, pt_init)
    3151 !
    3152 !--       Compute reference density
    3153           rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))
    3154        ENDIF
    3155 
    3156 !
    3157 !--    Compute length of time step
    3158        IF ( call_microphysics_at_all_substeps )  THEN
    3159           dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
    3160        ELSE
    3161           dt_micro = dt_3d
    3162        ENDIF
    3163 !
    3164 !--    Reset precipitation rate
    3165        IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0_wp
    3166 
    3167 !
    3168 !--    Compute cloud physics
    3169        IF( microphysics_kessler )  THEN
    3170 
    3171           CALL autoconversion_kessler( i,j )
    3172           IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
    3173 
    3174        ELSEIF ( microphysics_seifert )  THEN
    3175 
    3176           CALL adjust_cloud( i,j )
    3177           IF ( microphysics_morrison ) CALL activation( i,j )
    3178           IF ( microphysics_morrison ) CALL condensation( i,j )
    3179           CALL autoconversion( i,j )
    3180           CALL accretion( i,j )
    3181           CALL selfcollection_breakup( i,j )
    3182           CALL evaporation_rain( i,j )
    3183           CALL sedimentation_rain( i,j )
    3184           IF ( cloud_water_sedimentation )  CALL sedimentation_cloud( i,j )
    3185 
    3186        ENDIF
    3187 
    3188        CALL calc_precipitation_amount( i,j )
    3189 
    3190     END SUBROUTINE bcm_actions_ij
    3191 
    3192 !------------------------------------------------------------------------------!
    3193 ! Description:
    3194 ! ------------
    3195 !> Adjust number of raindrops to avoid nonlinear effects in
    3196 !> sedimentation and evaporation of rain drops due to too small or
    3197 !> too big weights of rain drops (Stevens and Seifert, 2008).
    3198 !> The same procedure is applied to cloud droplets if they are determined
    3199 !> prognostically. Call for grid point i,j
    3200 !------------------------------------------------------------------------------!
    3201     SUBROUTINE adjust_cloud_ij( i, j )
    3202 
    3203        IMPLICIT NONE
    3204 
    3205        INTEGER(iwp) ::  i                 !<
    3206        INTEGER(iwp) ::  j                 !<
    3207        INTEGER(iwp) ::  k                 !<
    3208 
    3209        REAL(wp) ::  flag                  !< flag to indicate first grid level above surface
    3210 
    3211        DO  k = nzb+1, nzt
    3212 !
    3213 !--       Predetermine flag to mask topography
    3214           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    3215 
    3216           IF ( qr(k,j,i) <= eps_sb )  THEN
    3217              qr(k,j,i) = 0.0_wp
    3218              nr(k,j,i) = 0.0_wp
    3219           ELSE
    3220 !
    3221 !--          Adjust number of raindrops to avoid nonlinear effects in
    3222 !--          sedimentation and evaporation of rain drops due to too small or
    3223 !--          too big weights of rain drops (Stevens and Seifert, 2008).
    3224              IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
    3225                 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag
    3226              ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
    3227                 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag
    3228              ENDIF
    3229 
    3230           ENDIF
    3231 
    3232           IF ( microphysics_morrison ) THEN
    3233              IF ( qc(k,j,i) <= eps_sb )  THEN
    3234                 qc(k,j,i) = 0.0_wp
    3235                 nc(k,j,i) = 0.0_wp
    3236              ELSE
    3237                 IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
    3238                    nc(k,j,i) = qc(k,j,i) * hyrho(k) / xamin * flag
    3239                 ENDIF
    3240              ENDIF
    3241           ENDIF
    3242 
    3243        ENDDO
    3244 
    3245     END SUBROUTINE adjust_cloud_ij
    3246 
    3247 !------------------------------------------------------------------------------!
    3248 ! Description:
    3249 ! ------------
    3250 !> Calculate number of activated condensation nucleii after simple activation
    3251 !> scheme of Twomey, 1959.
    3252 !------------------------------------------------------------------------------!
    3253     SUBROUTINE activation_ij( i, j )
    3254 
    3255        IMPLICIT NONE
    3256 
    3257        INTEGER(iwp) ::  i                 !<
    3258        INTEGER(iwp) ::  j                 !<
    3259        INTEGER(iwp) ::  k                 !<
    3260 
    3261        REAL(wp)     ::  activ             !<
    3262        REAL(wp)     ::  afactor           !<
    3263        REAL(wp)     ::  beta_act          !<
    3264        REAL(wp)     ::  bfactor           !<
    3265        REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
    3266        REAL(wp)     ::  k_act             !<
    3267        REAL(wp)     ::  n_act             !<
    3268        REAL(wp)     ::  n_ccn             !<
    3269        REAL(wp)     ::  s_0               !<
    3270        REAL(wp)     ::  sat_max           !<
    3271        REAL(wp)     ::  sigma             !<
    3272        REAL(wp)     ::  sigma_act         !<
    3273 
    3274        DO  k = nzb+1, nzt
    3275 !
    3276 !--       Predetermine flag to mask topography
    3277           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    3278 !
    3279 !--       Call calculation of supersaturation 
    3280           CALL supersaturation ( i, j, k )
    3281 !
    3282 !--       Prescribe parameters for activation
    3283 !--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
    3284           k_act  = 0.7_wp
    3285           activ  = 0.0_wp
    3286 
    3287           IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
    3288 !
    3289 !--          Compute the number of activated Aerosols
    3290 !--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
    3291              n_act     = na_init * sat**k_act
    3292 !
    3293 !--          Compute the number of cloud droplets
    3294 !--          (see: Morrison + Grabowski, 2007, JAS, 64)
    3295 !            activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
    3296 
    3297 !
    3298 !--          Compute activation rate after Khairoutdinov and Kogan
    3299 !--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)
    3300              sat_max = 0.8_wp / 100.0_wp
    3301              activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN             &
    3302                  ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /         &
    3303                   dt_micro
    3304 
    3305              nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
    3306           ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk )  THEN
    3307 !
    3308 !--          Curvature effect (afactor) with surface tension
    3309 !--          parameterization by Straka (2009)
    3310              sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )
    3311              afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )
    3312 !
    3313 !--          Solute effect (bfactor)
    3314              bfactor = vanthoff * molecular_weight_of_water *                  &
    3315                   rho_s / ( molecular_weight_of_solute * rho_l )
    3316 
    3317 !
    3318 !--          Prescribe power index that describes the soluble fraction
    3319 !--          of an aerosol particle (beta).
    3320 !--          (see: Morrison + Grabowski, 2007, JAS, 64)
    3321              beta_act  = 0.5_wp
    3322              sigma_act = sigma_bulk**( 1.0_wp + beta_act )
    3323 !
    3324 !--          Calculate mean geometric supersaturation (s_0) with
    3325 !--          parameterization by Khvorostyanov and Curry (2006)
    3326              s_0   = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) *         &
    3327                ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
    3328 
    3329 !
    3330 !--          Calculate number of activated CCN as a function of
    3331 !--          supersaturation and taking Koehler theory into account
    3332 !--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
    3333              n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                    &
    3334                 LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )
    3335              activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
    3336 
    3337              nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)
    3338           ENDIF
    3339 
    3340        ENDDO
    3341 
    3342     END SUBROUTINE activation_ij
    3343 
    3344 !------------------------------------------------------------------------------!
    3345 ! Description:
    3346 ! ------------
    3347 !> Calculate condensation rate for cloud water content (after Khairoutdinov and
    3348 !> Kogan, 2000).
    3349 !------------------------------------------------------------------------------!
    3350     SUBROUTINE condensation_ij( i, j )
    3351 
    3352        IMPLICIT NONE
    3353 
    3354        INTEGER(iwp) ::  i                 !<
    3355        INTEGER(iwp) ::  j                 !<
    3356        INTEGER(iwp) ::  k                 !<
    3357 
    3358        REAL(wp)     ::  cond              !<
    3359        REAL(wp)     ::  cond_max          !<
    3360        REAL(wp)     ::  dc                !<
    3361        REAL(wp)     ::  evap              !<
    3362        REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
    3363        REAL(wp)     ::  g_fac             !<
    3364        REAL(wp)     ::  nc_0              !<
    3365        REAL(wp)     ::  temp              !<
    3366        REAL(wp)     ::  xc                !<
    3367 
    3368 
    3369        DO  k = nzb+1, nzt
    3370 !
    3371 !--       Predetermine flag to mask topography
    3372           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    3373 !
    3374 !--       Call calculation of supersaturation 
    3375           CALL supersaturation ( i, j, k )
    3376 !
    3377 !--       Actual temperature:
    3378           temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )
    3379 
    3380           g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
    3381                               l_v / ( thermal_conductivity_l * temp )    &
    3382                               + r_v * temp / ( diff_coeff_l * e_s )      &
    3383                             )
    3384 !
    3385 !--       Mean weight of cloud drops
    3386           IF ( nc(k,j,i) <= 0.0_wp) CYCLE
    3387           xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)
    3388 !
    3389 !--       Weight averaged diameter of cloud drops:
    3390           dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
    3391 !
    3392 !--       Integral diameter of cloud drops
    3393           nc_0 = nc(k,j,i) * dc
    3394 !
    3395 !--       Condensation needs only to be calculated in supersaturated regions
    3396           IF ( sat > 0.0_wp )  THEN
    3397 !
    3398 !--          Condensation rate of cloud water content
    3399 !--          after KK scheme.
    3400 !--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
    3401              cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
    3402              cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
    3403              cond      = MIN( cond, cond_max / dt_micro )
    3404 
    3405              qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag
    3406           ELSEIF ( sat < 0.0_wp ) THEN
    3407              evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
    3408              evap      = MAX( evap, -qc(k,j,i) / dt_micro )
    3409 
    3410              qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag
    3411           ENDIF
    3412        ENDDO
    3413 
    3414     END SUBROUTINE condensation_ij
    3415 
    3416 
    3417 !------------------------------------------------------------------------------!
    3418 ! Description:
    3419 ! ------------
    3420 !> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j
    3421 !------------------------------------------------------------------------------!
    3422     SUBROUTINE autoconversion_ij( i, j )
    3423 
    3424        IMPLICIT NONE
    3425 
    3426        INTEGER(iwp) ::  i                 !<
    3427        INTEGER(iwp) ::  j                 !<
    3428        INTEGER(iwp) ::  k                 !<
    3429 
    3430        REAL(wp)     ::  alpha_cc          !<
    3431        REAL(wp)     ::  autocon           !<
    3432        REAL(wp)     ::  dissipation       !<
    3433        REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
    3434        REAL(wp)     ::  k_au              !<
    3435        REAL(wp)     ::  l_mix             !<
    3436        REAL(wp)     ::  nc_auto           !<
    3437        REAL(wp)     ::  nu_c              !<
    3438        REAL(wp)     ::  phi_au            !<
    3439        REAL(wp)     ::  r_cc              !<
    3440        REAL(wp)     ::  rc                !<
    3441        REAL(wp)     ::  re_lambda         !<
    3442        REAL(wp)     ::  sigma_cc          !<
    3443        REAL(wp)     ::  tau_cloud         !<
    3444        REAL(wp)     ::  xc                !<
    3445 
    3446        DO  k = nzb+1, nzt
    3447 !
    3448 !--       Predetermine flag to mask topography
    3449           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    3450           IF ( microphysics_morrison ) THEN
    3451              nc_auto = nc(k,j,i)
    3452           ELSE
    3453              nc_auto = nc_const
    3454           ENDIF
    3455 
    3456           IF ( qc(k,j,i) > eps_sb  .AND.  nc_auto > eps_mr )  THEN
    3457 
    3458              k_au = k_cc / ( 20.0_wp * x0 )
    3459 !
    3460 !--          Intern time scale of coagulation (Seifert and Beheng, 2006):
    3461 !--          (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))
    3462              tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ),     &
    3463                               0.0_wp )
    3464 !
    3465 !--          Universal function for autoconversion process
    3466 !--          (Seifert and Beheng, 2006):
    3467              phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
    3468 !
    3469 !--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
    3470 !--          (Use constant nu_c = 1.0_wp instead?)
    3471              nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )
    3472 !
    3473 !--          Mean weight of cloud droplets:
    3474              xc = hyrho(k) * qc(k,j,i) / nc_auto
    3475 !
    3476 !--          Parameterized turbulence effects on autoconversion (Seifert,
    3477 !--          Nuijens and Stevens, 2010)
    3478              IF ( collision_turbulence )  THEN
    3479 !
    3480 !--             Weight averaged radius of cloud droplets:
    3481                 rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
    3482 
    3483                 alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )
    3484                 r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )
    3485                 sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
    3486 !
    3487 !--             Mixing length (neglecting distance to ground and stratification)
    3488                 l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
    3489 !
    3490 !--             Limit dissipation rate according to Seifert, Nuijens and
    3491 !--             Stevens (2010)
    3492                 dissipation = MIN( 0.06_wp, diss(k,j,i) )
    3493 !
    3494 !--             Compute Taylor-microscale Reynolds number:
    3495                 re_lambda = 6.0_wp / 11.0_wp *                                 &
    3496                             ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *         &
    3497                             SQRT( 15.0_wp / kin_vis_air ) *                    &
    3498                             dissipation**( 1.0_wp / 6.0_wp )
    3499 !
    3500 !--             The factor of 1.0E4 is needed to convert the dissipation rate
    3501 !--             from m2 s-3 to cm2 s-3.
    3502                 k_au = k_au * ( 1.0_wp +                                       &
    3503                                 dissipation * 1.0E4_wp *                       &
    3504                                 ( re_lambda * 1.0E-3_wp )**0.25_wp *           &
    3505                                 ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) /  &
    3506                                                   sigma_cc )**2                &
    3507                                                 ) + beta_cc                    &
    3508                                 )                                              &
    3509                               )
    3510              ENDIF
    3511 !
    3512 !--          Autoconversion rate (Seifert and Beheng, 2006):
    3513              autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) /          &
    3514                        ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 *            &
    3515                        ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) *       &
    3516                        rho_surface
    3517              autocon = MIN( autocon, qc(k,j,i) / dt_micro )
    3518 
    3519              qr(k,j,i) = qr(k,j,i) + autocon * dt_micro                 * flag
    3520              qc(k,j,i) = qc(k,j,i) - autocon * dt_micro                 * flag
    3521              nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag
    3522              IF ( microphysics_morrison ) THEN
    3523                 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *                &
    3524                                     autocon / x0 * hyrho(k) * dt_micro * flag )
    3525              ENDIF
    3526 
    3527           ENDIF
    3528 
    3529        ENDDO
    3530 
    3531     END SUBROUTINE autoconversion_ij
    3532 
    3533 !------------------------------------------------------------------------------!
    3534 ! Description:
    3535 ! ------------
    3536 !> Autoconversion process (Kessler, 1969).
    3537 !------------------------------------------------------------------------------!
    3538     SUBROUTINE autoconversion_kessler_ij( i, j )
    3539 
    3540 
    3541        IMPLICIT NONE
    3542 
    3543        INTEGER(iwp) ::  i      !<
    3544        INTEGER(iwp) ::  j      !<
    3545        INTEGER(iwp) ::  k      !<
    3546        INTEGER(iwp) ::  k_wall !< topography top index
    3547 
    3548        REAL(wp)    ::  dqdt_precip !<
    3549        REAL(wp)    ::  flag              !< flag to indicate first grid level above surface
    3550 
    3551 !
    3552 !--    Determine vertical index of topography top
    3553        k_wall = get_topography_top_index_ji( j, i, 's' )
    3554        DO  k = nzb+1, nzt
    3555 !
    3556 !--       Predetermine flag to mask topography
    3557           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    3558 
    3559           IF ( qc(k,j,i) > ql_crit )  THEN
    3560              dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )
    3561           ELSE
    3562              dqdt_precip = 0.0_wp
    3563           ENDIF
    3564 
    3565           qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag
    3566           q(k,j,i)  = q(k,j,i)  - dqdt_precip * dt_micro * flag
    3567           pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k)  &
    3568                                                                  * flag
    3569 
    3570 !
    3571 !--       Compute the rain rate (stored on surface grid point)
    3572           prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag
    3573 
    3574        ENDDO
    3575 
    3576     END SUBROUTINE autoconversion_kessler_ij
    3577 
    3578 !------------------------------------------------------------------------------!
    3579 ! Description:
    3580 ! ------------
    3581 !> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j
    3582 !------------------------------------------------------------------------------!
    3583     SUBROUTINE accretion_ij( i, j )
    3584 
    3585        IMPLICIT NONE
    3586 
    3587        INTEGER(iwp) ::  i                 !<
    3588        INTEGER(iwp) ::  j                 !<
    3589        INTEGER(iwp) ::  k                 !<
    3590 
    3591        REAL(wp)     ::  accr              !<
    3592        REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
    3593        REAL(wp)     ::  k_cr              !<
    3594        REAL(wp)     ::  nc_accr           !<
    3595        REAL(wp)     ::  phi_ac            !<
    3596        REAL(wp)     ::  tau_cloud         !<
    3597        REAL(wp)     ::  xc                !<
    3598 
    3599 
    3600        DO  k = nzb+1, nzt
    3601 !
    3602 !--       Predetermine flag to mask topography
    3603           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    3604           IF ( microphysics_morrison ) THEN
    3605              nc_accr = nc(k,j,i)
    3606           ELSE
    3607              nc_accr = nc_const
    3608           ENDIF
    3609 
    3610           IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )  .AND.  &
    3611                ( nc_accr > eps_mr ) )  THEN
    3612 !
    3613 !--          Intern time scale of coagulation (Seifert and Beheng, 2006):
    3614              tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
    3615 !
    3616 !--          Universal function for accretion process
    3617 !--          (Seifert and Beheng, 2001):
    3618              phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
    3619 
    3620 !
    3621 !--          Mean weight of cloud drops
    3622              xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)
    3623 !
    3624 !--          Parameterized turbulence effects on autoconversion (Seifert,
    3625 !--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
    3626 !--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
    3627              IF ( collision_turbulence )  THEN
    3628                 k_cr = k_cr0 * ( 1.0_wp + 0.05_wp *                            &
    3629                                  MIN( 600.0_wp,                                &
    3630                                       diss(k,j,i) * 1.0E4_wp )**0.25_wp        &
    3631                                )
    3632              ELSE
    3633                 k_cr = k_cr0
    3634              ENDIF
    3635 !
    3636 !--          Accretion rate (Seifert and Beheng, 2006):
    3637              accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac *                      &
    3638                     SQRT( rho_surface * hyrho(k) )
    3639              accr = MIN( accr, qc(k,j,i) / dt_micro )
    3640 
    3641              qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
    3642              qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
    3643              IF ( microphysics_morrison )  THEN
    3644                 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc *               &
    3645                                              hyrho(k) * dt_micro * flag        &
    3646                                          )
    3647              ENDIF
    3648 
    3649 
    3650           ENDIF
    3651 
    3652        ENDDO
    3653 
    3654     END SUBROUTINE accretion_ij
    3655 
    3656 
    3657 !------------------------------------------------------------------------------!
    3658 ! Description:
    3659 ! ------------
    3660 !> Collisional breakup rate (Seifert, 2008). Call for grid point i,j
    3661 !------------------------------------------------------------------------------!
    3662     SUBROUTINE selfcollection_breakup_ij( i, j )
    3663 
    3664        IMPLICIT NONE
    3665 
    3666        INTEGER(iwp) ::  i                 !<
    3667        INTEGER(iwp) ::  j                 !<
    3668        INTEGER(iwp) ::  k                 !<
    3669 
    3670        REAL(wp)     ::  breakup           !<
    3671        REAL(wp)     ::  dr                !<
    3672        REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
    3673        REAL(wp)     ::  phi_br            !<
    3674        REAL(wp)     ::  selfcoll          !<
    3675 
    3676        DO  k = nzb+1, nzt
    3677 !
    3678 !--       Predetermine flag to mask topography
    3679           flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    3680 
    3681           IF ( qr(k,j,i) > eps_sb )  THEN
    3682 !
    3683 !--          Selfcollection rate (Seifert and Beheng, 2001):
    3684              selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * SQRT( hyrho(k) * rho_surface )
    3685 !
    3686 !--          Weight averaged diameter of rain drops:
    3687              dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
    3688 !
    3689 !--          Collisional breakup rate (Seifert, 2008):
    3690              IF ( dr >= 0.3E-3_wp )  THEN
    3691                 phi_br  = k_br * ( dr - 1.1E-3_wp )
    3692                 breakup = selfcoll * ( phi_br + 1.0_wp )
    3693              ELSE
    3694                 breakup = 0.0_wp
    3695              ENDIF
    3696 
    3697              selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )
    3698              nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag
    3699 
    3700           ENDIF
    3701        ENDDO
    3702 
    3703     END SUBROUTINE selfcollection_breakup_ij
    3704 
    3705 
    3706 !------------------------------------------------------------------------------!
    3707 ! Description:
    3708 ! ------------
    37093274!> Evaporation of precipitable water. Condensation is neglected for
    37103275!> precipitable water. Call for grid point i,j
     
    37393304          IF ( qr(k,j,i) > eps_sb )  THEN
    37403305!
    3741 !--          Call calculation of supersaturation 
     3306!--          Call calculation of supersaturation
    37423307             CALL supersaturation ( i, j, k )
    37433308!
     
    38223387! ------------
    38233388!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
     3389!------------------------------------------------------------------------------!
     3390    SUBROUTINE sedimentation_cloud
     3391
     3392
     3393       IMPLICIT NONE
     3394
     3395       INTEGER(iwp) ::  i             !<
     3396       INTEGER(iwp) ::  j             !<
     3397       INTEGER(iwp) ::  k             !<
     3398
     3399       REAL(wp) ::  flag              !< flag to mask topography grid points
     3400       REAL(wp) ::  nc_sedi           !<
     3401
     3402       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
     3403       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc !<
     3404
     3405
     3406       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
     3407
     3408       sed_qc(nzt+1) = 0.0_wp
     3409       sed_nc(nzt+1) = 0.0_wp
     3410
     3411       DO  i = nxlg, nxrg
     3412          DO  j = nysg, nyng
     3413             DO  k = nzt, nzb+1, -1
     3414!
     3415!--             Predetermine flag to mask topography
     3416                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     3417
     3418                IF ( microphysics_morrison ) THEN
     3419                   nc_sedi = nc(k,j,i)
     3420                ELSE
     3421                   nc_sedi = nc_const
     3422                ENDIF
     3423
     3424!
     3425!--             Sedimentation fluxes for number concentration are only calculated
     3426!--             for cloud_scheme = 'morrison'
     3427                IF ( microphysics_morrison ) THEN
     3428                   IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
     3429                      sed_nc(k) = sed_qc_const *                               &
     3430                              ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
     3431                              ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
     3432                   ELSE
     3433                      sed_nc(k) = 0.0_wp
     3434                   ENDIF
     3435
     3436                   sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *           &
     3437                                    nc(k,j,i) / dt_micro + sed_nc(k+1)         &
     3438                                  ) * flag
     3439
     3440                   nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *       &
     3441                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
     3442                ENDIF
     3443
     3444                IF ( qc(k,j,i) > eps_sb .AND.  nc_sedi > eps_mr )  THEN
     3445                   sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *  &
     3446                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
     3447                                                                           flag
     3448                ELSE
     3449                   sed_qc(k) = 0.0_wp
     3450                ENDIF
     3451
     3452                sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) /   &
     3453                                            dt_micro + sed_qc(k+1)             &
     3454                               ) * flag
     3455
     3456                q(k,j,i)  = q(k,j,i)  + ( sed_qc(k+1) - sed_qc(k) ) *          &
     3457                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
     3458                qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) *          &
     3459                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
     3460                pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) *          &
     3461                                        ddzu(k+1) / hyrho(k) * lv_d_cp *       &
     3462                                        d_exner(k) * dt_micro            * flag
     3463
     3464!
     3465!--             Compute the precipitation rate due to cloud (fog) droplets
     3466                IF ( call_microphysics_at_all_substeps )  THEN
     3467                   prr(k,j,i) = prr(k,j,i) +  sed_qc(k) / hyrho(k)             &
     3468                                * weight_substep(intermediate_timestep_count)  &
     3469                                * flag
     3470                ELSE
     3471                   prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag
     3472                ENDIF
     3473
     3474             ENDDO
     3475          ENDDO
     3476       ENDDO
     3477
     3478       CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )
     3479
     3480    END SUBROUTINE sedimentation_cloud
     3481
     3482
     3483!------------------------------------------------------------------------------!
     3484! Description:
     3485! ------------
     3486!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
    38243487!> Call for grid point i,j
    38253488!------------------------------------------------------------------------------!
     
    39023565
    39033566    END SUBROUTINE sedimentation_cloud_ij
     3567
     3568
     3569!------------------------------------------------------------------------------!
     3570! Description:
     3571! ------------
     3572!> Computation of sedimentation flux. Implementation according to Stevens
     3573!> and Seifert (2008). Code is based on UCLA-LES.
     3574!------------------------------------------------------------------------------!
     3575    SUBROUTINE sedimentation_rain
     3576
     3577       IMPLICIT NONE
     3578
     3579       INTEGER(iwp) ::  i             !< running index x direction
     3580       INTEGER(iwp) ::  j             !< running index y direction
     3581       INTEGER(iwp) ::  k             !< running index z direction
     3582       INTEGER(iwp) ::  k_run         !<
     3583       INTEGER(iwp) ::  m             !< running index surface elements
     3584       INTEGER(iwp) ::  surf_e        !< End index of surface elements at (j,i)-gridpoint
     3585       INTEGER(iwp) ::  surf_s        !< Start index of surface elements at (j,i)-gridpoint
     3586
     3587       REAL(wp)     ::  c_run                      !<
     3588       REAL(wp)     ::  d_max                      !<
     3589       REAL(wp)     ::  d_mean                     !<
     3590       REAL(wp)     ::  d_min                      !<
     3591       REAL(wp)     ::  dr                         !<
     3592       REAL(wp)     ::  flux                       !<
     3593       REAL(wp)     ::  flag                       !< flag to mask topography grid points
     3594       REAL(wp)     ::  lambda_r                   !<
     3595       REAL(wp)     ::  mu_r                       !<
     3596       REAL(wp)     ::  z_run                      !<
     3597
     3598       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_nr     !<
     3599       REAL(wp), DIMENSION(nzb:nzt+1) ::  c_qr     !<
     3600       REAL(wp), DIMENSION(nzb:nzt+1) ::  nr_slope !<
     3601       REAL(wp), DIMENSION(nzb:nzt+1) ::  qr_slope !<
     3602       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nr   !<
     3603       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qr   !<
     3604       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_nr     !<
     3605       REAL(wp), DIMENSION(nzb:nzt+1) ::  w_qr     !<
     3606
     3607       CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )
     3608
     3609!
     3610!--    Compute velocities
     3611       DO  i = nxlg, nxrg
     3612          DO  j = nysg, nyng
     3613             DO  k = nzb+1, nzt
     3614!
     3615!--             Predetermine flag to mask topography
     3616                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     3617
     3618                IF ( qr(k,j,i) > eps_sb )  THEN
     3619!
     3620!--                Weight averaged diameter of rain drops:
     3621                   dr = ( hyrho(k) * qr(k,j,i) /                               &
     3622                          nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )
     3623!
     3624!--                Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
     3625!--                Stevens and Seifert, 2008):
     3626                   mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *                &
     3627                                                     ( dr - 1.4E-3_wp ) ) )
     3628!
     3629!--                Slope parameter of gamma distribution (Seifert, 2008):
     3630                   lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) *        &
     3631                                ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr
     3632
     3633                   w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
     3634                                               a_term - b_term * ( 1.0_wp +    &
     3635                                                  c_term /                     &
     3636                                                  lambda_r )**( -1.0_wp *      &
     3637                                                  ( mu_r + 1.0_wp ) )          &
     3638                                              )                                &
     3639                                ) * flag
     3640
     3641                   w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp,                        &
     3642                                               a_term - b_term * ( 1.0_wp +    &
     3643                                                  c_term /                     &
     3644                                                  lambda_r )**( -1.0_wp *      &
     3645                                                  ( mu_r + 4.0_wp ) )          &
     3646                                             )                                 &
     3647                                ) * flag
     3648                ELSE
     3649                   w_nr(k) = 0.0_wp
     3650                   w_qr(k) = 0.0_wp
     3651                ENDIF
     3652             ENDDO
     3653!
     3654!--          Adjust boundary values using surface data type.
     3655!--          Upward-facing
     3656             surf_s = bc_h(0)%start_index(j,i)
     3657             surf_e = bc_h(0)%end_index(j,i)
     3658             DO  m = surf_s, surf_e
     3659                k         = bc_h(0)%k(m)
     3660                w_nr(k-1) = w_nr(k)
     3661                w_qr(k-1) = w_qr(k)
     3662             ENDDO
     3663!
     3664!--          Downward-facing
     3665             surf_s = bc_h(1)%start_index(j,i)
     3666             surf_e = bc_h(1)%end_index(j,i)
     3667             DO  m = surf_s, surf_e
     3668                k         = bc_h(1)%k(m)
     3669                w_nr(k+1) = w_nr(k)
     3670                w_qr(k+1) = w_qr(k)
     3671             ENDDO
     3672!
     3673!--          Model top boundary value
     3674             w_nr(nzt+1) = 0.0_wp
     3675             w_qr(nzt+1) = 0.0_wp
     3676!
     3677!--          Compute Courant number
     3678             DO  k = nzb+1, nzt
     3679!
     3680!--             Predetermine flag to mask topography
     3681                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     3682
     3683                c_nr(k) = 0.25_wp * ( w_nr(k-1) +                              &
     3684                                      2.0_wp * w_nr(k) + w_nr(k+1) ) *         &
     3685                          dt_micro * ddzu(k) * flag
     3686                c_qr(k) = 0.25_wp * ( w_qr(k-1) +                              &
     3687                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
     3688                          dt_micro * ddzu(k) * flag
     3689             ENDDO
     3690!
     3691!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
     3692             IF ( limiter_sedimentation )  THEN
     3693
     3694                DO k = nzb+1, nzt
     3695!
     3696!--                Predetermine flag to mask topography
     3697                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     3698
     3699                   d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )
     3700                   d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
     3701                   d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
     3702
     3703                   qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
     3704                                                              2.0_wp * d_max,  &
     3705                                                              ABS( d_mean ) )  &
     3706                                                      * flag
     3707
     3708                   d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )
     3709                   d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
     3710                   d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
     3711
     3712                   nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min,  &
     3713                                                              2.0_wp * d_max,  &
     3714                                                              ABS( d_mean ) )
     3715                ENDDO
     3716
     3717             ELSE
     3718
     3719                nr_slope = 0.0_wp
     3720                qr_slope = 0.0_wp
     3721
     3722             ENDIF
     3723
     3724             sed_nr(nzt+1) = 0.0_wp
     3725             sed_qr(nzt+1) = 0.0_wp
     3726!
     3727!--          Compute sedimentation flux
     3728             DO  k = nzt, nzb+1, -1
     3729!
     3730!--             Predetermine flag to mask topography
     3731                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     3732!
     3733!--             Sum up all rain drop number densities which contribute to the flux
     3734!--             through k-1/2
     3735                flux  = 0.0_wp
     3736                z_run = 0.0_wp ! height above z(k)
     3737                k_run = k
     3738                c_run = MIN( 1.0_wp, c_nr(k) )
     3739                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
     3740                   flux  = flux + hyrho(k_run) *                               &
     3741                           ( nr(k_run,j,i) + nr_slope(k_run) *                 &
     3742                           ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run)  &
     3743                                              * flag
     3744                   z_run = z_run + dzu(k_run) * flag
     3745                   k_run = k_run + 1          * flag
     3746                   c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) )    &
     3747                                              * flag
     3748                ENDDO
     3749!
     3750!--             It is not allowed to sediment more rain drop number density than
     3751!--             available
     3752                flux = MIN( flux,                                              &
     3753                            hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) *    &
     3754                            dt_micro                                           &
     3755                          )
     3756
     3757                sed_nr(k) = flux / dt_micro * flag
     3758                nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *          &
     3759                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
     3760!
     3761!--             Sum up all rain water content which contributes to the flux
     3762!--             through k-1/2
     3763                flux  = 0.0_wp
     3764                z_run = 0.0_wp ! height above z(k)
     3765                k_run = k
     3766                c_run = MIN( 1.0_wp, c_qr(k) )
     3767
     3768                DO WHILE ( c_run > 0.0_wp  .AND.  k_run <= nzt )
     3769
     3770                   flux  = flux + hyrho(k_run) * ( qr(k_run,j,i) +             &
     3771                                  qr_slope(k_run) * ( 1.0_wp - c_run ) *       &
     3772                                  0.5_wp ) * c_run * dzu(k_run) * flag
     3773                   z_run = z_run + dzu(k_run)                   * flag
     3774                   k_run = k_run + 1                            * flag
     3775                   c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) )    &
     3776                                                                * flag
     3777
     3778                ENDDO
     3779!
     3780!--             It is not allowed to sediment more rain water content than
     3781!--             available
     3782                flux = MIN( flux,                                              &
     3783                            hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) *      &
     3784                            dt_micro                                           &
     3785                          )
     3786
     3787                sed_qr(k) = flux / dt_micro * flag
     3788
     3789                qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *          &
     3790                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
     3791                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
     3792                                        ddzu(k+1) / hyrho(k) * dt_micro * flag
     3793                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
     3794                                        ddzu(k+1) / hyrho(k) * lv_d_cp *       &
     3795                                        d_exner(k) * dt_micro            * flag
     3796!
     3797!--             Compute the rain rate
     3798                IF ( call_microphysics_at_all_substeps )  THEN
     3799                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k)             &
     3800                                * weight_substep(intermediate_timestep_count) &
     3801                                * flag
     3802                ELSE
     3803                   prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag
     3804                ENDIF
     3805
     3806             ENDDO
     3807          ENDDO
     3808       ENDDO
     3809
     3810       CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )
     3811
     3812    END SUBROUTINE sedimentation_rain
    39043813
    39053814
     
    41234032
    41244033    END SUBROUTINE sedimentation_rain_ij
     4034
     4035
     4036!------------------------------------------------------------------------------!
     4037! Description:
     4038! ------------
     4039!> Computation of the precipitation amount due to gravitational settling of
     4040!> rain and cloud (fog) droplets
     4041!------------------------------------------------------------------------------!
     4042    SUBROUTINE calc_precipitation_amount
     4043
     4044       IMPLICIT NONE
     4045
     4046       INTEGER(iwp) ::  i             !< running index x direction
     4047       INTEGER(iwp) ::  j             !< running index y direction
     4048       INTEGER(iwp) ::  k             !< running index y direction
     4049       INTEGER(iwp) ::  m             !< running index surface elements
     4050
     4051       IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&
     4052            ( .NOT. call_microphysics_at_all_substeps .OR.                     &
     4053            intermediate_timestep_count == intermediate_timestep_count_max ) ) &
     4054       THEN
     4055!
     4056!--       Run over all upward-facing surface elements, i.e. non-natural,
     4057!--       natural and urban
     4058          DO  m = 1, bc_h(0)%ns
     4059             i = bc_h(0)%i(m)
     4060             j = bc_h(0)%j(m)
     4061             k = bc_h(0)%k(m)
     4062             precipitation_amount(j,i) = precipitation_amount(j,i) +           &
     4063                                               prr(k,j,i) * hyrho(k) * dt_3d
     4064          ENDDO
     4065
     4066       ENDIF
     4067
     4068    END SUBROUTINE calc_precipitation_amount
    41254069
    41264070
Note: See TracChangeset for help on using the changeset viewer.