Changeset 2292 for palm/trunk/SOURCE


Ignore:
Timestamp:
Jun 20, 2017 9:51:42 AM (7 years ago)
Author:
schwenkel
Message:

implementation of new bulk microphysics scheme

Location:
palm/trunk/SOURCE
Files:
26 edited

Legend:

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

    r2101 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2101 2017-01-05 16:42:31Z suehring
    2732!
    2833! 2000 2016-08-20 18:09:15Z knoop
     
    10121017          ENDIF
    10131018
    1014        ELSEIF ( sk_char == 'qr' )  THEN
    1015 
    1016 !
    1017 !--       Rain water content boundary condition at the bottom boundary:
     1019       ELSEIF ( sk_char == 'qc' )  THEN
     1020
     1021!
     1022!--       Cloud water content boundary condition at the bottom boundary:
    10181023!--       Dirichlet (fixed surface rain water content).
    10191024          DO  i = nxl, nxr
     
    10261031
    10271032!
     1033!--       Cloud water content boundary condition at the top boundary: Dirichlet
     1034          DO  i = nxl, nxr
     1035             DO  j = nys, nyn
     1036                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     1037                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
     1038             ENDDO
     1039          ENDDO
     1040
     1041       ELSEIF ( sk_char == 'qr' )  THEN
     1042
     1043!
     1044!--       Rain water content boundary condition at the bottom boundary:
     1045!--       Dirichlet (fixed surface rain water content).
     1046          DO  i = nxl, nxr
     1047             DO  j = nys, nyn
     1048                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
     1049                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     1050                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     1051             ENDDO
     1052          ENDDO
     1053
     1054!
    10281055!--       Rain water content boundary condition at the top boundary: Dirichlet
     1056          DO  i = nxl, nxr
     1057             DO  j = nys, nyn
     1058                sk_p(nzt+2,j,i)   = sk_p(nzt+1,j,i)
     1059                sk_p(nzt+3,j,i)   = sk_p(nzt+1,j,i)
     1060             ENDDO
     1061          ENDDO
     1062
     1063       ELSEIF ( sk_char == 'nc' )  THEN
     1064
     1065!
     1066!--       Cloud drop concentration boundary condition at the bottom boundary:
     1067!--       Dirichlet (fixed surface cloud drop concentration).
     1068          DO  i = nxl, nxr
     1069             DO  j = nys, nyn
     1070                sk_p(nzb,j,i)   = sk_p(nzb+1,j,i)
     1071                sk_p(nzb-1,j,i) = sk_p(nzb,j,i)
     1072                sk_p(nzb-2,j,i) = sk_p(nzb,j,i)
     1073             ENDDO
     1074          ENDDO
     1075
     1076!
     1077!--       Cloud drop concentration boundary condition at the top boundary: Dirichlet
    10291078          DO  i = nxl, nxr
    10301079             DO  j = nys, nyn
  • palm/trunk/SOURCE/advec_ws.f90

    r2233 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2233 2017-05-30 18:08:54Z suehring
    2732!
    2833! 2232 2017-05-30 17:47:52Z suehring
     
    239244
    240245       USE arrays_3d,                                                          &
    241            ONLY:  diss_l_e, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr,         &
    242                   diss_l_s, diss_l_sa, diss_l_u, diss_l_v, diss_l_w,  flux_l_e,&
    243                   flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_s,         &
    244                   flux_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_s_e, diss_s_nr,&
    245                   diss_s_pt, diss_s_q, diss_s_qr, diss_s_s, diss_s_sa,         &
    246                   diss_s_u, diss_s_v, diss_s_w, flux_s_e, flux_s_nr, flux_s_pt,&
    247                   flux_s_q, flux_s_qr, flux_s_s, flux_s_sa, flux_s_u, flux_s_v,&
    248                   flux_s_w
     246           ONLY:  diss_l_e, diss_l_nc, diss_l_nr, diss_l_pt, diss_l_q,         &
     247                  diss_l_qc, diss_l_qr, diss_l_s, diss_l_sa, diss_l_u,         &
     248                  diss_l_v, diss_l_w, flux_l_e, flux_l_nc, flux_l_nr,          &
     249                  flux_l_pt, flux_l_q, flux_l_qc, flux_l_qr, flux_l_s,         &
     250                  flux_l_sa, flux_l_u, flux_l_v, flux_l_w, diss_s_e,           &
     251                  diss_s_nc,  diss_s_nr, diss_s_pt, diss_s_q, diss_s_qc,       &
     252                  diss_s_qr, diss_s_s, diss_s_sa, diss_s_u, diss_s_v,          &
     253                  diss_s_w, flux_s_e, flux_s_nc, flux_s_nr, flux_s_pt,         &
     254                  flux_s_q, flux_s_qc, flux_s_qr, flux_s_s, flux_s_sa,         &
     255                  flux_s_u, flux_s_v, flux_s_w
    249256
    250257       USE constants,                                                          &
     
    254261       USE control_parameters,                                                 &
    255262           ONLY:  cloud_physics, humidity, loop_optimization,                  &
    256                   passive_scalar, microphysics_seifert,                        &
     263                  passive_scalar, microphysics_morrison, microphysics_seifert, &
    257264                  ocean, ws_scheme_mom, ws_scheme_sca
    258265
     
    265272
    266273       USE statistics,                                                         &
    267            ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,&
    268                   sums_wspts_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l,            &
    269                   sums_wsss_ws_l, sums_wssas_ws_l,  sums_wsss_ws_l,            &
    270                   sums_wsus_ws_l, sums_wsvs_ws_l 
     274           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsncs_ws_l,&
     275                  sums_wsnrs_ws_l,sums_wspts_ws_l, sums_wsqcs_ws_l,            &
     276                  sums_wsqrs_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l,             &
     277                  sums_wssas_ws_l,  sums_wsss_ws_l, sums_wsus_ws_l,            &
     278                  sums_wsvs_ws_l
     279 
    271280
    272281!
     
    309318             ALLOCATE( sums_wsss_ws_l(nzb:nzt+1,0:threads_per_task-1) )
    310319             sums_wsss_ws_l = 0.0_wp
     320          ENDIF
     321
     322          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     323             ALLOCATE( sums_wsqcs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     324             ALLOCATE( sums_wsncs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
     325             sums_wsqcs_ws_l = 0.0_wp
     326             sums_wsncs_ws_l = 0.0_wp
    311327          ENDIF
    312328
     
    372388                ALLOCATE( flux_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1),    &
    373389                          diss_l_s(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
    374              ENDIF             
     390             ENDIF
     391
     392             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     393                ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1),           &
     394                          diss_s_qc(nzb+1:nzt,0:threads_per_task-1),           &
     395                          flux_s_nc(nzb+1:nzt,0:threads_per_task-1),           &
     396                          diss_s_nc(nzb+1:nzt,0:threads_per_task-1) )
     397                ALLOCATE( flux_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
     398                          diss_l_qc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
     399                          flux_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1),   &
     400                          diss_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
     401             ENDIF                 
    375402
    376403             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    10131040       USE control_parameters,                                                 &
    10141041           ONLY:  cloud_physics, humidity, passive_scalar, ocean,              &
    1015                   microphysics_seifert, ws_scheme_mom, ws_scheme_sca
     1042                  microphysics_morrison, microphysics_seifert, ws_scheme_mom,  &
     1043                  ws_scheme_sca
    10161044
    10171045       USE kinds
    10181046
    10191047       USE statistics,                                                         &
    1020            ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsnrs_ws_l,&
    1021                   sums_wspts_ws_l, sums_wsqrs_ws_l, sums_wsqs_ws_l,            &
    1022                   sums_wsss_ws_l, sums_wssas_ws_l,  sums_wsus_ws_l,            &
    1023                   sums_wsvs_ws_l 
     1048           ONLY:  sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wsncs_ws_l,&
     1049                  sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqcs_ws_l,           &
     1050                  sums_wsqrs_ws_l, sums_wsqs_ws_l, sums_wsss_ws_l,             &
     1051                  sums_wssas_ws_l, sums_wsus_ws_l, sums_wsvs_ws_l     
     1052                   
    10241053
    10251054       IMPLICIT NONE
     
    10401069          IF ( humidity       )  sums_wsqs_ws_l = 0.0_wp
    10411070          IF ( passive_scalar )  sums_wsss_ws_l = 0.0_wp
     1071          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     1072             sums_wsqcs_ws_l = 0.0_wp
     1073             sums_wsncs_ws_l = 0.0_wp
     1074          ENDIF
    10421075          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    10431076             sums_wsqrs_ws_l = 0.0_wp
     
    10811114
    10821115       USE statistics,                                                         &
    1083            ONLY:  hom, sums_wsnrs_ws_l, sums_wspts_ws_l, sums_wsqrs_ws_l,      &
    1084                   sums_wsqs_ws_l, sums_wssas_ws_l, sums_wsss_ws_l,             &
    1085                   weight_substep
     1116           ONLY:  hom, sums_wsncs_ws_l, sums_wsnrs_ws_l, sums_wspts_ws_l,      &
     1117                  sums_wsqcs_ws_l,  sums_wsqrs_ws_l, sums_wsqs_ws_l,           &
     1118                  sums_wssas_ws_l, sums_wsss_ws_l, weight_substep   
     1119                 
    10861120
    10871121       IMPLICIT NONE
     
    15531587             ENDDO
    15541588
     1589          CASE ( 'qc' )
     1590
     1591             DO  k = nzb, nzt
     1592                sums_wsqcs_ws_l(k,tn)  = sums_wsqcs_ws_l(k,tn) +               &
     1593                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
     1594                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
     1595                    + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
     1596                                *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
     1597                    ) * weight_substep(intermediate_timestep_count)
     1598             ENDDO
     1599
     1600
    15551601          CASE ( 'qr' )
    15561602
    15571603             DO  k = nzb, nzt
    15581604                sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn) +               &
     1605                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
     1606                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
     1607                    + diss_t(k) / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
     1608                                *   ABS( w(k,j,i) - hom(k,1,3,0)            )  &
     1609                    ) * weight_substep(intermediate_timestep_count)
     1610             ENDDO
     1611
     1612          CASE ( 'nc' )
     1613
     1614             DO  k = nzb, nzt
     1615                sums_wsncs_ws_l(k,tn)  = sums_wsncs_ws_l(k,tn) +               &
    15591616                    ( flux_t(k) / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
    15601617                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
     
    30893146       USE statistics,                                                        &
    30903147           ONLY:  hom, sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l,      &
    3091                   sums_wsqrs_ws_l, sums_wsnrs_ws_l, sums_wsss_ws_l,           &
    3092                   weight_substep
     3148                  sums_wsqcs_ws_l, sums_wsqrs_ws_l, sums_wsncs_ws_l,          &
     3149                  sums_wsnrs_ws_l, sums_wsss_ws_l, weight_substep 
     3150                 
     3151
    30933152
    30943153       IMPLICIT NONE
     
    35503609                            ) * weight_substep(intermediate_timestep_count)
    35513610                    ENDDO
     3611                 CASE ( 'qc' )
     3612                    DO  k = nzb, nzt
     3613                       sums_wsqcs_ws_l(k,tn)  = sums_wsqcs_ws_l(k,tn)          &
     3614                          + ( flux_t(k)                                        &
     3615                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
     3616                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
     3617                            + diss_t(k)                                        &
     3618                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
     3619                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
     3620                            ) * weight_substep(intermediate_timestep_count)
     3621                    ENDDO
    35523622                 CASE ( 'qr' )
    35533623                    DO  k = nzb, nzt
    35543624                       sums_wsqrs_ws_l(k,tn)  = sums_wsqrs_ws_l(k,tn)          &
     3625                          + ( flux_t(k)                                        &
     3626                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
     3627                                * ( w(k,j,i) - hom(k,1,3,0)                 )  &
     3628                            + diss_t(k)                                        &
     3629                                / ( ABS(w(k,j,i)) + 1.0E-20_wp              )  &
     3630                                *   ABS(w(k,j,i) - hom(k,1,3,0)             )  &
     3631                            ) * weight_substep(intermediate_timestep_count)
     3632                    ENDDO
     3633                 CASE ( 'nc' )
     3634                    DO  k = nzb, nzt
     3635                       sums_wsncs_ws_l(k,tn)  = sums_wsncs_ws_l(k,tn)          &
    35553636                          + ( flux_t(k)                                        &
    35563637                                / ( w(k,j,i) + SIGN( 1.0E-20_wp, w(k,j,i) ) )  &
  • palm/trunk/SOURCE/average_3d_data.f90

    r2233 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2233 2017-05-30 18:08:54Z suehring
    2732!
    2833! 2232 2017-05-30 17:47:52Z suehring
     
    194199             ENDDO
    195200
     201          CASE ( 'nc' )
     202             DO  i = nxlg, nxrg
     203                DO  j = nysg, nyng
     204                   DO  k = nzb, nzt+1
     205                      nc_av(k,j,i) = nc_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     206                   ENDDO
     207                ENDDO
     208             ENDDO
     209
    196210          CASE ( 'nr' )
    197211             DO  i = nxlg, nxrg
  • palm/trunk/SOURCE/boundary_conds.f90

    r2233 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2233 2017-05-30 18:08:54Z suehring
    2732!
    2833! 2232 2017-05-30 17:47:52Z suehring
     
    151156    USE arrays_3d,                                                             &
    152157        ONLY:  c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l,  &
    153                dzu, e_p, nr_p, pt, pt_p, q, q_p, qr_p, s, s_p, sa, sa_p,       &
    154                u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p,                 &
     158               dzu, e_p, nc_p, nr_p, pt, pt_p, q, q_p, qc_p, qr_p, s, s_p, sa, &
     159               sa_p, u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p,           &
    155160               v, vg, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p,                 &
    156161               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s, pt_init
    157 
     162               
    158163    USE control_parameters,                                                    &
    159164        ONLY:  bc_pt_t_val, bc_q_t_val, bc_s_t_val, constant_diffusion,        &
     
    162167               ibc_sa_t, ibc_uv_b, ibc_uv_t, inflow_l, inflow_n, inflow_r,     &
    163168               inflow_s, intermediate_timestep_count, large_scale_forcing,     &
    164                microphysics_seifert, nest_domain, nest_bound_l, nest_bound_s,  &
    165                nudging, ocean, outflow_l, outflow_n, outflow_r, outflow_s,     &
    166                passive_scalar, tsc, use_cmax
     169               microphysics_morrison, microphysics_seifert, nest_domain,       &
     170               nest_bound_l, nest_bound_s, nudging, ocean, outflow_l,          &
     171               outflow_n, outflow_r, outflow_s, passive_scalar, tsc, use_cmax
    167172
    168173    USE grid_variables,                                                        &
     
    393398       ELSEIF ( ibc_q_t == 1 ) THEN
    394399          q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
     400       ENDIF
     401
     402       IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     403!             
     404!--       Surface conditions cloud water (Dirichlet)
     405!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
     406!--       the k coordinate belongs to the atmospheric grid point, therefore, set
     407!--       qr_p and nr_p at k-1
     408          !$OMP PARALLEL DO PRIVATE( i, j, k )
     409          DO  m = 1, bc_h(0)%ns
     410             i = bc_h(0)%i(m)           
     411             j = bc_h(0)%j(m)
     412             k = bc_h(0)%k(m)
     413             qc_p(k-1,j,i) = 0.0_wp
     414             nc_p(k-1,j,i) = 0.0_wp
     415          ENDDO
     416!
     417!--       Top boundary condition for cloud water (Dirichlet)
     418          qc_p(nzt+1,:,:) = 0.0_wp
     419          nc_p(nzt+1,:,:) = 0.0_wp
     420           
    395421       ENDIF
    396422
     
    514540       IF ( humidity )  THEN
    515541          q_p(:,nys-1,:) = q_p(:,nys,:)
     542          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     543             qc_p(:,nys-1,:) = qc_p(:,nys,:)
     544             nc_p(:,nys-1,:) = nc_p(:,nys,:)
     545          ENDIF
    516546          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    517547             qr_p(:,nys-1,:) = qr_p(:,nys,:)
     
    525555       IF ( humidity )  THEN
    526556          q_p(:,nyn+1,:) = q_p(:,nyn,:)
     557          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     558             qc_p(:,nyn+1,:) = qc_p(:,nyn,:)
     559             nc_p(:,nyn+1,:) = nc_p(:,nyn,:)
     560          ENDIF
    527561          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    528562             qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
     
    536570       IF ( humidity )  THEN
    537571          q_p(:,:,nxl-1) = q_p(:,:,nxl)
     572          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     573             qc_p(:,:,nxl-1) = qc_p(:,:,nxl)
     574             nc_p(:,:,nxl-1) = nc_p(:,:,nxl)
     575          ENDIF
    538576          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    539577             qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
     
    547585       IF ( humidity )  THEN
    548586          q_p(:,:,nxr+1) = q_p(:,:,nxr)
     587          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     588             qc_p(:,:,nxr+1) = qc_p(:,:,nxr)
     589             nc_p(:,:,nxr+1) = nc_p(:,:,nxr)
     590          ENDIF
    549591          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    550592             qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
  • palm/trunk/SOURCE/calc_liquid_water_content.f90

    r2233 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2233 2017-05-30 18:08:54Z suehring
    2732!
    2833! 2232 2017-05-30 17:47:52Z suehring
     
    8388
    8489    USE control_parameters,                                                    &
    85         ONLY:  microphysics_seifert
     90        ONLY:  microphysics_morrison, microphysics_seifert
    8691
    8792    USE indices,                                                               &
     
    132137!
    133138!--          Compute the liquid water content
    134              IF ( microphysics_seifert )  THEN
    135                 IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp ) THEN
     139             IF ( microphysics_seifert  .AND.  .NOT. microphysics_morrison )   &
     140             THEN
     141                IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp )  THEN
    136142                   qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) )                  &
    137143                                      * MERGE( 1.0_wp, 0.0_wp,                 &
     
    147153                                               BTEST( wall_flags_0(k,j,i), 0 ) )
    148154                ENDIF
     155             ELSEIF ( microphysics_morrison )  THEN
     156                ql(k,j,i) = qc(k,j,i) + qr(k,j,i)                              &
     157                                      * MERGE( 1.0_wp, 0.0_wp,                 &
     158                                               BTEST( wall_flags_0(k,j,i), 0 ) )
    149159             ELSE
    150                 IF ( ( q(k,j,i) - q_s ) > 0.0_wp ) THEN
     160                IF ( ( q(k,j,i) - q_s ) > 0.0_wp )  THEN
    151161                   qc(k,j,i) = ( q(k,j,i) - q_s )                              &
    152162                                      * MERGE( 1.0_wp, 0.0_wp,                 &
  • palm/trunk/SOURCE/check_parameters.f90

    r2274 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2274 2017-06-09 13:27:48Z Giersch
    2732! Changed error messages
    2833!
     
    983988       microphysics_seifert    = .TRUE.
    984989       microphysics_kessler    = .FALSE.
     990       microphysics_morrison  = .FALSE.
    985991       precipitation           = .TRUE.
    986992    ELSEIF ( cloud_scheme == 'kessler' )  THEN
     
    988994       microphysics_seifert    = .FALSE.
    989995       microphysics_kessler    = .TRUE.
     996       microphysics_morrison   = .FALSE.
     997       precipitation           = .TRUE.
     998    ELSEIF ( cloud_scheme == 'morrison' )  THEN
     999       microphysics_sat_adjust = .FALSE.
     1000       microphysics_seifert    = .TRUE.
     1001       microphysics_kessler    = .FALSE.
     1002       microphysics_morrison   = .TRUE.
    9901003       precipitation           = .TRUE.
    9911004    ELSE
     
    27662779             hom(:,2,122,:) = SPREAD( zw, 2, statistic_regions+1 )
    27672780
     2781          CASE ( 'nc' )
     2782             IF (  .NOT.  cloud_physics )  THEN
     2783                message_string = 'data_output_pr = ' //                        &
     2784                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     2785                                 'lemented for cloud_physics = .FALSE.'
     2786                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
     2787             ELSEIF ( .NOT.  microphysics_morrison )  THEN
     2788                message_string = 'data_output_pr = ' //                        &
     2789                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
     2790                                 'lemented for cloud_scheme /= morrison'
     2791                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
     2792             ELSE
     2793                dopr_index(i) = 89
     2794                dopr_unit(i)  = '1/m3'
     2795                hom(:,2,89,:)  = SPREAD( zu, 2, statistic_regions+1 )
     2796             ENDIF
     2797
    27682798          CASE ( 'nr' )
    27692799             IF (  .NOT.  cloud_physics )  THEN
     
    30613091             ENDIF
    30623092             unit = 'K'
     3093
     3094          CASE ( 'nc' )
     3095             IF (  .NOT.  cloud_physics )  THEN
     3096                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3097                         'res cloud_physics = .TRUE.'
     3098                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
     3099             ELSEIF ( .NOT.  microphysics_morrison )  THEN
     3100                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     3101                         'res = microphysics morrison '
     3102                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
     3103             ENDIF
     3104             unit = '1/m3'
    30633105
    30643106          CASE ( 'nr' )
  • palm/trunk/SOURCE/data_output_2d.f90

    r2277 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2277 2017-06-12 10:47:51Z kanani
    2732! Removed unused variables do2d_xy_n, do2d_xz_n, do2d_yz_n
    2833!
     
    169174
    170175    USE arrays_3d,                                                             &
    171         ONLY:  dzw, e, nr, p, pt, precipitation_amount, precipitation_rate,    &
     176        ONLY:  dzw, e, nc, nr, p, pt, precipitation_amount, precipitation_rate,&
    172177               prr, q, qc, ql, ql_c, ql_v, ql_vp, qr, rho_ocean, s, sa,        &
    173178               tend, u, v, vpt, w, zu, zw
     
    492497                two_d = .TRUE.
    493498                level_z(nzb+1) = zu(nzb+1)
     499
     500             CASE ( 'nc_xy', 'nc_xz', 'nc_yz' )
     501                IF ( av == 0 )  THEN
     502                   to_be_resorted => nc
     503                ELSE
     504                   to_be_resorted => nc_av
     505                ENDIF
     506                IF ( mode == 'xy' )  level_z = zu
    494507
    495508             CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
  • palm/trunk/SOURCE/data_output_3d.f90

    r2233 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2233 2017-05-30 18:08:54Z suehring
    2732!
    2833! 2232 2017-05-30 17:47:52Z suehring
     
    155160
    156161    USE arrays_3d,                                                             &
    157         ONLY:  e, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa, &
    158                tend, u, v, vpt, w
     162        ONLY:  e, nc, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, &
     163               sa, tend, u, v, vpt, w
    159164       
    160165    USE averaging
     
    325330             ELSE
    326331                to_be_resorted => lpt_av
     332             ENDIF
     333
     334          CASE ( 'nc' )
     335             IF ( av == 0 )  THEN
     336                to_be_resorted => nc
     337             ELSE
     338                to_be_resorted => nc_av
    327339             ENDIF
    328340
  • palm/trunk/SOURCE/data_output_mask.f90

    r2101 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2101 2017-01-05 16:42:31Z suehring
    2732!
    2833! 2031 2016-10-21 15:11:58Z knoop
     
    105110#if defined( __netcdf )
    106111    USE arrays_3d,                                                             &
    107         ONLY:  e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa, tend, u,   &
    108                v, vpt, w
     112        ONLY:  e, nc, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa,  &
     113               tend, u, v, vpt, w
    109114   
    110115    USE averaging,                                                             &
    111         ONLY:  e_av, lpt_av, nr_av, p_av, pc_av, pr_av, pt_av, q_av, qc_av,    &
    112                ql_av, ql_c_av, ql_v_av, ql_vp_av, qv_av, qr_av, rho_ocean_av, s_av,  &
    113                sa_av, u_av, v_av, vpt_av, w_av
     116        ONLY:  e_av, lpt_av, nc_av, nr_av, p_av, pc_av, pr_av, pt_av, q_av,    &
     117               qc_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qv_av, qr_av,         &
     118               rho_ocean_av, s_av, sa_av, u_av, v_av, vpt_av, w_av
    114119   
    115120    USE cloud_parameters,                                                      &
     
    243248             ELSE
    244249                to_be_resorted => lpt_av
     250             ENDIF
     251
     252          CASE ( 'nc' )
     253             IF ( av == 0 )  THEN
     254                to_be_resorted => nc
     255             ELSE
     256                to_be_resorted => nc_av
    245257             ENDIF
    246258
  • palm/trunk/SOURCE/flow_statistics.f90

    r2270 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2270 2017-06-09 12:18:47Z maronga
    2732! Revised numbering (removed 2 timeseries)
    2833!
     
    233238    USE arrays_3d,                                                             &
    234239        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
    235                momentumflux_output_conversion, nr, p, prho, prr, pt, q,        &
     240               momentumflux_output_conversion, nc, nr, p, prho, prr, pt, q,    &
    236241               qc, ql, qr, rho_air, rho_air_zw, rho_ocean, s,                  &
    237242               sa, td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert, u,   &
     
    245250                dt_3d, g, humidity, kappa, land_surface, large_scale_forcing,  &
    246251                large_scale_subsidence, max_pr_user, message_string, neutral,  &
    247                 microphysics_seifert, ocean, passive_scalar, simulated_time,   &
    248                 use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
    249                 ws_scheme_mom, ws_scheme_sca
     252                microphysics_morrison, microphysics_seifert, ocean,            &
     253                passive_scalar, simulated_time, use_subsidence_tendencies,     &
     254                use_surface_fluxes, use_top_fluxes, ws_scheme_mom,             &
     255                ws_scheme_sca
    250256       
    251257    USE cpulog,                                                                &
     
    12171223                                                             rmask(j,i,sr) *   &
    12181224                                                             flag
     1225                         IF ( microphysics_morrison )  THEN
     1226                            sums_l(k,123,tn) = sums_l(k,123,tn) + nc(k,j,i) *  &
     1227                                                                rmask(j,i,sr) *&
     1228                                                                flag
     1229                         ENDIF
    12191230                         IF ( microphysics_seifert )  THEN
    12201231                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
     
    17591770       hom(:,1,71,sr) = sums(:,71)     ! prho
    17601771       hom(:,1,72,sr) = hyp * 1E-2_wp  ! hyp in hPa
     1772       hom(:,1,123,sr) = sums(:,123)   ! nc
    17611773       hom(:,1,73,sr) = sums(:,73)     ! nr
    17621774       hom(:,1,74,sr) = sums(:,74)     ! qr
  • palm/trunk/SOURCE/init_3d_model.f90

    r2277 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2277 2017-06-12 10:47:51Z kanani
    2732! Removed unused variable sums_up_fraction_l
    2833!
     
    560565
    561566          IF ( cloud_physics )  THEN
    562 
    563567!
    564568!--          Liquid water content
     
    568572             ALLOCATE ( ql_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    569573#endif
     574
     575!
     576!--          3D-cloud water content
     577             IF ( .NOT. microphysics_morrison )  THEN
     578#if defined( __nopointer )
     579                ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     580#else
     581                ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     582#endif
     583             ENDIF
    570584!
    571585!--          Precipitation amount and rate (only needed if output is switched)
     
    574588
    575589!
    576 !--          3D-cloud water content
    577 #if defined( __nopointer )
    578              ALLOCATE( qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    579 #else
    580              ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    581 #endif
    582 !
    583590!--          3d-precipitation rate
    584591             ALLOCATE( prr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     592
     593             IF ( microphysics_morrison )  THEN
     594!
     595!--             3D-cloud drop water content, cloud drop concentration arrays
     596#if defined( __nopointer )
     597                ALLOCATE( nc(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
     598                          nc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     599                          qc(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                &
     600                          qc_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     601                          tnc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)
     602                          tqc_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     603#else
     604                ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     605                          nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     606                          nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     607                          qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     608                          qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg),              &
     609                          qc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     610#endif
     611             ENDIF
    585612
    586613             IF ( microphysics_seifert )  THEN
     
    880907          IF ( cloud_physics )  THEN
    881908             ql => ql_1
    882              qc => qc_1
     909             IF ( .NOT. microphysics_morrison )  THEN
     910                qc => qc_1
     911             ENDIF
     912             IF ( microphysics_morrison )  THEN
     913                qc => qc_1;  qc_p  => qc_2;  tqc_m  => qc_3
     914                nc => nc_1;  nc_p  => nc_2;  tnc_m  => nc_3
     915             ENDIF
    883916             IF ( microphysics_seifert )  THEN
    884917                qr => qr_1;  qr_p  => qr_2;  tqr_m  => qr_3
     
    9791012                ENDDO
    9801013             ENDDO
     1014             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     1015                DO  i = nxlg, nxrg
     1016                   DO  j = nysg, nyng
     1017                      qc(:,j,i) = 0.0_wp
     1018                      nc(:,j,i) = 0.0_wp
     1019                   ENDDO
     1020                ENDDO
     1021             ENDIF
    9811022             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    9821023                DO  i = nxlg, nxrg
     
    9861027                   ENDDO
    9871028                ENDDO
    988 
    9891029             ENDIF
    9901030          ENDIF
     1031
    9911032          IF ( passive_scalar )  THEN
    9921033             DO  i = nxlg, nxrg
     
    11071148                ENDDO
    11081149             ENDDO
     1150             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     1151                DO  i = nxlg, nxrg
     1152                   DO  j = nysg, nyng
     1153                      qc(:,j,i) = 0.0_wp
     1154                      nc(:,j,i) = 0.0_wp
     1155                   ENDDO
     1156                ENDDO
     1157             ENDIF
     1158
    11091159             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    1110 
    11111160                DO  i = nxlg, nxrg
    11121161                   DO  j = nysg, nyng
     
    11151164                   ENDDO
    11161165                ENDDO
    1117 
    11181166             ENDIF
     1167
    11191168          ENDIF
    11201169         
     
    13191368          tq_m = 0.0_wp
    13201369          q_p = q
     1370          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     1371             tqc_m = 0.0_wp
     1372             qc_p  = qc
     1373             tnc_m = 0.0_wp
     1374             nc_p  = nc
     1375          ENDIF
    13211376          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    13221377             tqr_m = 0.0_wp
     
    15321587       IF ( humidity )  THEN
    15331588          q_p = q
     1589          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     1590             qc_p = qc
     1591             nc_p = nc
     1592          ENDIF
    15341593          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    15351594             qr_p = qr
     
    15471606       IF ( humidity )  THEN
    15481607          tq_m = 0.0_wp
     1608          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     1609             tqc_m = 0.0_wp
     1610             tnc_m = 0.0_wp
     1611          ENDIF
    15491612          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    15501613             tqr_m = 0.0_wp
  • palm/trunk/SOURCE/init_masks.f90

    r2271 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2271 2017-06-09 12:34:55Z sward
    2732! Changed error message
    2833!
     
    115120               mask_x_loop, mask_xyz_dimension, mask_y, mask_y_loop, mask_z,   &
    116121               mask_z_loop, max_masks,  message_string, mid,                   &
    117                microphysics_seifert, passive_scalar, ocean
     122               microphysics_morrison, microphysics_seifert, passive_scalar,    &
     123               ocean
    118124
    119125    USE grid_variables,                                                        &
     
    278284                unit = 'K'
    279285
     286             CASE ( 'nc' )
     287                IF ( .NOT. cloud_physics )  THEN
     288                   WRITE ( message_string, * ) 'output of "', TRIM( var ), &
     289                        '" requires cloud_physics = .TRUE.'
     290                   CALL message( 'init_masks', 'PA0108', 1, 2, 0, 6, 0 )
     291                 ELSEIF ( .NOT. microphysics_morrison ) THEN
     292                   message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     293                         'res  = morrison'
     294                   CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
     295                ENDIF
     296                unit = '1/m3'
     297
    280298             CASE ( 'nr' )
    281299                IF ( .NOT. cloud_physics )  THEN
  • palm/trunk/SOURCE/microphysics_mod.f90

    r2233 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30! - The process of activation is parameterized with a simple Twomey
     31!   activion scheme or with considering solution and curvature
     32!   effects (Khvorostyanov and Curry ,2006).
     33! - The saturation adjustment scheme is replaced by the parameterization
     34!   of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128).
     35! - All other microphysical processes of Seifert and Beheng are used.
     36!   Additionally, in those processes the reduction of cloud number concentration
     37!   is considered.
     38!
     39! 2233 2017-05-30 18:08:54Z suehring
    2740!
    2841! 2232 2017-05-30 17:47:52Z suehring
     
    123136    IMPLICIT NONE
    124137
    125     LOGICAL ::  cloud_water_sedimentation = .FALSE.  !< cloud water sedimentation
    126     LOGICAL ::  limiter_sedimentation = .TRUE.       !< sedimentation limiter
    127     LOGICAL ::  collision_turbulence = .FALSE.       !< turbulence effects
    128     LOGICAL ::  ventilation_effect = .TRUE.          !< ventilation effect
     138    LOGICAL ::  cloud_water_sedimentation = .FALSE.       !< cloud water sedimentation
     139    LOGICAL ::  curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory
     140    LOGICAL ::  limiter_sedimentation = .TRUE.            !< sedimentation limiter
     141    LOGICAL ::  collision_turbulence = .FALSE.            !< turbulence effects
     142    LOGICAL ::  ventilation_effect = .TRUE.               !< ventilation effect
    129143
    130144    REAL(wp) ::  a_1 = 8.69E-4_wp          !< coef. in turb. parametrization (cm-2 s3)
     
    147161    REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
    148162    REAL(wp) ::  eps_sb = 1.0E-10_wp       !< threshold in two-moments scheme
     163    REAL(wp) ::  eps_mr = 0.0_wp           !< threshold for morrison scheme
    149164    REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
    150165    REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
     
    161176    REAL(wp) ::  w_precipitation = 9.65_wp  !< maximum terminal velocity (m s-1)
    162177    REAL(wp) ::  x0 = 2.6E-10_wp           !< separating drop mass (kg)
     178    REAL(wp) ::  xamin = 5.24E-19_wp       !< average aerosol mass (kg) (~ 0.05µm)
     179    REAL(wp) ::  xcmin = 4.18E-15_wp       !< minimum cloud drop size (kg) (~ 1µm)
     180    REAL(wp) ::  xcmax = 2.6E-10_wp        !< maximum cloud drop size (kg) (~ 40µm)
    163181    REAL(wp) ::  xrmin = 2.6E-10_wp        !< minimum rain drop size (kg)
    164182    REAL(wp) ::  xrmax = 5.0E-6_wp         !< maximum rain drop site (kg)
     
    167185    REAL(wp) ::  dpirho_l                  !< 6.0 / ( pi * rho_l )
    168186    REAL(wp) ::  dt_micro                  !< microphysics time step
     187    REAL(wp) ::  na_init = 100.0E6_wp      !< Total particle/aerosol concentration (cm-3)
    169188    REAL(wp) ::  nc_const = 70.0E6_wp      !< cloud droplet concentration
    170189    REAL(wp) ::  dt_precipitation = 100.0_wp !< timestep precipitation (s)
     
    184203    PUBLIC microphysics_control, microphysics_init
    185204
    186     PUBLIC cloud_water_sedimentation, collision_turbulence, c_sedimentation,   &
    187            dt_precipitation, limiter_sedimentation, nc_const, sigma_gc,        &
     205    PUBLIC cloud_water_sedimentation, collision_turbulence,                    &
     206           curvature_solution_effects_bulk, c_sedimentation, dt_precipitation, &
     207           limiter_sedimentation, na_init, nc_const, sigma_gc,                 &
    188208           ventilation_effect
     209           
    189210
    190211    INTERFACE microphysics_control
     
    198219    END INTERFACE adjust_cloud
    199220
     221    INTERFACE activation
     222       MODULE PROCEDURE activation
     223       MODULE PROCEDURE activation_ij
     224    END INTERFACE activation
     225
     226    INTERFACE condensation
     227       MODULE PROCEDURE condensation
     228       MODULE PROCEDURE condensation_ij
     229    END INTERFACE condensation
     230
    200231    INTERFACE autoconversion
    201232       MODULE PROCEDURE autoconversion
     
    256287
    257288       USE control_parameters,                                                 &
    258            ONLY:  microphysics_seifert
     289           ONLY:  microphysics_morrison, microphysics_seifert
    259290
    260291       USE indices,                                                            &
     
    283314!
    284315!--    Allocate 1D microphysics arrays
    285        ALLOCATE ( nc_1d(nzb:nzt+1), pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1),         &
     316       ALLOCATE ( pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1),         &
    286317                  qc_1d(nzb:nzt+1) )
     318
     319       IF ( microphysics_morrison .OR. microphysics_seifert ) THEN
     320          ALLOCATE ( nc_1d(nzb:nzt+1) )
     321       ENDIF
    287322
    288323       IF ( microphysics_seifert )  THEN
     
    290325       ENDIF
    291326
    292 !
    293 !--    Initialze nc_1d with nc_const
    294        nc_1d = nc_const
    295 
    296327    END SUBROUTINE microphysics_init
    297328
     
    313344           ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
    314345                  intermediate_timestep_count, large_scale_forcing,            &
    315                   lsf_surf, microphysics_kessler, microphysics_seifert,        &
    316                   pt_surface, rho_surface,surface_pressure
     346                  lsf_surf, microphysics_kessler, microphysics_morrison,       &
     347                  microphysics_seifert, pt_surface,                            &
     348                  rho_surface,surface_pressure
    317349
    318350       USE indices,                                                            &
     
    373405
    374406          CALL adjust_cloud
     407          IF ( microphysics_morrison )  CALL activation
     408          IF ( microphysics_morrison )  CALL condensation
    375409          CALL autoconversion
    376410          CALL accretion
     
    396430
    397431       USE arrays_3d,                                                          &
    398            ONLY:  qr, nr
     432           ONLY:  qc, nc, qr, nr
    399433
    400434       USE cloud_parameters,                                                   &
    401435           ONLY:  hyrho
     436
     437       USE control_parameters,                                                 &
     438           ONLY:  microphysics_morrison
    402439
    403440       USE cpulog,                                                             &
     
    434471                   ENDIF
    435472                ENDIF
     473                IF ( microphysics_morrison ) THEN
     474                   IF ( qc(k,j,i) <= eps_sb )  THEN
     475                      qc(k,j,i) = 0.0_wp
     476                      nc(k,j,i) = 0.0_wp
     477                   ELSE
     478                      IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
     479                         nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin *            &
     480                                          MERGE( 1.0_wp, 0.0_wp,               &           
     481                                              BTEST( wall_flags_0(k,j,i), 0 ) )
     482                      ENDIF
     483                   ENDIF
     484                ENDIF
    436485             ENDDO
    437486          ENDDO
     
    442491    END SUBROUTINE adjust_cloud
    443492
     493!------------------------------------------------------------------------------!
     494! Description:
     495! ------------
     496!> Calculate number of activated condensation nucleii after simple activation
     497!> scheme of Twomey, 1959.
     498!------------------------------------------------------------------------------!
     499    SUBROUTINE activation
     500
     501       USE arrays_3d,                                                          &
     502           ONLY:  hyp, nc, nr, pt, q,  qc, qr
     503
     504       USE cloud_parameters,                                                   &
     505           ONLY:  hyrho, l_d_cp, l_d_r, l_v, rho_l, r_v, t_d_pt
     506
     507       USE constants,                                                          &
     508           ONLY:  pi
     509
     510       USE cpulog,                                                             &
     511           ONLY:  cpu_log, log_point_s
     512
     513       USE indices,                                                            &
     514           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
     515
     516       USE kinds
     517
     518       USE control_parameters,                                                 &
     519          ONLY: simulated_time
     520
     521       USE particle_attributes,                                                &
     522          ONLY:  molecular_weight_of_solute, molecular_weight_of_water, rho_s, &
     523              s1, s2, s3, vanthoff
     524
     525       IMPLICIT NONE
     526
     527       INTEGER(iwp) ::  i                 !<
     528       INTEGER(iwp) ::  j                 !<
     529       INTEGER(iwp) ::  k                 !<
     530
     531       REAL(wp)     ::  activ             !<
     532       REAL(wp)     ::  afactor           !<
     533       REAL(wp)     ::  alpha             !<
     534       REAL(wp)     ::  beta_act          !<
     535       REAL(wp)     ::  bfactor           !<
     536       REAL(wp)     ::  e_s               !<
     537       REAL(wp)     ::  k_act             !<
     538       REAL(wp)     ::  n_act             !<
     539       REAL(wp)     ::  n_ccn             !<
     540       REAL(wp)     ::  q_s               !<
     541       REAL(wp)     ::  rd0               !<
     542       REAL(wp)     ::  s_0               !<
     543       REAL(wp)     ::  sat               !<
     544       REAL(wp)     ::  sat_max           !<
     545       REAL(wp)     ::  sigma             !<
     546       REAL(wp)     ::  t_int             !<
     547       REAL(wp)     ::  t_l               !<
     548
     549       CALL cpu_log( log_point_s(65), 'activation', 'start' )
     550
     551       DO  i = nxlg, nxrg
     552          DO  j = nysg, nyng
     553             DO  k = nzb+1, nzt
     554
     555!
     556!--             Actual liquid water temperature:
     557                t_l = t_d_pt(k) * pt(k,j,i)
     558
     559!
     560!--             Calculate actual temperature
     561                t_int = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp
     562!
     563!--             Saturation vapor pressure at t_l:
     564                e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
     565                                       ( t_l - 35.86_wp )                   &
     566                                     )
     567!
     568!--             Computation of saturation humidity:
     569                q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
     570                alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
     571                q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /               &
     572                        ( 1.0_wp + alpha * q_s )
     573
     574!--             Supersaturation:
     575                sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
     576
     577!
     578!--             Prescribe parameters for activation
     579!--             (see: Bott + Trautmann, 2002, Atm. Res., 64)
     580                k_act  = 0.7_wp
     581
     582                IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN
     583!
     584!--                Compute the number of activated Aerosols
     585!--                (see: Twomey, 1959, Pure and applied Geophysics, 43)
     586                   n_act     = na_init * sat**k_act
     587!
     588!--                Compute the number of cloud droplets
     589!--                (see: Morrison + Grabowski, 2007, JAS, 64)
     590!                  activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro
     591
     592!
     593!--                Compute activation rate after Khairoutdinov and Kogan
     594!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)   
     595                   sat_max = 1.0_wp / 100.0_wp                 
     596                   activ   = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN       &
     597                      ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) /    &
     598                       dt_micro
     599                ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk ) THEN
     600!
     601!--                Curvature effect (afactor) with surface tension
     602!--                parameterization by Straka (2009)
     603                   sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
     604                   afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
     605!
     606!--                Solute effect (bfactor)
     607                   bfactor = vanthoff * molecular_weight_of_water *            &
     608                       rho_s / ( molecular_weight_of_solute * rho_l )
     609
     610!
     611!--                Prescribe power index that describes the soluble fraction
     612!--                of an aerosol particle (beta) and mean geometric radius of
     613!--                dry aerosol spectrum
     614!--                (see: Morrison + Grabowski, 2007, JAS, 64)
     615                   beta_act = 0.5_wp
     616                   rd0      = 0.05E-6_wp
     617!
     618!--                Calculate mean geometric supersaturation (s_0) with
     619!--                parameterization by Khvorostyanov and Curry (2006)
     620                   s_0   = rd0 **(- ( 1.0_wp + beta_act ) ) *                  &
     621                       ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
     622
     623!
     624!--                Calculate number of activated CCN as a function of
     625!--                supersaturation and taking Koehler theory into account
     626!--                (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)       
     627                   n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(              &
     628                      LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(s1) ) ) )     
     629                   activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp )
     630                ENDIF
     631
     632                nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)
     633
     634             ENDDO
     635          ENDDO
     636       ENDDO
     637
     638       CALL cpu_log( log_point_s(65), 'activation', 'stop' )
     639
     640    END SUBROUTINE activation
     641
     642
     643!------------------------------------------------------------------------------!
     644! Description:
     645! ------------
     646!> Calculate condensation rate for cloud water content (after Khairoutdinov and 
     647!> Kogan, 2000).
     648!------------------------------------------------------------------------------!
     649    SUBROUTINE condensation
     650
     651       USE arrays_3d,                                                          &
     652           ONLY:  hyp, nr, pt, q,  qc, qr, nc
     653
     654       USE cloud_parameters,                                                   &
     655           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
     656
     657       USE constants,                                                          &
     658           ONLY:  pi
     659
     660       USE cpulog,                                                             &
     661           ONLY:  cpu_log, log_point_s
     662
     663       USE indices,                                                            &
     664           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
     665
     666       USE kinds
     667
     668       USE control_parameters,                                                 &
     669          ONLY: simulated_time
     670
     671       IMPLICIT NONE
     672
     673       INTEGER(iwp) ::  i                 !<
     674       INTEGER(iwp) ::  j                 !<
     675       INTEGER(iwp) ::  k                 !<
     676
     677       REAL(wp)     ::  alpha             !<
     678       REAL(wp)     ::  cond              !<
     679       REAL(wp)     ::  cond_max          !<
     680       REAL(wp)     ::  dc                !<
     681       REAL(wp)     ::  e_s               !<
     682       REAL(wp)     ::  evap              !<
     683       REAL(wp)     ::  evap_nc           !<
     684       REAL(wp)     ::  g_fac             !<
     685       REAL(wp)     ::  nc_0              !<
     686       REAL(wp)     ::  q_s               !<
     687       REAL(wp)     ::  sat               !<
     688       REAL(wp)     ::  t_l               !<
     689       REAL(wp)     ::  temp              !<
     690       REAL(wp)     ::  xc                !<
     691
     692       CALL cpu_log( log_point_s(66), 'condensation', 'start' )
     693
     694       DO  i = nxlg, nxrg
     695          DO  j = nysg, nyng
     696             DO  k = nzb+1, nzt
     697!
     698!--             Actual liquid water temperature:
     699                t_l = t_d_pt(k) * pt(k,j,i)
     700!
     701!--             Saturation vapor pressure at t_l:
     702                e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /       &
     703                                       ( t_l - 35.86_wp )                      &
     704                                     )
     705!
     706!--             Computation of saturation humidity:
     707                q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
     708                alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
     709                q_s   = q_s * ( 1.0_wp + alpha * q(k,j,i) ) /                  &
     710                        ( 1.0_wp + alpha * q_s )
     711
     712!--             Supersaturation:
     713                sat   = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp
     714
     715!
     716!--             Actual temperature:
     717                temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
     718
     719                g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
     720                                    l_v / ( thermal_conductivity_l * temp )    &
     721                                    + r_v * temp / ( diff_coeff_l * e_s )      &
     722                                    )
     723!
     724!--             Mean weight of cloud drops
     725                IF ( nc(k,j,i) <= 0.0_wp) CYCLE
     726                xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)               
     727!
     728!--             Weight averaged diameter of cloud drops:
     729                dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
     730!
     731!--             Integral diameter of cloud drops
     732                nc_0 = nc(k,j,i) * dc               
     733!
     734!--             Condensation needs only to be calculated in supersaturated regions
     735                IF ( sat > 0.0_wp )  THEN
     736!
     737!--                Condensation rate of cloud water content
     738!--                after KK scheme.
     739!--                (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
     740                   cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
     741                   cond_max  = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)
     742                   cond      = MIN( cond, cond_max / dt_micro )
     743               
     744                   qc(k,j,i) = qc(k,j,i) + cond * dt_micro
     745                ELSEIF ( sat < 0.0_wp ) THEN
     746                   evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
     747                   evap      = MAX( evap, -qc(k,j,i) / dt_micro )
     748
     749                   qc(k,j,i) = qc(k,j,i) + evap * dt_micro
     750                ENDIF
     751                IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) )  THEN
     752                   nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin
     753                ENDIF
     754             ENDDO
     755          ENDDO
     756       ENDDO
     757
     758       CALL cpu_log( log_point_s(66), 'condensation', 'stop' )
     759
     760    END SUBROUTINE condensation
     761
    444762
    445763!------------------------------------------------------------------------------!
     
    451769
    452770       USE arrays_3d,                                                          &
    453            ONLY:  diss, dzu, nr, qc, qr
     771           ONLY:  diss, dzu, nc, nr, qc, qr
    454772
    455773       USE cloud_parameters,                                                   &
     
    457775
    458776       USE control_parameters,                                                 &
    459            ONLY:  rho_surface
     777           ONLY:  microphysics_morrison, rho_surface
    460778
    461779       USE cpulog,                                                             &
     
    482800       REAL(wp)     ::  k_au              !<
    483801       REAL(wp)     ::  l_mix             !<
     802       REAL(wp)     ::  nc_auto           !<
    484803       REAL(wp)     ::  nu_c              !<
    485804       REAL(wp)     ::  phi_au            !<
     
    500819                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    501820
     821                nc_auto = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
     822
    502823                IF ( qc(k,j,i) > eps_sb )  THEN
    503824
     
    518839!
    519840!--                Mean weight of cloud droplets:
    520                    xc = hyrho(k) * qc(k,j,i) / nc_const
     841                   xc = hyrho(k) * qc(k,j,i) / nc_auto
    521842!
    522843!--                Parameterized turbulence effects on autoconversion (Seifert,
     
    569890                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro  &
    570891                                                              * flag
     892                   IF ( microphysics_morrison ) THEN
     893                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp *         &
     894                                    autocon / x0 * hyrho(k) * dt_micro * flag )
     895                   ENDIF
    571896
    572897                ENDIF
     
    654979
    655980       USE arrays_3d,                                                          &
    656            ONLY:  diss, qc, qr
     981           ONLY:  diss, qc, qr, nc
    657982
    658983       USE cloud_parameters,                                                   &
     
    660985
    661986       USE control_parameters,                                                 &
    662            ONLY:  rho_surface
     987           ONLY:  microphysics_morrison, rho_surface
    663988
    664989       USE cpulog,                                                             &
     
    6791004       REAL(wp)     ::  flag              !< flag to mask topography grid points
    6801005       REAL(wp)     ::  k_cr              !<
     1006       REAL(wp)     ::  nc_accr           !<
    6811007       REAL(wp)     ::  phi_ac            !<
    6821008       REAL(wp)     ::  tau_cloud         !<
     1009       REAL(wp)     ::  xc                !<
     1010
    6831011
    6841012       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
     
    6911019                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    6921020
    693                 IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb ) )  THEN
     1021                nc_accr = MERGE( nc(k,j,i), nc_const, microphysics_morrison )
     1022
     1023                IF ( ( qc(k,j,i) > eps_sb )  .AND.  ( qr(k,j,i) > eps_sb )    &
     1024                                             .AND.  ( nc_accr > eps_mr ) ) THEN
    6941025!
    6951026!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
     
    6991030!--                Beheng, 2001):
    7001031                   phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
     1032
     1033!
     1034!--                Mean weight of cloud drops
     1035                   xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)     
    7011036!
    7021037!--                Parameterized turbulence effects on autoconversion (Seifert,
     
    7191054                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag
    7201055                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag
     1056                   IF ( microphysics_morrison )  THEN
     1057                      nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i),                  &
     1058                                         accr / xc * hyrho(k) * dt_micro * flag)
     1059                   ENDIF
    7211060
    7221061                ENDIF
     
    9761315
    9771316       USE arrays_3d,                                                          &
    978            ONLY:  ddzu, dzu, pt, prr, q, qc
     1317           ONLY:  ddzu, dzu, nc, pt, prr, q, qc
    9791318
    9801319       USE cloud_parameters,                                                   &
     
    9821321
    9831322       USE control_parameters,                                                 &
    984            ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
     1323           ONLY:  call_microphysics_at_all_substeps,                           &
     1324                  intermediate_timestep_count, microphysics_morrison
    9851325
    9861326       USE cpulog,                                                             &
     
    10021342       INTEGER(iwp) ::  k             !<
    10031343
    1004        REAL(wp)                       ::  flag   !< flag to mask topography grid points
     1344       REAL(wp) ::  flag              !< flag to mask topography grid points
     1345       REAL(wp) ::  nc_sedi           !<
     1346
    10051347       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc !<
     1348       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc !<
     1349
    10061350
    10071351       CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )
     
    10151359!--             Predetermine flag to mask topography
    10161360                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     1361                nc_sedi = MERGE ( nc(k,j,i), nc_const, microphysics_morrison )
     1362
     1363!
     1364!--             Sedimentation fluxes for number concentration are only calculated
     1365!--             for cloud_scheme = 'morrison'
     1366                IF ( microphysics_morrison ) THEN
     1367                   IF ( qc(k,j,i) > eps_sb  .AND.  nc(k,j,i) > eps_mr )  THEN
     1368                      sed_nc(k) = sed_qc_const *                               &
     1369                              ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *  &
     1370                              ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )
     1371                   ELSE
     1372                      sed_nc(k) = 0.0_wp
     1373                   ENDIF
     1374
     1375                   sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *           &
     1376                                    nc(k,j,i) / dt_micro + sed_nc(k+1)         &
     1377                                  ) * flag
     1378
     1379                   nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) *       &
     1380                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
     1381                ENDIF
    10171382
    10181383                IF ( qc(k,j,i) > eps_sb )  THEN
    1019                    sed_qc(k) = sed_qc_const * nc_const**( -2.0_wp / 3.0_wp ) * &
     1384                   sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) * &
    10201385                               ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &
    10211386                                                                           flag
     
    13891754
    13901755       USE arrays_3d,                                                          &
    1391            ONLY:  hyp, nr, pt, pt_init, prr, q, qc, qr, zu
     1756           ONLY:  hyp, nc, nr, pt, pt_init, prr, q, qc, qr, zu
    13921757
    13931758       USE cloud_parameters,                                                   &
     
    13971762           ONLY:  call_microphysics_at_all_substeps, dt_3d, g,                 &
    13981763                  intermediate_timestep_count, large_scale_forcing,            &
    1399                   lsf_surf, microphysics_seifert, microphysics_kessler,        &
    1400                   pt_surface, rho_surface, surface_pressure
     1764                  lsf_surf, microphysics_morrison, microphysics_seifert,       &
     1765                  microphysics_kessler, pt_surface, rho_surface,               &
     1766                  surface_pressure
    14011767
    14021768       USE indices,                                                            &
     
    14481814       pt_1d(:) = pt(:,j,i)
    14491815       qc_1d(:) = qc(:,j,i)
    1450        nc_1d(:) = nc_const
    14511816       IF ( microphysics_seifert )  THEN
    14521817          qr_1d(:) = qr(:,j,i)
    14531818          nr_1d(:) = nr(:,j,i)
    14541819       ENDIF
     1820       IF ( microphysics_morrison )  THEN
     1821          nc_1d(:) = nc(:,j,i)
     1822       ENDIF
     1823
    14551824
    14561825!
     
    14681837
    14691838          CALL adjust_cloud( i,j )
     1839          IF ( microphysics_morrison ) CALL activation( i,j )
     1840          IF ( microphysics_morrison ) CALL condensation( i,j )
    14701841          CALL autoconversion( i,j )
    14711842          CALL accretion( i,j )
     
    14831854       q(:,j,i)  = q_1d(:)
    14841855       pt(:,j,i) = pt_1d(:)
     1856       IF ( microphysics_morrison ) THEN
     1857          qc(:,j,i) = qc_1d(:)
     1858          nc(:,j,i) = nc_1d(:)
     1859       ENDIF
    14851860       IF ( microphysics_seifert )  THEN
    14861861          qr(:,j,i) = qr_1d(:)
     
    15031878       USE cloud_parameters,                                                   &
    15041879           ONLY:  hyrho
     1880
     1881       USE control_parameters,                                                 &
     1882           ONLY: microphysics_morrison
    15051883
    15061884       USE indices,                                                            &
     
    15381916          ENDIF
    15391917
     1918          IF ( microphysics_morrison ) THEN
     1919             IF ( qc_1d(k) <= eps_sb )  THEN
     1920                qc_1d(k) = 0.0_wp
     1921                nc_1d(k) = 0.0_wp
     1922             ELSE
     1923                IF ( nc_1d(k) * xcmin > qc_1d(k) * hyrho(k) )  THEN
     1924                   nc_1d(k) = qc_1d(k) * hyrho(k) / xamin * flag
     1925                ENDIF
     1926             ENDIF
     1927          ENDIF
     1928
    15401929       ENDDO
    15411930
    15421931    END SUBROUTINE adjust_cloud_ij
     1932
     1933!------------------------------------------------------------------------------!
     1934! Description:
     1935! ------------
     1936!> Calculate number of activated condensation nucleii after simple activation
     1937!> scheme of Twomey, 1959.
     1938!------------------------------------------------------------------------------!
     1939    SUBROUTINE activation_ij( i, j )
     1940
     1941       USE arrays_3d,                                                          &
     1942           ONLY:  hyp, nr, pt, q,  qc, qr, nc
     1943
     1944       USE cloud_parameters,                                                   &
     1945           ONLY:  hyrho, l_d_cp, l_d_r, l_v, rho_l, r_v, t_d_pt
     1946
     1947       USE constants,                                                          &
     1948           ONLY:  pi
     1949
     1950       USE cpulog,                                                             &
     1951           ONLY:  cpu_log, log_point_s
     1952
     1953       USE indices,                                                            &
     1954           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
     1955
     1956       USE kinds
     1957
     1958       USE control_parameters,                                                 &
     1959          ONLY: simulated_time
     1960
     1961       USE particle_attributes,                                                &
     1962          ONLY:  molecular_weight_of_solute, molecular_weight_of_water, rho_s, &
     1963              s1, s2, s3, vanthoff
     1964
     1965       IMPLICIT NONE
     1966
     1967       INTEGER(iwp) ::  i                 !<
     1968       INTEGER(iwp) ::  j                 !<
     1969       INTEGER(iwp) ::  k                 !<
     1970
     1971       REAL(wp)     ::  activ             !<
     1972       REAL(wp)     ::  afactor           !<
     1973       REAL(wp)     ::  alpha             !<
     1974       REAL(wp)     ::  beta_act          !<
     1975       REAL(wp)     ::  bfactor           !<
     1976       REAL(wp)     ::  e_s               !<
     1977       REAL(wp)     ::  k_act             !<
     1978       REAL(wp)     ::  n_act             !<
     1979       REAL(wp)     ::  n_ccn             !<
     1980       REAL(wp)     ::  q_s               !<
     1981       REAL(wp)     ::  rd0               !<
     1982       REAL(wp)     ::  s_0               !<
     1983       REAL(wp)     ::  sat               !<
     1984       REAL(wp)     ::  sat_max           !<
     1985       REAL(wp)     ::  sigma             !<
     1986       REAL(wp)     ::  t_int             !<
     1987       REAL(wp)     ::  t_l               !<
     1988
     1989       DO  k = nzb+1, nzt
     1990!
     1991!--       Actual liquid water temperature:
     1992          t_l = t_d_pt(k) * pt_1d(k)
     1993
     1994!
     1995!--       Calculate actual temperature
     1996          t_int = pt_1d(k) * ( hyp(k) / 100000.0_wp )**0.286_wp
     1997!
     1998!--             Saturation vapor pressure at t_l:
     1999          e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /             &
     2000                                 ( t_l - 35.86_wp )                            &
     2001                                 )
     2002!
     2003!--       Computation of saturation humidity:
     2004          q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
     2005          alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
     2006          q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) /                         &
     2007                        ( 1.0_wp + alpha * q_s )
     2008
     2009!--       Supersaturation:
     2010          sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
     2011
     2012!
     2013!--       Prescribe parameters for activation
     2014!--       (see: Bott + Trautmann, 2002, Atm. Res., 64)
     2015          k_act  = 0.7_wp
     2016
     2017          IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk )  THEN
     2018!
     2019!--          Compute the number of activated Aerosols
     2020!--          (see: Twomey, 1959, Pure and applied Geophysics, 43)
     2021             n_act     = na_init * sat**k_act
     2022!
     2023!--          Compute the number of cloud droplets
     2024!--          (see: Morrison + Grabowski, 2007, JAS, 64)
     2025!            activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro
     2026
     2027!
     2028!--          Compute activation rate after Khairoutdinov and Kogan
     2029!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)   
     2030             sat_max = 0.8_wp / 100.0_wp                 
     2031             activ   = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN              &
     2032                 ( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) /          &
     2033                  dt_micro
     2034
     2035             nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)
     2036          ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk )  THEN
     2037!
     2038!--          Curvature effect (afactor) with surface tension
     2039!--          parameterization by Straka (2009)
     2040             sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp )
     2041             afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int )
     2042!
     2043!--          Solute effect (bfactor)
     2044             bfactor = vanthoff * molecular_weight_of_water *                  &
     2045                  rho_s / ( molecular_weight_of_solute * rho_l )
     2046
     2047!
     2048!--          Prescribe power index that describes the soluble fraction
     2049!--          of an aerosol particle (beta) and mean geometric radius of
     2050!--          dry aerosol spectrum
     2051!--          (see: Morrison + Grabowski, 2007, JAS, 64)
     2052             beta_act = 0.5_wp
     2053             rd0      = 0.05E-6_wp
     2054!
     2055!--          Calculate mean geometric supersaturation (s_0) with
     2056!--          parameterization by Khvorostyanov and Curry (2006)
     2057             s_0   = rd0 **(- ( 1.0_wp + beta_act ) ) *                        &
     2058               ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp
     2059
     2060!
     2061!--          Calculate number of activated CCN as a function of
     2062!--          supersaturation and taking Koehler theory into account
     2063!--          (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)
     2064             n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF(                     &
     2065                LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(s1) ) ) )     
     2066             activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )
     2067                                     
     2068             nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init)                   
     2069          ENDIF
     2070
     2071       ENDDO
     2072
     2073    END SUBROUTINE activation_ij
     2074
     2075!------------------------------------------------------------------------------!
     2076! Description:
     2077! ------------
     2078!> Calculate condensation rate for cloud water content (after Khairoutdinov and 
     2079!> Kogan, 2000).
     2080!------------------------------------------------------------------------------!
     2081    SUBROUTINE condensation_ij( i, j )
     2082
     2083       USE arrays_3d,                                                          &
     2084           ONLY:  hyp, nr, pt, q,  qc, qr, nc
     2085
     2086       USE cloud_parameters,                                                   &
     2087           ONLY:  hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt
     2088
     2089       USE constants,                                                          &
     2090           ONLY:  pi
     2091
     2092       USE cpulog,                                                             &
     2093           ONLY:  cpu_log, log_point_s
     2094
     2095       USE indices,                                                            &
     2096           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
     2097
     2098       USE kinds
     2099
     2100       USE control_parameters,                                                 &
     2101          ONLY: simulated_time
     2102
     2103       IMPLICIT NONE
     2104
     2105       INTEGER(iwp) ::  i                 !<
     2106       INTEGER(iwp) ::  j                 !<
     2107       INTEGER(iwp) ::  k                 !<
     2108
     2109       REAL(wp)     ::  alpha             !<
     2110       REAL(wp)     ::  cond              !<
     2111       REAL(wp)     ::  cond_max          !<
     2112       REAL(wp)     ::  dc                !<
     2113       REAL(wp)     ::  e_s               !<
     2114       REAL(wp)     ::  evap              !<
     2115       REAL(wp)     ::  evap_nc           !<
     2116       REAL(wp)     ::  g_fac             !<
     2117       REAL(wp)     ::  nc_0              !<
     2118       REAL(wp)     ::  q_s               !<
     2119       REAL(wp)     ::  sat               !<
     2120       REAL(wp)     ::  t_l               !<
     2121       REAL(wp)     ::  temp              !<
     2122       REAL(wp)     ::  xc                !<
     2123
     2124
     2125       DO  k = nzb+1, nzt
     2126!
     2127!--       Actual liquid water temperature:
     2128          t_l = t_d_pt(k) * pt_1d(k)
     2129!
     2130!--             Saturation vapor pressure at t_l:
     2131          e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) /    &
     2132                                 ( t_l - 35.86_wp )                   &
     2133                                 )
     2134!
     2135!--       Computation of saturation humidity:
     2136          q_s   = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s )
     2137          alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l )
     2138          q_s   = q_s * ( 1.0_wp + alpha * q_1d(k) ) /               &
     2139                        ( 1.0_wp + alpha * q_s )
     2140
     2141!--       Supersaturation:
     2142          sat   = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp
     2143
     2144
     2145!
     2146!--       Actual temperature:
     2147          temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
     2148
     2149          g_fac  = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *        &
     2150                              l_v / ( thermal_conductivity_l * temp )    &
     2151                              + r_v * temp / ( diff_coeff_l * e_s )      &
     2152                            )
     2153!
     2154!--       Mean weight of cloud drops
     2155          IF ( nc_1d(k) <= 0.0_wp) CYCLE
     2156          xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin)               
     2157!
     2158!--       Weight averaged diameter of cloud drops:
     2159          dc   = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
     2160!
     2161!--       Integral diameter of cloud drops
     2162          nc_0 = nc_1d(k) * dc               
     2163!
     2164!--       Condensation needs only to be calculated in supersaturated regions
     2165          IF ( sat > 0.0_wp )  THEN
     2166!
     2167!--          Condensation rate of cloud water content
     2168!--          after KK scheme.
     2169!--          (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)
     2170             cond      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
     2171             cond_max  = q_1d(k) - q_s - qc_1d(k) - qr_1d(k)
     2172             cond      = MIN( cond, cond_max / dt_micro )
     2173               
     2174             qc_1d(k) = qc_1d(k) + cond * dt_micro
     2175          ELSEIF ( sat < 0.0_wp ) THEN
     2176             evap      = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)
     2177             evap      = MAX( evap, -qc_1d(k) / dt_micro )
     2178
     2179             qc_1d(k) = qc_1d(k) + evap * dt_micro
     2180          ENDIF
     2181       ENDDO
     2182
     2183    END SUBROUTINE condensation_ij
    15432184
    15442185
     
    15572198
    15582199       USE control_parameters,                                                 &
    1559            ONLY:  rho_surface
     2200           ONLY:  microphysics_morrison, rho_surface
    15602201
    15612202       USE grid_variables,                                                     &
     
    15792220       REAL(wp)     ::  k_au              !<
    15802221       REAL(wp)     ::  l_mix             !<
     2222       REAL(wp)     ::  nc_auto           !<
    15812223       REAL(wp)     ::  nu_c              !<
    15822224       REAL(wp)     ::  phi_au            !<
     
    15922234!--       Predetermine flag to mask topography
    15932235          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     2236          nc_auto = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
    15942237
    15952238          IF ( qc_1d(k) > eps_sb )  THEN
     
    16102253!
    16112254!--          Mean weight of cloud droplets:
    1612              xc = hyrho(k) * qc_1d(k) / nc_1d(k)
     2255             xc = hyrho(k) * qc_1d(k) / nc_auto
    16132256!
    16142257!--          Parameterized turbulence effects on autoconversion (Seifert,
     
    16582301             qc_1d(k) = qc_1d(k) - autocon * dt_micro                 * flag
    16592302             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag
     2303             IF ( microphysics_morrison ) THEN
     2304                nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp *                &
     2305                                    autocon / x0 * hyrho(k) * dt_micro * flag )
     2306             ENDIF
    16602307
    16612308          ENDIF
     
    17382385
    17392386       USE control_parameters,                                                 &
    1740            ONLY:  rho_surface
     2387           ONLY:  microphysics_morrison, rho_surface
    17412388
    17422389       USE indices,                                                            &
     
    17542401       REAL(wp)     ::  flag              !< flag to indicate first grid level above surface
    17552402       REAL(wp)     ::  k_cr              !<
     2403       REAL(wp)     ::  nc_accr           !<
    17562404       REAL(wp)     ::  phi_ac            !<
    17572405       REAL(wp)     ::  tau_cloud         !<
     2406       REAL(wp)     ::  xc                !<
     2407
    17582408
    17592409       DO  k = nzb+1, nzt
     
    17612411!--       Predetermine flag to mask topography
    17622412          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    1763 
    1764           IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb ) )  THEN
     2413          nc_accr = MERGE ( nc_1d(k), nc_const, microphysics_morrison )
     2414
     2415          IF ( ( qc_1d(k) > eps_sb )  .AND.  ( qr_1d(k) > eps_sb )  .AND.  &
     2416               ( nc_accr > eps_mr ) )  THEN
    17652417!
    17662418!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
     
    17702422!--          (Seifert and Beheng, 2001):
    17712423             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
     2424
     2425!
     2426!--          Mean weight of cloud drops
     2427             xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin)     
    17722428!
    17732429!--          Parameterized turbulence effects on autoconversion (Seifert,
     
    17842440!
    17852441!--          Accretion rate (Seifert and Beheng, 2006):
    1786              accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) )
     2442             accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac *                      &
     2443                    SQRT( rho_surface * hyrho(k) )
    17872444             accr = MIN( accr, qc_1d(k) / dt_micro )
    17882445
    17892446             qr_1d(k) = qr_1d(k) + accr * dt_micro * flag
    17902447             qc_1d(k) = qc_1d(k) - accr * dt_micro * flag
     2448             IF ( microphysics_morrison )  THEN
     2449                nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), accr / xc *               &
     2450                                             hyrho(k) * dt_micro * flag        &
     2451                                         )
     2452             ENDIF
     2453
    17912454
    17922455          ENDIF
     
    20182681
    20192682       USE control_parameters,                                                 &
    2020            ONLY:  call_microphysics_at_all_substeps, intermediate_timestep_count
     2683           ONLY:  call_microphysics_at_all_substeps,                           &
     2684                  intermediate_timestep_count, microphysics_morrison
    20212685
    20222686       USE indices,                                                            &
     
    20342698       INTEGER(iwp) ::  k             !<
    20352699
    2036        REAL(wp)                       ::  flag    !< flag to indicate first grid level above surface
     2700       REAL(wp)     ::  flag    !< flag to indicate first grid level above surface
     2701       REAL(wp)     ::  nc_sedi !<
     2702
     2703       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_nc  !<
    20372704       REAL(wp), DIMENSION(nzb:nzt+1) ::  sed_qc  !<
    20382705
     
    20432710!--       Predetermine flag to mask topography
    20442711          flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     2712          nc_sedi = MERGE( nc_1d(k), nc_const, microphysics_morrison )
     2713!
     2714!--       Sedimentation fluxes for number concentration are only calculated
     2715!--       for cloud_scheme = 'morrison'
     2716          IF ( microphysics_morrison ) THEN
     2717             IF ( qc_1d(k) > eps_sb  .AND.  nc_1d(k) > eps_mr )  THEN
     2718                sed_nc(k) = sed_qc_const *                                     &
     2719                            ( qc_1d(k) * hyrho(k) )**( 2.0_wp / 3.0_wp ) *     &
     2720                            ( nc_1d(k) )**( 1.0_wp / 3.0_wp )
     2721             ELSE
     2722                sed_nc(k) = 0.0_wp
     2723             ENDIF
     2724
     2725             sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) *                 &
     2726                              nc_1d(k) / dt_micro + sed_nc(k+1)                &
     2727                            ) * flag
     2728
     2729             nc_1d(k) = nc_1d(k) + ( sed_nc(k+1) - sed_nc(k) ) *               &
     2730                                         ddzu(k+1) / hyrho(k) * dt_micro * flag
     2731          ENDIF
    20452732
    20462733          IF ( qc_1d(k) > eps_sb )  THEN
    2047              sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) *       &
     2734             sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) *        &
    20482735                         ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )  * flag
    20492736          ELSE
  • palm/trunk/SOURCE/modules.f90

    r2277 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2277 2017-06-12 10:47:51Z kanani
    2732! Added doxygen comments for variables/parameters,
    2833! removed unused variables dissipation_control, do2d_xy_n, do2d_xz_n, do2d_yz_n,
     
    543548    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  c_w                   !< phase speed of w-velocity component
    544549    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_e              !< artificial numerical dissipation flux at south face of grid box - subgrid-scale TKE
     550    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_nc             !< artificial numerical dissipation flux at south face of grid box - clouddrop-number concentration   
    545551    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_nr             !< artificial numerical dissipation flux at south face of grid box - raindrop-number concentration   
    546552    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_pt             !< artificial numerical dissipation flux at south face of grid box - potential temperature
    547553    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_q              !< artificial numerical dissipation flux at south face of grid box - mixing ratio
     554    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qc             !< artificial numerical dissipation flux at south face of grid box - cloudwater
    548555    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_qr             !< artificial numerical dissipation flux at south face of grid box - rainwater
    549556    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  diss_s_s              !< artificial numerical dissipation flux at south face of grid box - passive scalar
     
    555562    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  dzw_mg                !< vertical grid size (w-grid) for multigrid pressure solver
    556563    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_e              !< 6th-order advective flux at south face of grid box - subgrid-scale TKE
     564    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_nc             !< 6th-order advective flux at south face of grid box - clouddrop-number concentration
    557565    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_nr             !< 6th-order advective flux at south face of grid box - raindrop-number concentration
    558566    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_pt             !< 6th-order advective flux at south face of grid box - potential temperature
    559567    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_q              !< 6th-order advective flux at south face of grid box - mixing ratio
     568    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_qc             !< 6th-order advective flux at south face of grid box - cloudwater
    560569    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_qr             !< 6th-order advective flux at south face of grid box - rainwater
    561570    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  flux_s_s              !< 6th-order advective flux at south face of grid box - passive scalar
     
    594603    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss       !< TKE dissipation rate
    595604    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_e   !< artificial numerical dissipation flux at left face of grid box - subgrid-scale TKE
     605    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_nc  !< artificial numerical dissipation flux at left face of grid box - clouddrop-number concentration
    596606    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_nr  !< artificial numerical dissipation flux at left face of grid box - raindrop-number concentration
    597607    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_pt  !< artificial numerical dissipation flux at left face of grid box - potential temperature
    598608    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_q   !< artificial numerical dissipation flux at left face of grid box - mixing ratio
     609    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_qc  !< artificial numerical dissipation flux at left face of grid box - cloudwater
    599610    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_qr  !< artificial numerical dissipation flux at left face of grid box - rainwater
    600611    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_s   !< artificial numerical dissipation flux at left face of grid box - passive scalar
     
    604615    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  diss_l_w   !< artificial numerical dissipation flux at left face of grid box - w-component
    605616    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_e   !< 6th-order advective flux at south face of grid box - subgrid-scale TKE
     617    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_nc  !< 6th-order advective flux at south face of grid box - clouddrop-number concentration
    606618    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_nr  !< 6th-order advective flux at south face of grid box - raindrop-number concentration
    607619    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_pt  !< 6th-order advective flux at south face of grid box - potential temperature
    608620    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_q   !< 6th-order advective flux at south face of grid box - mixing ratio
     621    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_qc  !< 6th-order advective flux at south face of grid box - cloudwater
    609622    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_qr  !< 6th-order advective flux at south face of grid box - rainwater
    610623    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  flux_l_s   !< 6th-order advective flux at south face of grid box - passive scalar
     
    636649    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  e          !< subgrid-scale turbulence kinetic energy (sgs tke)
    637650    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  e_p        !< prognostic value of sgs tke
     651    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc         !< cloud drop number density
     652    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_p       !< prognostic value of cloud drop number density
    638653    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr         !< rain drop number density
    639654    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_p       !< prognostic value of rain drop number density
     
    646661    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  q_p        !< prognostic value of specific humidity
    647662    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc         !< cloud water content
     663    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_p       !< cloud water content
    648664    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql         !< liquid water content
    649665    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_c       !< change in liquid water content due to
     
    659675    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  sa_p       !< prognostic value of ocean salinity
    660676    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  te_m       !< weighted tendency of e for previous sub-timestep (Runge-Kutta)
     677    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  tnc_m      !< weighted tendency of nc for previous sub-timestep (Runge-Kutta)
    661678    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  tnr_m      !< weighted tendency of nr for previous sub-timestep (Runge-Kutta)
    662679    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  tpt_m      !< weighted tendency of pt for previous sub-timestep (Runge-Kutta)
    663680    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  tq_m       !< weighted tendency of q for previous sub-timestep (Runge-Kutta)
     681    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  tqc_m      !< weighted tendency of qc for previous sub-timestep (Runge-Kutta)
    664682    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  tqr_m      !< weighted tendency of qr for previous sub-timestep (Runge-Kutta)
    665683    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ts_m       !< weighted tendency of s for previous sub-timestep (Runge-Kutta)
     
    681699    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  p       !< pointer: perturbation pressure
    682700    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  prho_1  !< pointer for swapping of timelevels for respective quantity
     701    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_1    !< pointer for swapping of timelevels for respective quantity
     702    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_2    !< pointer for swapping of timelevels for respective quantity
     703    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_3    !< pointer for swapping of timelevels for respective quantity
    683704    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_1    !< pointer for swapping of timelevels for respective quantity
    684705    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_2    !< pointer for swapping of timelevels for respective quantity
     
    691712    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  q_3     !< pointer for swapping of timelevels for respective quantity
    692713    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_1    !< pointer for swapping of timelevels for respective quantity
     714    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_2    !< pointer for swapping of timelevels for respective quantity
     715    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qc_3    !< pointer for swapping of timelevels for respective quantity
    693716    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_v    !< pointer: volume of liquid water
    694717    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ql_vp   !< pointer: liquid water weighting factor
     
    718741    REAL(wp), DIMENSION(:,:,:), POINTER ::  e          !< pointer: subgrid-scale turbulence kinetic energy (sgs tke)
    719742    REAL(wp), DIMENSION(:,:,:), POINTER ::  e_p        !< pointer: prognostic value of sgs tke
     743    REAL(wp), DIMENSION(:,:,:), POINTER ::  nc         !< pointer: cloud drop number density
     744    REAL(wp), DIMENSION(:,:,:), POINTER ::  nc_p       !< pointer: prognostic value of cloud drop number density
    720745    REAL(wp), DIMENSION(:,:,:), POINTER ::  nr         !< pointer: rain drop number density
    721746    REAL(wp), DIMENSION(:,:,:), POINTER ::  nr_p       !< pointer: prognostic value of rain drop number density
     
    726751    REAL(wp), DIMENSION(:,:,:), POINTER ::  q_p        !< pointer: prognostic value of specific humidity
    727752    REAL(wp), DIMENSION(:,:,:), POINTER ::  qc         !< pointer: cloud water content
     753    REAL(wp), DIMENSION(:,:,:), POINTER ::  qc_p       !< pointer: cloud water content
    728754    REAL(wp), DIMENSION(:,:,:), POINTER ::  ql         !< pointer: liquid water content
    729755    REAL(wp), DIMENSION(:,:,:), POINTER ::  ql_c       !< pointer: change in liquid water content due to
     
    737763    REAL(wp), DIMENSION(:,:,:), POINTER ::  sa_p       !< pointer: prognostic value of ocean salinity
    738764    REAL(wp), DIMENSION(:,:,:), POINTER ::  te_m       !< pointer: weighted tendency of e for previous sub-timestep (Runge-Kutta)
     765    REAL(wp), DIMENSION(:,:,:), POINTER ::  tnc_m      !< pointer: weighted tendency of nc for previous sub-timestep (Runge-Kutta)
    739766    REAL(wp), DIMENSION(:,:,:), POINTER ::  tnr_m      !< pointer: weighted tendency of nr for previous sub-timestep (Runge-Kutta)
    740767    REAL(wp), DIMENSION(:,:,:), POINTER ::  tpt_m      !< pointer: weighted tendency of pt for previous sub-timestep (Runge-Kutta)
    741768    REAL(wp), DIMENSION(:,:,:), POINTER ::  tq_m       !< pointer: weighted tendency of q for previous sub-timestep (Runge-Kutta)
     769    REAL(wp), DIMENSION(:,:,:), POINTER ::  tqc_m      !< pointer: weighted tendency of qc for previous sub-timestep (Runge-Kutta)
    742770    REAL(wp), DIMENSION(:,:,:), POINTER ::  tqr_m      !< pointer: weighted tendency of qr for previous sub-timestep (Runge-Kutta)
    743771    REAL(wp), DIMENSION(:,:,:), POINTER ::  ts_m       !< pointer: weighted tendency of s for previous sub-timestep (Runge-Kutta)
     
    800828    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  e_av          !< avg. subgrid-scale tke
    801829    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  lpt_av        !< avg. liquid water potential temperature
     830    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nc_av         !< avg. cloud drop number density
    802831    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nr_av         !< avg. rain drop number density
    803832    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  p_av          !< avg. perturbation pressure
     
    11441173    LOGICAL ::  microphysics_sat_adjust = .FALSE.            !< use saturation adjust bulk scheme?
    11451174    LOGICAL ::  microphysics_kessler = .FALSE.               !< use kessler bulk scheme?
    1146     LOGICAL ::  microphysics_seifert = .FALSE.               !< use 2-moment Seifert and Beheng scheme?
     1175    LOGICAL ::  microphysics_morrison = .FALSE.              !< use 2-moment Morrison (add. prog. eq. for nc and qc)
     1176    LOGICAL ::  microphysics_seifert = .FALSE.               !< use 2-moment Seifert and Beheng scheme
    11471177    LOGICAL ::  mg_switch_to_pe0 = .FALSE.                   !< (Siggi add short description)
    11481178    LOGICAL ::  nest_bound_l = .FALSE.                       !< nested boundary on left side?
     
    19802010    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_vs2_ws_l    !< subdomain sum of horizontal momentum flux v'v' (5th-order advection scheme only)
    19812011    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_ws2_ws_l    !< subdomain sum of vertical momentum flux w'w' (5th-order advection scheme only)
     2012    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsncs_ws_l  !< subdomain sum of vertical clouddrop-number concentration flux w'nc' (5th-order advection scheme only)
    19822013    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsnrs_ws_l  !< subdomain sum of vertical raindrop-number concentration flux w'nr' (5th-order advection scheme only)
    19832014    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wspts_ws_l  !< subdomain sum of vertical sensible heat flux w'pt' (5th-order advection scheme only)
    19842015    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqs_ws_l   !< subdomain sum of vertical latent heat flux w'q' (5th-order advection scheme only)
     2016    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqcs_ws_l  !< subdomain sum of vertical cloudwater flux w'qc' (5th-order advection scheme only)
    19852017    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wsqrs_ws_l  !< subdomain sum of vertical rainwater flux w'qr' (5th-order advection scheme only)
    19862018    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sums_wssas_ws_l  !< subdomain sum of vertical salinity flux w'sa' (5th-order advection scheme only)
  • palm/trunk/SOURCE/netcdf_interface_mod.f90

    r2270 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2270 2017-06-09 12:18:47Z maronga
    2732! Removed 2 timeseries (shf_eb + qsws_eb). Removed _eb suffixes
    2833!
     
    844849!
    845850!--             Most variables are defined on the scalar grid
    846                 CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
    847                        'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',        &
     851                CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
     852                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
    848853                       'rho_ocean', 's', 'sa', 'vpt'  )
    849854
     
    13761381!
    13771382!--             Most variables are defined on the scalar grid
    1378                 CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
    1379                        'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv', 'rho_ocean', &
    1380                        's', 'sa', 'vpt' )
     1383                CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
     1384                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',  &
     1385                       'rho_ocean', 's', 'sa', 'vpt' )
    13811386
    13821387                   grid_x = 'x'
     
    20272032!
    20282033!--                   Most variables are defined on the zu grid
    2029                       CASE ( 'e_xy', 'lpt_xy', 'nr_xy', 'p_xy', 'pc_xy',       &
    2030                              'pr_xy', 'prr_xy', 'pt_xy', 'q_xy', 'qc_xy',      &
    2031                              'ql_xy', 'ql_c_xy', 'ql_v_xy', 'ql_vp_xy',        &
    2032                              'qr_xy', 'qv_xy', 'rho_ocean_xy', 's_xy',         &
    2033                              'sa_xy', 'vpt_xy' )
     2034                      CASE ( 'e_xy', 'lpt_xy', 'nc_xy', 'nr_xy', 'p_xy',       &
     2035                             'pc_xy', 'pr_xy', 'prr_xy', 'pt_xy', 'q_xy',      &
     2036                             'qc_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy',           &
     2037                             'ql_vp_xy', 'qr_xy', 'qv_xy', 'rho_ocean_xy',     &
     2038                             's_xy', 'sa_xy', 'vpt_xy' )
    20342039
    20352040                         grid_x = 'x'
     
    27262731!
    27272732!--                Most variables are defined on the zu grid
    2728                    CASE ( 'e_xz', 'lpt_xz', 'nr_xz', 'p_xz', 'pc_xz', 'pr_xz', &
    2729                           'prr_xz', 'pt_xz', 'q_xz', 'qc_xz', 'ql_xz',         &
    2730                           'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz', 'qv_xz',  &
    2731                           'rho_ocean_xz', 's_xz', 'sa_xz', 'vpt_xz' )
     2733                   CASE ( 'e_xz', 'lpt_xz', 'nc_xz', 'nr_xz', 'p_xz', 'pc_xz', &
     2734                          'pr_xz', 'prr_xz', 'pt_xz', 'q_xz', 'qc_xz',         &
     2735                          'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz',  &
     2736                          'qv_xz', 'rho_ocean_xz', 's_xz', 'sa_xz', 'vpt_xz' )
    27322737
    27332738                      grid_x = 'x'
     
    33803385!
    33813386!--                Most variables are defined on the zu grid
    3382                    CASE ( 'e_yz', 'lpt_yz', 'nr_yz', 'p_yz', 'pc_yz', 'pr_yz', &
    3383                           'prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz',        &
     3387                   CASE ( 'e_yz', 'lpt_yz', 'nc_yz', 'nr_yz', 'p_yz', 'pc_yz', &
     3388                          'pr_yz','prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz', &
    33843389                          'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz',  &
    33853390                          'rho_ocean_yz', 's_yz', 'sa_yz', 'vpt_yz' )
     
    45574562!
    45584563!--             Most variables are defined on the zu levels
    4559                 CASE ( 'e', 'lpt', 'nr', 'p', 'pc', 'pr', 'prr', 'pt', 'q',    &
    4560                        'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',        &
     4564                CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
     4565                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
    45614566                       'rho_ocean', 's', 'sa', 'u', 'v', 'vpt' )
    45624567
  • palm/trunk/SOURCE/palm.f90

    r2261 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2261 2017-06-08 14:25:57Z raasch
    2732! output of run number for mrun to create unified cycle numbers
    2833!
     
    177182               do2d_at_begin, do3d_at_begin, humidity, initializing_actions,   &
    178183               io_blocks, io_group, land_surface, large_scale_forcing,         &
    179                message_string, microphysics_seifert, nest_domain, neutral,     &
    180                nudging, passive_scalar, runnr, simulated_time,                 &
    181                simulated_time_chr, urban_surface,                              &
     184               message_string, microphysics_morrison, microphysics_seifert,    &
     185               nest_domain, neutral, nudging, passive_scalar, runnr,           &
     186               simulated_time, simulated_time_chr, urban_surface,              &
    182187               user_interface_current_revision,                                &
    183188               user_interface_required_revision, version, wall_heatflux,       &
     
    364369             IF ( humidity )  THEN
    365370                CALL exchange_horiz( q, nbgp )
     371                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     372                  CALL exchange_horiz( qc, nbgp )
     373                  CALL exchange_horiz( nc, nbgp )
     374                ENDIF
    366375                IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    367 !                   CALL exchange_horiz( qc, nbgp )
    368376                   CALL exchange_horiz( qr, nbgp )
    369 !                   CALL exchange_horiz( nc, nbgp )
    370377                   CALL exchange_horiz( nr, nbgp )
    371378                ENDIF
  • palm/trunk/SOURCE/parin.f90

    r2267 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2267 2017-06-09 09:33:25Z gronemeier
    2732! Bugfix: removed skipping of reading namelists in case of omitted d3par
    2833!
     
    292297    USE microphysics_mod,                                                      &
    293298        ONLY:  c_sedimentation, cloud_water_sedimentation,                     &
    294                collision_turbulence, limiter_sedimentation, nc_const,          &
    295                ventilation_effect
     299               collision_turbulence, curvature_solution_effects_bulk,          &
     300               limiter_sedimentation ,na_init, nc_const, ventilation_effect
    296301
    297302    USE model_1d,                                                              &
     
    350355             collective_wait, collision_turbulence, conserve_volume_flow,      &
    351356             conserve_volume_flow_mode, constant_flux_layer,                   &
    352              coupling_start_time,                                              &
     357             coupling_start_time, curvature_solution_effects_bulk,             &
    353358             cycle_mg, damp_level_1d,                                          &
    354359             dissipation_1d,                                                   &
     
    366371             loop_optimization, lsf_exception, masking_method, mg_cycles,      &
    367372             mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
    368              most_method, nc_const, netcdf_precision, neutral, ngsrb,          &
     373             most_method, na_init, nc_const, netcdf_precision, neutral, ngsrb, &
    369374             nsor, nsor_ini, nudging, nx, ny, nz, ocean, omega, omega_sor,     &
    370375             outflow_source_plane, passive_scalar, phi,                        &
  • palm/trunk/SOURCE/pmc_interface_mod.f90

    r2285 r2292  
    2626! -----------------
    2727! $Id$
     28! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     29! includes two more prognostic equations for cloud drop concentration (nc) 
     30! and cloud water content (qc).
     31!
     32! 2285 2017-06-15 13:15:41Z suehring
    2833! Consider multiple pure-vertical nest domains in overlap check
    2934!
     
    154159#if defined( __nopointer )
    155160    USE arrays_3d,                                                             &
    156         ONLY:  dzu, dzw, e, e_p, nr, pt, pt_p, q, q_p, qr, u, u_p, v, v_p,     &
    157                w, w_p, zu, zw
     161        ONLY:  dzu, dzw, e, e_p, nc, nr, pt, pt_p, q, q_p, qc, qr, u, u_p, v,  &
     162               v_p, w, w_p, zu, zw
    158163#else
    159164   USE arrays_3d,                                                              &
    160         ONLY:  dzu, dzw, e, e_p, e_1, e_2, nr, nr_2, nr_p, pt, pt_p, pt_1,    &
    161                pt_2, q, q_p, q_1, q_2, qr, qr_2, s, s_2, u, u_p, u_1, u_2, v, &
    162                v_p, v_1, v_2, w, w_p, w_1, w_2, zu, zw
     165        ONLY:  dzu, dzw, e, e_p, e_1, e_2, nc, nc_2, nc_p, nr, nr_2, nr_p, pt, &
     166               pt_p, pt_1, pt_2, q, q_p, q_1, q_2, qc, qc_2, qr, qr_2, s, s_2, &
     167               u, u_p, u_1, u_2, v, v_p, v_1, v_2, w, w_p, w_1, w_2, zu, zw
    163168#endif
    164169
    165170    USE control_parameters,                                                     &
    166171        ONLY:  cloud_physics, coupling_char, dt_3d, dz, humidity,               &
    167                message_string, microphysics_seifert, nest_bound_l, nest_bound_r,&
    168                nest_bound_s, nest_bound_n, nest_domain, neutral, passive_scalar,&
    169                roughness_length, simulated_time, topography, volume_flow
     172               message_string, microphysics_morrison, microphysics_seifert,     &
     173               nest_bound_l, nest_bound_r, nest_bound_s, nest_bound_n,          &
     174               nest_domain, neutral, passive_scalar, roughness_length,          &
     175               simulated_time, topography, volume_flow
    170176
    171177    USE cpulog,                                                                 &
     
    281287    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wc   !:
    282288    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  q_c  !:
    283 !     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qcc  !:
     289    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qcc  !:
    284290    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  qrc  !:
    285291    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  nrc  !:
    286 !     REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ncc  !:
     292    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  ncc  !:
    287293    REAL(wp), SAVE, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  sc   !:
    288294
     
    966972          CALL pmc_set_dataarray_name( 'coarse', 'q'  ,'fine', 'q',  ierr )
    967973
     974          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     975            CALL pmc_set_dataarray_name( 'coarse', 'qc'  ,'fine', 'qc',  ierr ) 
     976            CALL pmc_set_dataarray_name( 'coarse', 'nc'  ,'fine', 'nc',  ierr )
     977          ENDIF
     978
    968979          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    969 !              CALL pmc_set_dataarray_name( 'coarse', 'qc'  ,'fine', 'qc',  ierr ) 
    970980             CALL pmc_set_dataarray_name( 'coarse', 'qr'  ,'fine', 'qr',  ierr )
    971 !             CALL pmc_set_dataarray_name( 'coarse', 'nc'  ,'fine', 'nc',  ierr )
    972981             CALL pmc_set_dataarray_name( 'coarse', 'nr'  ,'fine', 'nr',  ierr )
    973 
    974982          ENDIF
    975983     
     
    29902998    IF ( TRIM(name) == "pt" )  p_3d => pt
    29912999    IF ( TRIM(name) == "q"  )  p_3d => q
    2992 !     IF ( TRIM(name) == "qc" )  p_3d => qc
     3000    IF ( TRIM(name) == "qc" )  p_3d => qc
    29933001    IF ( TRIM(name) == "qr" )  p_3d => qr
    29943002    IF ( TRIM(name) == "nr" )  p_3d => nr
    2995 !     IF ( TRIM(name) == "nc" )  p_3d => nc
     3003    IF ( TRIM(name) == "nc" )  p_3d => nc
    29963004    IF ( TRIM(name) == "s"  )  p_3d => s
    29973005!
     
    30263034    IF ( TRIM(name) == "pt" )  p_3d_sec => pt_2
    30273035    IF ( TRIM(name) == "q"  )  p_3d_sec => q_2
    3028 !     IF ( TRIM(name) == "qc" )  p_3d_sec => qc_2
     3036    IF ( TRIM(name) == "qc" )  p_3d_sec => qc_2
    30293037    IF ( TRIM(name) == "qr" )  p_3d_sec => qr_2
    30303038    IF ( TRIM(name) == "nr" )  p_3d_sec => nr_2
    3031 !     IF ( TRIM(name) == "nc" )  p_3d_sec => nc_2
     3039    IF ( TRIM(name) == "nc" )  p_3d_sec => nc_2
    30323040    IF ( TRIM(name) == "s"  )  p_3d_sec => s_2
    30333041
     
    31023110       IF ( .NOT. ALLOCATED( q_c ) ) ALLOCATE( q_c(0:nzc+1, js:je, is:ie) )
    31033111       p_3d => q_c
    3104 !     ELSEIF ( TRIM( name ) == "qc")  THEN
    3105 !        IF ( .NOT. ALLOCATED( qcc ) ) ALLOCATE( qcc(0:nzc+1, js:je, is:ie) )
    3106 !        p_3d => qcc
     3112    ELSEIF ( TRIM( name ) == "qc")  THEN
     3113       IF ( .NOT. ALLOCATED( qcc ) ) ALLOCATE( qcc(0:nzc+1, js:je, is:ie) )
     3114       p_3d => qcc
    31073115    ELSEIF ( TRIM( name ) == "qr")  THEN
    31083116       IF ( .NOT. ALLOCATED( qrc ) ) ALLOCATE( qrc(0:nzc+1, js:je, is:ie) )
     
    31113119       IF ( .NOT. ALLOCATED( nrc ) ) ALLOCATE( nrc(0:nzc+1, js:je, is:ie) )
    31123120       p_3d => nrc
    3113 !     ELSEIF ( TRIM( name ) == "nc")  THEN
    3114 !        IF ( .NOT. ALLOCATED( ncc ) ) ALLOCATE( ncc(0:nzc+1, js:je, is:ie) )
    3115 !        p_3d => ncc
     3121    ELSEIF ( TRIM( name ) == "nc")  THEN
     3122       IF ( .NOT. ALLOCATED( ncc ) ) ALLOCATE( ncc(0:nzc+1, js:je, is:ie) )
     3123       p_3d => ncc
    31163124    ELSEIF ( TRIM( name ) == "s")  THEN
    31173125       IF ( .NOT. ALLOCATED( sc ) ) ALLOCATE( sc(0:nzc+1, js:je, is:ie) )
     
    32243232                                      r2yo, r1zo, r2zo, 's' )
    32253233
     3234          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     3235             CALL pmci_interp_tril_all ( qc, qcc, ico, jco, kco, r1xo, r2xo,    &
     3236                                          r1yo, r2yo, r1zo, r2zo, 's' )
     3237             CALL pmci_interp_tril_all ( nc, ncc, ico, jco, kco, r1xo, r2xo,    &
     3238                                         r1yo, r2yo, r1zo, r2zo, 's' )   
     3239          ENDIF
     3240
    32263241          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    3227 !              CALL pmci_interp_tril_all ( qc, qcc, ico, jco, kco, r1xo, r2xo,    &
    3228 !                                          r1yo, r2yo, r1zo, r2zo, 's' )
    32293242             CALL pmci_interp_tril_all ( qr, qrc, ico, jco, kco, r1xo, r2xo,    &
    32303243                                         r1yo, r2yo, r1zo, r2zo, 's' )
    3231 !             CALL pmci_interp_tril_all ( nc, ncc, ico, jco, kco, r1xo, r2xo,    &
    3232 !                                         r1yo, r2yo, r1zo, r2zo, 's' )   
    32333244             CALL pmci_interp_tril_all ( nr, nrc, ico, jco, kco, r1xo, r2xo,    &
    32343245                                         r1yo, r2yo, r1zo, r2zo, 's' )
     
    39553966                                          nzt_topo_nestbc_l, 'l', 's' )
    39563967
     3968                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     3969                   CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo, r2xo,&
     3970                                             r1yo, r2yo, r1zo, r2zo,            &
     3971                                             logc_u_l,                          &
     3972                                             logc_ratio_u_l, nzt_topo_nestbc_l, &
     3973                                             'l', 's' ) 
     3974
     3975                   CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo, r2xo,&
     3976                                             r1yo, r2yo, r1zo, r2zo,            &
     3977                                             logc_u_l,                          &
     3978                                             logc_ratio_u_l, nzt_topo_nestbc_l, &
     3979                                             'l', 's' )         
     3980                ENDIF
    39573981
    39583982                IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    3959 !                    CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo, r2xo,&
    3960 !                                              r1yo, r2yo, r1zo, r2zo,            &
    3961 !                                              logc_u_l,                          &
    3962 !                                              logc_ratio_u_l, nzt_topo_nestbc_l, &
    3963 !                                              'l', 's' ) 
    3964 
    39653983                   CALL pmci_interp_tril_lr( qr, qrc, ico, jco, kco, r1xo, r2xo,&
    39663984                                             r1yo, r2yo, r1zo, r2zo,            &
     
    39683986                                             logc_ratio_u_l, nzt_topo_nestbc_l, &
    39693987                                             'l', 's' )
    3970 
    3971 !                    CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo, r2xo,&
    3972 !                                              r1yo, r2yo, r1zo, r2zo,            &
    3973 !                                              logc_u_l,                          &
    3974 !                                              logc_ratio_u_l, nzt_topo_nestbc_l, &
    3975 !                                              'l', 's' )
    39763988
    39773989                   CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo, r2xo,&
     
    40044016                   CALL pmci_extrap_ifoutflow_lr( q, 'l', 's' )
    40054017
     4018                   IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     4019
     4020                      CALL pmci_extrap_ifoutflow_lr( qc, 'l', 's' )
     4021                      CALL pmci_extrap_ifoutflow_lr( nc, 'l', 's' )
     4022
     4023                   ENDIF
     4024
    40064025                   IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    40074026
    4008 !                       CALL pmci_extrap_ifoutflow_lr( qc, 'l', 's' )
    40094027                      CALL pmci_extrap_ifoutflow_lr( qr, 'l', 's' )
    4010 !                      CALL pmci_extrap_ifoutflow_lr( nc, 'l', 's' )
    40114028                      CALL pmci_extrap_ifoutflow_lr( nr, 'l', 's' )
    40124029
     
    40614078                                          nzt_topo_nestbc_r, 'r', 's' )
    40624079
     4080                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     4081
     4082                   CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo,     &
     4083                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
     4084                                             logc_u_r,                         &
     4085                                             logc_ratio_u_r, nzt_topo_nestbc_r,&
     4086                                             'r', 's' )
     4087     
     4088                   CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo,     &
     4089                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
     4090                                             logc_u_r,                         &
     4091                                             logc_ratio_u_r, nzt_topo_nestbc_r,&
     4092                                             'r', 's' )
     4093
     4094
     4095                ENDIF
    40634096
    40644097                IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    40654098
    4066 !                    CALL pmci_interp_tril_lr( qc, qcc, ico, jco, kco, r1xo,     &
    4067 !                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
    4068 !                                              logc_u_r,                         &
    4069 !                                              logc_ratio_u_r, nzt_topo_nestbc_r,&
    4070 !                                              'r', 's' )
    40714099     
    40724100                   CALL pmci_interp_tril_lr( qr, qrc, ico, jco, kco, r1xo,     &
     
    40764104                                             'r', 's' )
    40774105
    4078 !                    CALL pmci_interp_tril_lr( nc, ncc, ico, jco, kco, r1xo,     &
    4079 !                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
    4080 !                                              logc_u_r,                         &
    4081 !                                              logc_ratio_u_r, nzt_topo_nestbc_r,&
    4082 !                                              'r', 's' )
    4083 
    40844106                   CALL pmci_interp_tril_lr( nr, nrc, ico, jco, kco, r1xo,     &
    40854107                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
     
    41134135                   CALL pmci_extrap_ifoutflow_lr( q, 'r', 's' )
    41144136
     4137                   IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     4138                      CALL pmci_extrap_ifoutflow_lr( qc, 'r', 's' )
     4139                      CALL pmci_extrap_ifoutflow_lr( nc, 'r', 's' ) 
     4140                   ENDIF
     4141
    41154142                   IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    4116 !                       CALL pmci_extrap_ifoutflow_lr( qc, 'r', 's' )
    41174143                      CALL pmci_extrap_ifoutflow_lr( qr, 'r', 's' )
    4118 !                      CALL pmci_extrap_ifoutflow_lr( nc, 'r', 's' ) 
    41194144                      CALL pmci_extrap_ifoutflow_lr( nr, 'r', 's' )
    41204145                   ENDIF
     
    41624187                                          nzt_topo_nestbc_s, 's', 's' )
    41634188
     4189                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     4190
     4191                   CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
     4192                                             r2xo, r1yo,r2yo, r1zo, r2zo,      &
     4193                                             logc_u_s,                         &
     4194                                             logc_ratio_u_s, nzt_topo_nestbc_s,&
     4195                                             's', 's' )
     4196
     4197                   CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
     4198                                             r2xo, r1yo,r2yo, r1zo, r2zo,      &
     4199                                             logc_u_s,                         &
     4200                                             logc_ratio_u_s, nzt_topo_nestbc_s,&
     4201                                             's', 's' )
     4202
     4203                ENDIF
     4204
    41644205                IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    4165 
    4166 !                    CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
    4167 !                                              r2xo, r1yo,r2yo, r1zo, r2zo,      &
    4168 !                                              logc_u_s,                         &
    4169 !                                              logc_ratio_u_s, nzt_topo_nestbc_s,&
    4170 !                                              's', 's' )
    41714206
    41724207                   CALL pmci_interp_tril_sn( qr, qrc, ico, jco, kco, r1xo,     &
     
    41754210                                             logc_ratio_u_s, nzt_topo_nestbc_s,&
    41764211                                             's', 's' )
    4177 
    4178 !                    CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
    4179 !                                              r2xo, r1yo,r2yo, r1zo, r2zo,      &
    4180 !                                              logc_u_s,                         &
    4181 !                                              logc_ratio_u_s, nzt_topo_nestbc_s,&
    4182 !                                              's', 's' )
    41834212
    41844213                   CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo,     &
     
    42124241                   CALL pmci_extrap_ifoutflow_sn( q,  's', 's' )
    42134242
     4243                   IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     4244                      CALL pmci_extrap_ifoutflow_sn( qc, 's', 's' )
     4245                      CALL pmci_extrap_ifoutflow_sn( nc, 's', 's' )
     4246                   ENDIF
     4247
    42144248                   IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    4215 !                       CALL pmci_extrap_ifoutflow_sn( qc, 's', 's' )
    42164249                      CALL pmci_extrap_ifoutflow_sn( qr, 's', 's' )     
    4217 !                       CALL pmci_extrap_ifoutflow_sn( nc, 's', 's' )
    42184250                      CALL pmci_extrap_ifoutflow_sn( nr, 's', 's' )
    42194251
     
    42664298                                          nzt_topo_nestbc_n, 'n', 's' )
    42674299
     4300                IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     4301
     4302                   CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
     4303                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
     4304                                             logc_u_n,                         &
     4305                                             logc_ratio_u_n, nzt_topo_nestbc_n,&
     4306                                             'n', 's' )
     4307
     4308                   CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
     4309                                             r2xo, r1yo, r2yo, r1zo, r2zo,     &
     4310                                             logc_u_n,                         &
     4311                                             logc_ratio_u_n, nzt_topo_nestbc_n,&
     4312                                             'n', 's' )
     4313
     4314                ENDIF
     4315
    42684316                IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    4269 
    4270 !                    CALL pmci_interp_tril_sn( qc, qcc, ico, jco, kco, r1xo,     &
    4271 !                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
    4272 !                                              logc_u_n,                         &
    4273 !                                              logc_ratio_u_n, nzt_topo_nestbc_n,&
    4274 !                                              'n', 's' )
    42754317
    42764318                   CALL pmci_interp_tril_sn( qr, qrc, ico, jco, kco, r1xo,     &
     
    42794321                                             logc_ratio_u_n, nzt_topo_nestbc_n,&
    42804322                                             'n', 's' )
    4281 
    4282 !                    CALL pmci_interp_tril_sn( nc, ncc, ico, jco, kco, r1xo,     &
    4283 !                                              r2xo, r1yo, r2yo, r1zo, r2zo,     &
    4284 !                                              logc_u_n,                         &
    4285 !                                              logc_ratio_u_n, nzt_topo_nestbc_n,&
    4286 !                                              'n', 's' )
    42874323
    42884324                   CALL pmci_interp_tril_sn( nr, nrc, ico, jco, kco, r1xo,     &
     
    43174353
    43184354                   IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    4319 !                       CALL pmci_extrap_ifoutflow_sn( qc, 'n', 's' )
     4355                      CALL pmci_extrap_ifoutflow_sn( qc, 'n', 's' )
     4356                      CALL pmci_extrap_ifoutflow_sn( nc, 'n', 's' )
     4357                   ENDIF
     4358
     4359                   IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    43204360                      CALL pmci_extrap_ifoutflow_sn( qr, 'n', 's' )
    4321 !                       CALL pmci_extrap_ifoutflow_sn( nc, 'n', 's' )
    43224361                      CALL pmci_extrap_ifoutflow_sn( nr, 'n', 's' )
    43234362                   ENDIF
     
    43564395                                   r2yo, r1zo, r2zo, 's' )
    43574396
     4397          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     4398
     4399             CALL pmci_interp_tril_t( qc, qcc, ico, jco, kco, r1xo, r2xo, r1yo,&
     4400                                      r2yo, r1zo, r2zo, 's' )
     4401
     4402             CALL pmci_interp_tril_t( nc, ncc, ico, jco, kco, r1xo, r2xo, r1yo,&
     4403                                      r2yo, r1zo, r2zo, 's' )
     4404
     4405          ENDIF
     4406
    43584407          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    43594408
    4360 !              CALL pmci_interp_tril_t( qc, qcc, ico, jco, kco, r1xo, r2xo, r1yo,&
    4361 !                                       r2yo, r1zo, r2zo, 's' )
    43624409
    43634410             CALL pmci_interp_tril_t( qr, qrc, ico, jco, kco, r1xo, r2xo, r1yo,&
    43644411                                      r2yo, r1zo, r2zo, 's' )
    4365 
    4366 !              CALL pmci_interp_tril_t( nc, ncc, ico, jco, kco, r1xo, r2xo, r1yo,&
    4367 !                                       r2yo, r1zo, r2zo, 's' )
    43684412
    43694413             CALL pmci_interp_tril_t( nr, nrc, ico, jco, kco, r1xo, r2xo, r1yo,&
     
    43944438             CALL pmci_extrap_ifoutflow_t( q, 's' )
    43954439
     4440             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     4441                CALL pmci_extrap_ifoutflow_t( qc, 's' )
     4442                CALL pmci_extrap_ifoutflow_t( nc, 's' )
     4443             ENDIF
     4444
    43964445             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    4397 !                 CALL pmci_extrap_ifoutflow_t( qc, 's' )
    43984446                CALL pmci_extrap_ifoutflow_t( qr, 's' )
    4399 !                 CALL pmci_extrap_ifoutflow_t( nc, 's' )
    44004447                CALL pmci_extrap_ifoutflow_t( nr, 's' )
    44014448
     
    44384485                                  kfuo, ijfc_s, kfc_s, 's' )
    44394486
     4487         IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     4488
     4489            CALL pmci_anterp_tophat( qc, qcc, kctu, iflo, ifuo, jflo, jfuo,    &
     4490                                     kflo, kfuo, ijfc_s, kfc_s, 's' )
     4491
     4492            CALL pmci_anterp_tophat( nc, ncc, kctu, iflo, ifuo, jflo, jfuo,    &
     4493                                     kflo, kfuo, ijfc_s, kfc_s,  's' )
     4494
     4495         ENDIF
     4496
    44404497         IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    4441 
    4442 !             CALL pmci_anterp_tophat( qc, qcc, kctu, iflo, ifuo, jflo, jfuo,    &
    4443 !                                      kflo, kfuo, ijfc_s, kfc_s, 's' )
    44444498
    44454499            CALL pmci_anterp_tophat( qr, qrc, kctu, iflo, ifuo, jflo, jfuo,    &
    44464500                                     kflo, kfuo, ijfc_s, kfc_s, 's' )
    4447 
    4448 !             CALL pmci_anterp_tophat( nc, ncc, kctu, iflo, ifuo, jflo, jfuo,    &
    4449 !                                      kflo, kfuo, ijfc_s, kfc_s,  's' )
    44504501
    44514502            CALL pmci_anterp_tophat( nr, nrc, kctu, iflo, ifuo, jflo, jfuo,    &
  • palm/trunk/SOURCE/prognostic_equations.f90

    r2261 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2261 2017-06-08 14:25:57Z raasch
    2732! bugfix for r2232: openmp directives removed
    2833!
     
    227232
    228233    USE arrays_3d,                                                             &
    229         ONLY:  diss_l_e, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr,            &
    230                diss_l_s, diss_l_sa, diss_s_e, diss_s_nr, diss_s_pt, diss_s_q,  &
    231                diss_s_qr, diss_s_s, diss_s_sa, e, e_p, flux_s_e, flux_s_nr,    &
    232                flux_s_pt, flux_s_q, flux_s_qr, flux_s_s, flux_s_sa, flux_l_e,  &
    233                flux_l_nr, flux_l_pt, flux_l_q, flux_l_qr, flux_l_s, flux_l_sa, &
    234                nr, nr_p, pt, ptdf_x, ptdf_y, pt_init, pt_p,                    &
    235                prho, q, q_init, q_p, qr, qr_p, rdf,                            &
    236                rdf_sc, ref_state, rho_ocean, s, s_init, s_p, sa, sa_init, sa_p,&
    237                tend, te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m, tsa_m, tu_m, tv_m, &
    238                tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
     234        ONLY:  diss_l_e, diss_l_nc, diss_l_nr, diss_l_pt, diss_l_q, diss_l_qc, &
     235               diss_l_qr, diss_l_s, diss_l_sa, diss_s_e, diss_s_nc, diss_s_nr, &
     236               diss_s_pt, diss_s_q, diss_s_qc, diss_s_qr, diss_s_s, diss_s_sa, &
     237               e, e_p, flux_s_e, flux_s_nc, flux_s_nr, flux_s_pt, flux_s_q,    &
     238               flux_s_qc, flux_s_qr, flux_s_s, flux_s_sa, flux_l_e, flux_l_nc, &
     239               flux_l_nr, flux_l_pt, flux_l_q, flux_l_qc, flux_l_qr, flux_l_s, &
     240               flux_l_sa, nc, nc_p, nr, nr_p, pt, ptdf_x, ptdf_y, pt_init,     &
     241               pt_p, prho, q, q_init, q_p, qc, qc_p, qr, qr_p, rdf, rdf_sc,    &
     242               ref_state, rho_ocean, s,  s_init, s_p, sa, sa_init, sa_p, tend, &
     243               te_m, tnc_m,  tnr_m, tpt_m, tq_m, tqc_m, tqr_m, ts_m, tsa_m,    &
     244               tu_m, tv_m, tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p,  &
     245               w, w_p
    239246
    240247    USE control_parameters,                                                    &
     
    244251               inflow_l, intermediate_timestep_count,                          &
    245252               intermediate_timestep_count_max, large_scale_forcing,           &
    246                large_scale_subsidence, microphysics_seifert,                  &
     253               large_scale_subsidence, microphysics_morrison, microphysics_seifert, &
    247254               microphysics_sat_adjust, neutral, nudging, ocean, outflow_l,    &
    248255               outflow_s, passive_scalar, prho_reference, prho_reference,      &
     
    892899
    893900!
     901!--          If required, calculate prognostic equations for cloud water content
     902!--          and cloud drop concentration
     903             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     904!
     905!--             Calculate prognostic equation for cloud water content
     906                tend(:,j,i) = 0.0_wp
     907                IF ( timestep_scheme(1:5) == 'runge' ) &
     908                THEN
     909                   IF ( ws_scheme_sca )  THEN
     910                      CALL advec_s_ws( i, j, qc, 'qc', flux_s_qc,       &
     911                                       diss_s_qc, flux_l_qc, diss_l_qc, &
     912                                       i_omp_start, tn )
     913                   ELSE
     914                      CALL advec_s_pw( i, j, qc )
     915                   ENDIF
     916                ELSE
     917                   CALL advec_s_up( i, j, qc )
     918                ENDIF
     919                CALL diffusion_s( i, j, qc,                                   &
     920                                  surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,   &
     921                                  surf_def_h(2)%qcsws,                        &
     922                                  surf_lsm_h%qcsws,    surf_usm_h%qcsws,      & 
     923                                  surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,   &
     924                                  surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,   &
     925                                  surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,   &
     926                                  surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,   &
     927                                  surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,   &
     928                                  surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
     929
     930!
     931!--             Prognostic equation for cloud water content
     932                DO  k = nzb+1, nzt
     933                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d *                         &
     934                                                      ( tsc(2) * tend(k,j,i) + &
     935                                                        tsc(3) * tqc_m(k,j,i) )&
     936                                                      - tsc(5) * rdf_sc(k)     &
     937                                                               * qc(k,j,i)     &
     938                                             )                                 &
     939                                       * MERGE( 1.0_wp, 0.0_wp,                &
     940                                                BTEST( wall_flags_0(k,j,i), 0 )&
     941                                              )
     942                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
     943                ENDDO
     944!
     945!--             Calculate tendencies for the next Runge-Kutta step
     946                IF ( timestep_scheme(1:5) == 'runge' )  THEN
     947                   IF ( intermediate_timestep_count == 1 )  THEN
     948                      DO  k = nzb+1, nzt
     949                         tqc_m(k,j,i) = tend(k,j,i)
     950                      ENDDO
     951                   ELSEIF ( intermediate_timestep_count < &
     952                            intermediate_timestep_count_max )  THEN
     953                      DO  k = nzb+1, nzt
     954                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
     955                                           5.3125_wp * tqc_m(k,j,i)
     956                      ENDDO
     957                   ENDIF
     958                ENDIF
     959
     960!
     961!--             Calculate prognostic equation for cloud drop concentration.
     962                tend(:,j,i) = 0.0_wp
     963                IF ( timestep_scheme(1:5) == 'runge' )  THEN
     964                   IF ( ws_scheme_sca )  THEN
     965                      CALL advec_s_ws( i, j, nc, 'nc', flux_s_nc,    &
     966                                    diss_s_nc, flux_l_nc, diss_l_nc, &
     967                                    i_omp_start, tn )
     968                   ELSE
     969                      CALL advec_s_pw( i, j, nc )
     970                   ENDIF
     971                ELSE
     972                   CALL advec_s_up( i, j, nc )
     973                ENDIF
     974                CALL diffusion_s( i, j, nc,                                    &
     975                                  surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,    &
     976                                  surf_def_h(2)%ncsws,                         &
     977                                  surf_lsm_h%ncsws,    surf_usm_h%ncsws,       &
     978                                  surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,    &
     979                                  surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,    &
     980                                  surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,    &
     981                                  surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,    &
     982                                  surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,    &
     983                                  surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
     984
     985!
     986!--             Prognostic equation for cloud drop concentration
     987                DO  k = nzb+1, nzt
     988                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d *                         &
     989                                                      ( tsc(2) * tend(k,j,i) + &
     990                                                        tsc(3) * tnc_m(k,j,i) )&
     991                                                      - tsc(5) * rdf_sc(k)     &
     992                                                               * nc(k,j,i)     &
     993                                             )                                 &
     994                                       * MERGE( 1.0_wp, 0.0_wp,                &
     995                                                BTEST( wall_flags_0(k,j,i), 0 )&
     996                                              )
     997                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
     998                ENDDO
     999!
     1000!--             Calculate tendencies for the next Runge-Kutta step
     1001                IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1002                   IF ( intermediate_timestep_count == 1 )  THEN
     1003                      DO  k = nzb+1, nzt
     1004                         tnc_m(k,j,i) = tend(k,j,i)
     1005                      ENDDO
     1006                   ELSEIF ( intermediate_timestep_count < &
     1007                            intermediate_timestep_count_max )  THEN
     1008                      DO  k = nzb+1, nzt
     1009                         tnc_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +           &
     1010                                           5.3125_wp * tnc_m(k,j,i)
     1011                      ENDDO
     1012                   ENDIF
     1013                ENDIF
     1014
     1015             ENDIF
     1016!
    8941017!--          If required, calculate prognostic equations for rain water content
    8951018!--          and rain drop concentration
     
    18081931
    18091932!
     1933!--    If required, calculate prognostic equations for cloud water content
     1934!--    and cloud drop concentration
     1935       IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     1936
     1937          CALL cpu_log( log_point(67), 'qc-equation', 'start' )
     1938
     1939!
     1940!--       Calculate prognostic equation for cloud water content
     1941          sbt = tsc(2)
     1942          IF ( scalar_advec == 'bc-scheme' )  THEN
     1943
     1944             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     1945!
     1946!--             Bott-Chlond scheme always uses Euler time step. Thus:
     1947                sbt = 1.0_wp
     1948             ENDIF
     1949             tend = 0.0_wp
     1950             CALL advec_s_bc( qc, 'qc' )
     1951
     1952          ENDIF
     1953
     1954!
     1955!--       qc-tendency terms with no communication
     1956          IF ( scalar_advec /= 'bc-scheme' )  THEN
     1957             tend = 0.0_wp
     1958             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1959                IF ( ws_scheme_sca )  THEN
     1960                   CALL advec_s_ws( qc, 'qc' )
     1961                ELSE
     1962                   CALL advec_s_pw( qc )
     1963                ENDIF
     1964             ELSE
     1965                CALL advec_s_up( qc )
     1966             ENDIF
     1967          ENDIF
     1968
     1969          CALL diffusion_s( qc,                                                &
     1970                            surf_def_h(0)%qcsws, surf_def_h(1)%qcsws,          &
     1971                            surf_def_h(2)%qcsws,                               &
     1972                            surf_lsm_h%qcsws,    surf_usm_h%qcsws,             &
     1973                            surf_def_v(0)%qcsws, surf_def_v(1)%qcsws,          &
     1974                            surf_def_v(2)%qcsws, surf_def_v(3)%qcsws,          &
     1975                            surf_lsm_v(0)%qcsws, surf_lsm_v(1)%qcsws,          &
     1976                            surf_lsm_v(2)%qcsws, surf_lsm_v(3)%qcsws,          &
     1977                            surf_usm_v(0)%qcsws, surf_usm_v(1)%qcsws,          &
     1978                            surf_usm_v(2)%qcsws, surf_usm_v(3)%qcsws )
     1979
     1980!
     1981!--       Prognostic equation for cloud water content
     1982          DO  i = nxl, nxr
     1983             DO  j = nys, nyn
     1984                DO  k = nzb+1, nzt
     1985                   qc_p(k,j,i) = qc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
     1986                                                      tsc(3) * tqc_m(k,j,i) )  &
     1987                                                    - tsc(5) * rdf_sc(k) *     &
     1988                                                               qc(k,j,i)       &
     1989                                             )                                 &
     1990                                    * MERGE( 1.0_wp, 0.0_wp,                   &
     1991                                             BTEST( wall_flags_0(k,j,i), 0 )   &
     1992                                          )
     1993                   IF ( qc_p(k,j,i) < 0.0_wp )  qc_p(k,j,i) = 0.0_wp
     1994                ENDDO
     1995             ENDDO
     1996          ENDDO
     1997
     1998!
     1999!--       Calculate tendencies for the next Runge-Kutta step
     2000          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2001             IF ( intermediate_timestep_count == 1 )  THEN
     2002                DO  i = nxl, nxr
     2003                   DO  j = nys, nyn
     2004                      DO  k = nzb+1, nzt
     2005                         tqc_m(k,j,i) = tend(k,j,i)
     2006                      ENDDO
     2007                   ENDDO
     2008                ENDDO
     2009             ELSEIF ( intermediate_timestep_count < &
     2010                      intermediate_timestep_count_max )  THEN
     2011                DO  i = nxl, nxr
     2012                   DO  j = nys, nyn
     2013                      DO  k = nzb+1, nzt
     2014                         tqc_m(k,j,i) =   -9.5625_wp * tend(k,j,i)             &
     2015                                         + 5.3125_wp * tqc_m(k,j,i)
     2016                      ENDDO
     2017                   ENDDO
     2018                ENDDO
     2019             ENDIF
     2020          ENDIF
     2021
     2022          CALL cpu_log( log_point(67), 'qc-equation', 'stop' )
     2023          CALL cpu_log( log_point(68), 'nc-equation', 'start' )
     2024
     2025!
     2026!--       Calculate prognostic equation for cloud drop concentration
     2027          sbt = tsc(2)
     2028          IF ( scalar_advec == 'bc-scheme' )  THEN
     2029
     2030             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
     2031!
     2032!--             Bott-Chlond scheme always uses Euler time step. Thus:
     2033                sbt = 1.0_wp
     2034             ENDIF
     2035             tend = 0.0_wp
     2036             CALL advec_s_bc( nc, 'nc' )
     2037
     2038          ENDIF
     2039
     2040!
     2041!--       nc-tendency terms with no communication
     2042          IF ( scalar_advec /= 'bc-scheme' )  THEN
     2043             tend = 0.0_wp
     2044             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2045                IF ( ws_scheme_sca )  THEN
     2046                   CALL advec_s_ws( nc, 'nc' )
     2047                ELSE
     2048                   CALL advec_s_pw( nc )
     2049                ENDIF
     2050             ELSE
     2051                CALL advec_s_up( nc )
     2052             ENDIF
     2053          ENDIF
     2054
     2055          CALL diffusion_s( nc,                                                &
     2056                            surf_def_h(0)%ncsws, surf_def_h(1)%ncsws,          &
     2057                            surf_def_h(2)%ncsws,                               &
     2058                            surf_lsm_h%ncsws,    surf_usm_h%ncsws,             &
     2059                            surf_def_v(0)%ncsws, surf_def_v(1)%ncsws,          &
     2060                            surf_def_v(2)%ncsws, surf_def_v(3)%ncsws,          &
     2061                            surf_lsm_v(0)%ncsws, surf_lsm_v(1)%ncsws,          &
     2062                            surf_lsm_v(2)%ncsws, surf_lsm_v(3)%ncsws,          &
     2063                            surf_usm_v(0)%ncsws, surf_usm_v(1)%ncsws,          &
     2064                            surf_usm_v(2)%ncsws, surf_usm_v(3)%ncsws )
     2065
     2066!
     2067!--       Prognostic equation for cloud drop concentration
     2068          DO  i = nxl, nxr
     2069             DO  j = nys, nyn
     2070                DO  k = nzb+1, nzt
     2071                   nc_p(k,j,i) = nc(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +   &
     2072                                                      tsc(3) * tnc_m(k,j,i) )  &
     2073                                                    - tsc(5) * rdf_sc(k) *     &
     2074                                                               nc(k,j,i)       &
     2075                                             )                                 &
     2076                                   * MERGE( 1.0_wp, 0.0_wp,                    &
     2077                                             BTEST( wall_flags_0(k,j,i), 0 )   &
     2078                                          )
     2079                   IF ( nc_p(k,j,i) < 0.0_wp )  nc_p(k,j,i) = 0.0_wp
     2080                ENDDO
     2081             ENDDO
     2082          ENDDO
     2083
     2084!
     2085!--       Calculate tendencies for the next Runge-Kutta step
     2086          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2087             IF ( intermediate_timestep_count == 1 )  THEN
     2088                DO  i = nxl, nxr
     2089                   DO  j = nys, nyn
     2090                      DO  k = nzb+1, nzt
     2091                         tnc_m(k,j,i) = tend(k,j,i)
     2092                      ENDDO
     2093                   ENDDO
     2094                ENDDO
     2095             ELSEIF ( intermediate_timestep_count < &
     2096                      intermediate_timestep_count_max )  THEN
     2097                DO  i = nxl, nxr
     2098                   DO  j = nys, nyn
     2099                      DO  k = nzb+1, nzt
     2100                         tnc_m(k,j,i) =  -9.5625_wp * tend(k,j,i)             &
     2101                                         + 5.3125_wp * tnc_m(k,j,i)
     2102                      ENDDO
     2103                   ENDDO
     2104                ENDDO
     2105             ENDIF
     2106          ENDIF
     2107
     2108          CALL cpu_log( log_point(68), 'nc-equation', 'stop' )
     2109
     2110       ENDIF
     2111!
    18102112!--    If required, calculate prognostic equations for rain water content
    18112113!--    and rain drop concentration
  • palm/trunk/SOURCE/read_3d_binary.f90

    r2269 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2269 2017-06-09 11:57:32Z suehring
    2732! Enable restart runs for urban_surface_mod
    2833!
     
    127132
    128133    USE arrays_3d,                                                             &
    129         ONLY:  e, kh, km, p, pt, q, ql, qc, nr, prr, precipitation_amount, qr, &
    130                s, sa, u, u_m_l, u_m_n, u_m_r, u_m_s, v, v_m_l, v_m_n, v_m_r,   &
    131                v_m_s, vpt, w, w_m_l, w_m_n, w_m_r, w_m_s
     134        ONLY:  e, kh, km, p, pt, q, ql, qc, nc, nr, prr, precipitation_amount, &
     135               qr, s, sa, u, u_m_l, u_m_n, u_m_r, u_m_s, v, v_m_l, v_m_n,      &
     136               v_m_r, v_m_s, vpt, w, w_m_l, w_m_n, w_m_r, w_m_s
    132137
    133138    USE averaging
     
    496501                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    497502
     503                CASE ( 'nc' )
     504                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     505                   nc(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     506                                   tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     507
     508                CASE ( 'nc_av' )
     509                   IF ( .NOT. ALLOCATED( nc_av ) )  THEN
     510                      ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     511                   ENDIF
     512                   IF ( k == 1 )  READ ( 13 )  tmp_3d
     513                   nc_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     514                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     515
     516
    498517                CASE ( 'nr' )
    499518                   IF ( k == 1 )  READ ( 13 )  tmp_3d
  • palm/trunk/SOURCE/sum_up_3d_data.f90

    r2233 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2233 2017-05-30 18:08:54Z suehring
    2732!
    2833! 2232 2017-05-30 17:47:52Z suehring
     
    135140
    136141    USE arrays_3d,                                                             &
    137         ONLY:  dzw, e, nr, p, pt, precipitation_rate, q, qc, ql, ql_c, ql_v,   &
    138                qr, rho_ocean, s, sa, u, v, vpt, w
     142        ONLY:  dzw, e, nc, nr, p, pt, precipitation_rate, q, qc, ql, ql_c,     &
     143               ql_v, qr, rho_ocean, s, sa, u, v, vpt, w
    139144
    140145    USE averaging,                                                             &
    141         ONLY:  e_av, lpt_av, lwp_av, nr_av, ol_av, p_av, pc_av, pr_av, prr_av, &
    142                precipitation_rate_av, pt_av, q_av, qc_av, ql_av, ql_c_av,      &
    143                ql_v_av, ql_vp_av, qr_av, qsws_av, qv_av, rho_ocean_av, s_av, sa_av,  &
    144                shf_av, ssws_av, ts_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, &
    145                z0h_av, z0q_av
     146        ONLY:  e_av, lpt_av, lwp_av, nc_av, nr_av, ol_av, p_av, pc_av, pr_av, &
     147               prr_av, precipitation_rate_av, pt_av, q_av, qc_av, ql_av,       &
     148               ql_c_av, ql_v_av, ql_vp_av, qr_av, qsws_av, qv_av, rho_ocean_av,&
     149               s_av, sa_av, shf_av, ssws_av, ts_av, u_av, us_av, v_av, vpt_av, &
     150               w_av, z0_av, z0h_av, z0q_av
    146151
    147152    USE cloud_parameters,                                                      &
     
    230235                lwp_av = 0.0_wp
    231236
     237             CASE ( 'nc' )
     238                IF ( .NOT. ALLOCATED( nc_av ) )  THEN
     239                   ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     240                ENDIF
     241                nc_av = 0.0_wp
     242
    232243             CASE ( 'nr' )
    233244                IF ( .NOT. ALLOCATED( nr_av ) )  THEN
     
    483494                   lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i)            &
    484495                                               * dzw(1:nzt+1) ) * rho_surface
     496                ENDDO
     497             ENDDO
     498
     499          CASE ( 'nc' )
     500             DO  i = nxlg, nxrg
     501                DO  j = nysg, nyng
     502                   DO  k = nzb, nzt+1
     503                      nc_av(k,j,i) = nc_av(k,j,i) + nc(k,j,i)
     504                   ENDDO
    485505                ENDDO
    486506             ENDDO
  • palm/trunk/SOURCE/surface_layer_fluxes_mod.f90

    r2281 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2281 2017-06-13 11:34:50Z suehring
    2732! Clean-up unnecessary index access to surface type
    2833!
     
    178183
    179184    USE arrays_3d,                                                             &
    180         ONLY:  e, kh, nr, pt, q, ql, qr, s, u, v, vpt, w, zu, zw, drho_air_zw, &
    181                rho_air_zw
     185        ONLY:  e, kh, nc, nr, pt, q, ql, qc, qr, s, u, v, vpt, w, zu, zw,      &
     186               drho_air_zw, rho_air_zw
    182187
    183188    USE cloud_parameters,                                                      &
     
    195200               intermediate_timestep_count, intermediate_timestep_count_max,   &
    196201               land_surface, large_scale_forcing, lsf_surf,                    &
    197                message_string, microphysics_seifert, most_method, neutral,     &
    198                passive_scalar, pt_surface, q_surface, run_coupled,             &
    199                surface_pressure, simulated_time, terminate_run,                &
     202               message_string, microphysics_morrison, microphysics_seifert,    &
     203               most_method, neutral, passive_scalar, pt_surface, q_surface,    &
     204               run_coupled, surface_pressure, simulated_time, terminate_run,   &
    200205               urban_surface, zeta_max, zeta_min
    201206
     
    15261531
    15271532!
    1528 !-- Calculate the other MOST scaling parameters theta*, q*, (qr*, nr*)
     1533!-- Calculate the other MOST scaling parameters theta*, q*, (qc*, qr*, nc*, nr*)
    15291534    SUBROUTINE calc_scaling_parameters
    15301535
     
    17491754       ENDIF
    17501755
     1756!
     1757!--    If required compute qc* and nc*
     1758       IF ( cloud_physics  .AND.  microphysics_morrison  .AND.                 &
     1759            .NOT. surf_vertical )  THEN
     1760          !$OMP PARALLEL DO PRIVATE( i, j, k, z_mo )
     1761          DO  m = 1, surf%ns   
     1762
     1763             i    = surf%i(m)           
     1764             j    = surf%j(m)
     1765             k    = surf%k(m)
     1766
     1767             z_mo = surf%z_mo(m)
     1768
     1769             surf%qcs(m) = kappa * ( qc(k,j,i) - qc(k-1,j,i) )                 &
     1770                                 / ( LOG( z_mo / surf%z0q(m) )                 &
     1771                                 - psi_h( z_mo / surf%ol(m) )                  &
     1772                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
     1773
     1774             surf%ncs(m) = kappa * ( nc(k,j,i) - nc(k-1,j,i) )                 &
     1775                                 / ( LOG( z_mo / surf%z0q(m) )                 &
     1776                                 - psi_h( z_mo / surf%ol(m) )                  &
     1777                                 + psi_h( surf%z0q(m) / surf%ol(m) ) )
     1778          ENDDO
     1779
     1780       ENDIF
    17511781
    17521782!
     
    17811811
    17821812!
    1783 !-- Calculate surface fluxes usws, vsws, shf, qsws, (qrsws, nrsws)
     1813!-- Calculate surface fluxes usws, vsws, shf, qsws, (qcsws, qrsws, ncsws, nrsws)
    17841814    SUBROUTINE calc_surface_fluxes
    17851815
     
    19271957
    19281958             ENDDO
    1929           ENDIF       
     1959          ENDIF   
     1960!
     1961!--       Compute (turbulent) fluxes of cloud water content and cloud drop conc.
     1962          IF ( cloud_physics  .AND.  microphysics_morrison  .AND.              &
     1963               .NOT. downward)  THEN
     1964             !$OMP PARALLEL DO PRIVATE( i, j )
     1965             DO  m = 1, surf%ns   
     1966
     1967                i    = surf%i(m)           
     1968                j    = surf%j(m)
     1969
     1970                surf%qcsws(m) = -surf%qcs(m) * surf%us(m)
     1971                surf%ncsws(m) = -surf%ncs(m) * surf%us(m)
     1972             ENDDO
     1973          ENDIF   
    19301974!
    19311975!--       Compute (turbulent) fluxes of rain water content and rain drop conc.
    1932           IF ( cloud_physics  .AND.  microphysics_seifert  .AND.              &
     1976          IF ( cloud_physics  .AND.  microphysics_seifert  .AND.               &
    19331977               .NOT. downward)  THEN
    19341978             !$OMP PARALLEL DO PRIVATE( i, j )
  • palm/trunk/SOURCE/surface_mod.f90

    r2270 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2270 2017-06-09 12:18:47Z maronga
    2732! Parameters removed/added due to changes in the LSM
    2833!
     
    104109       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qs        !< scaling parameter humidity
    105110       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ss        !< scaling parameter passive scalar
     111       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qcs       !< scaling parameter qc
     112       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ncs       !< scaling parameter nc
    106113       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qrs       !< scaling parameter qr
    107114       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrs       !< scaling parameter nr
     
    124131       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qsws      !< surface flux latent heat
    125132       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ssws      !< surface flux passive scalar
     133       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qcsws     !< surface flux qc
     134       REAL(wp), DIMENSION(:), ALLOCATABLE ::  ncsws     !< surface flux nc
    126135       REAL(wp), DIMENSION(:), ALLOCATABLE ::  qrsws     !< surface flux qr
    127136       REAL(wp), DIMENSION(:), ALLOCATABLE ::  nrsws     !< surface flux nr
     
    719728!
    720729!--       
     730       IF ( cloud_physics .AND. microphysics_morrison)  THEN
     731          ALLOCATE ( surfaces%qcs(1:surfaces%ns)   )
     732          ALLOCATE ( surfaces%ncs(1:surfaces%ns)   )
     733          ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
     734          ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
     735       ENDIF
     736!
     737!--       
    721738       IF ( cloud_physics .AND. microphysics_seifert)  THEN
    722739          ALLOCATE ( surfaces%qrs(1:surfaces%ns)   )
     
    775792          ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
    776793       ENDIF
     794!
     795!--       
     796       IF ( cloud_physics .AND. microphysics_morrison)  THEN
     797          ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
     798          ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
     799       ENDIF
    777800!
    778801!--       
     
    877900          ALLOCATE ( surfaces%ssws(1:surfaces%ns) )
    878901       ENDIF
     902
     903       IF ( cloud_physics .AND. microphysics_seifert)  THEN
     904          ALLOCATE ( surfaces%qcs(1:surfaces%ns)   )
     905          ALLOCATE ( surfaces%ncs(1:surfaces%ns)   )
     906          ALLOCATE ( surfaces%qcsws(1:surfaces%ns) )
     907          ALLOCATE ( surfaces%ncsws(1:surfaces%ns) )
     908       ENDIF
    879909
    880910       IF ( cloud_physics .AND. microphysics_seifert)  THEN
     
    13081338             IF ( humidity )  THEN
    13091339                surf%qs(num_h)   = 0.0_wp
     1340                IF ( cloud_physics .AND. microphysics_morrison)  THEN
     1341                   surf%qcs(num_h) = 0.0_wp
     1342                   surf%ncs(num_h) = 0.0_wp
     1343   
     1344                   surf%qcsws(num_h) = 0.0_wp
     1345                   surf%ncsws(num_h) = 0.0_wp
     1346
     1347                ENDIF
    13101348                IF ( cloud_physics .AND. microphysics_seifert)  THEN
    13111349                   surf%qrs(num_h) = 0.0_wp
     
    14471485             IF ( humidity )  THEN
    14481486             surf%qsws = 0.0_wp
     1487                IF ( cloud_physics  .AND.  microphysics_morrison ) THEN
     1488                   surf%ncsws = 0.0_wp
     1489                   surf%qcsws = 0.0_wp
     1490                ENDIF
    14491491                IF ( cloud_physics  .AND.  microphysics_seifert ) THEN
    14501492                   surf%nrsws = 0.0_wp
     
    15681610!
    15691611!--             Following wall fluxes are assumed to be zero
     1612                IF ( cloud_physics .AND. microphysics_morrison)  THEN
     1613                   surf%qcs(num_v) = 0.0_wp
     1614                   surf%ncs(num_v) = 0.0_wp
     1615   
     1616                   surf%qcsws(num_v) = 0.0_wp
     1617                   surf%ncsws(num_v) = 0.0_wp
     1618                ENDIF
    15701619                IF ( cloud_physics .AND. microphysics_seifert)  THEN
    15711620                   surf%qrs(num_v) = 0.0_wp
     
    16561705                   IF ( ALLOCATED( surf_def_h(l)%ss ) )                        &
    16571706                      surf_h(l)%ss(mm(l))      = surf_def_h(l)%ss(m)
     1707                   IF ( ALLOCATED( surf_def_h(l)%qcs ) )                       &
     1708                      surf_h(l)%qcs(mm(l))     = surf_def_h(l)%qcs(m)
     1709                   IF ( ALLOCATED( surf_def_h(l)%ncs ) )                       &
     1710                      surf_h(l)%ncs(mm(l))     = surf_def_h(l)%ncs(m)
    16581711                   IF ( ALLOCATED( surf_def_h(l)%qrs ) )                       &
    16591712                      surf_h(l)%qrs(mm(l))     = surf_def_h(l)%qrs(m)
     
    16741727                   IF ( ALLOCATED( surf_def_h(l)%ssws ) )                      &
    16751728                      surf_h(l)%qsws(mm(l))    = surf_def_h(l)%ssws(m)
     1729                   IF ( ALLOCATED( surf_def_h(l)%ncsws ) )                     &
     1730                      surf_h(l)%ncsws(mm(l))   = surf_def_h(l)%ncsws(m)
    16761731                   IF ( ALLOCATED( surf_def_h(l)%nrsws ) )                     &
    16771732                      surf_h(l)%nrsws(mm(l))   = surf_def_h(l)%nrsws(m)
     
    16931748                      IF ( ALLOCATED( surf_lsm_h%ss ) )                        &
    16941749                         surf_h(0)%ss(mm(0))      = surf_lsm_h%ss(m)
     1750                      IF ( ALLOCATED( surf_lsm_h%qcs ) )                       &
     1751                         surf_h(0)%qcs(mm(0))     = surf_lsm_h%qcs(m)
     1752                      IF ( ALLOCATED( surf_lsm_h%ncs ) )                       &
     1753                         surf_h(0)%ncs(mm(0))     = surf_lsm_h%ncs(m)
    16951754                      IF ( ALLOCATED( surf_lsm_h%qrs ) )                       &
    16961755                         surf_h(0)%qrs(mm(0))     = surf_lsm_h%qrs(m)
     
    17111770                      IF ( ALLOCATED( surf_lsm_h%ssws ) )                      &
    17121771                         surf_h(0)%qsws(mm(0))    = surf_lsm_h%ssws(m)
     1772                      IF ( ALLOCATED( surf_lsm_h%ncsws ) )                     &
     1773                         surf_h(0)%ncsws(mm(0))   = surf_lsm_h%ncsws(m)
    17131774                      IF ( ALLOCATED( surf_lsm_h%nrsws ) )                     &
    17141775                         surf_h(0)%nrsws(mm(0))   = surf_lsm_h%nrsws(m)
     
    17301791                      IF ( ALLOCATED( surf_usm_h%ss ) )                        &
    17311792                         surf_h(0)%ss(mm(0))      = surf_usm_h%ss(m)
     1793                      IF ( ALLOCATED( surf_usm_h%qcs ) )                       &
     1794                         surf_h(0)%qcs(mm(0))     = surf_usm_h%qcs(m)
     1795                      IF ( ALLOCATED( surf_usm_h%ncs ) )                       &
     1796                         surf_h(0)%ncs(mm(0))     = surf_usm_h%ncs(m)
    17321797                      IF ( ALLOCATED( surf_usm_h%qrs ) )                       &
    17331798                         surf_h(0)%qrs(mm(0))     = surf_usm_h%qrs(m)
     
    17481813                      IF ( ALLOCATED( surf_usm_h%ssws ) )                      &
    17491814                         surf_h(0)%qsws(mm(0))    = surf_usm_h%ssws(m)
     1815                      IF ( ALLOCATED( surf_usm_h%ncsws ) )                     &
     1816                         surf_h(0)%ncsws(mm(0))   = surf_usm_h%ncsws(m)
    17501817                      IF ( ALLOCATED( surf_usm_h%nrsws ) )                     &
    17511818                         surf_h(0)%nrsws(mm(0))   = surf_usm_h%nrsws(m)
     
    17911858                   IF ( ALLOCATED( surf_def_v(l)%ss ) )                        &
    17921859                      surf_v(l)%ss(mm(l))      = surf_def_v(l)%ss(m)
     1860                   IF ( ALLOCATED( surf_def_v(l)%qcs ) )                       &
     1861                      surf_v(l)%qcs(mm(l))     = surf_def_v(l)%qcs(m)
     1862                   IF ( ALLOCATED( surf_def_v(l)%ncs ) )                       &
     1863                      surf_v(l)%ncs(mm(l))     = surf_def_v(l)%ncs(m)
    17931864                   IF ( ALLOCATED( surf_def_v(l)%qrs ) )                       &
    17941865                      surf_v(l)%qrs(mm(l))     = surf_def_v(l)%qrs(m)
     
    18051876                   IF ( ALLOCATED( surf_def_v(l)%ssws ) )                      &
    18061877                      surf_v(l)%qsws(mm(l))    = surf_def_v(l)%ssws(m)
     1878                   IF ( ALLOCATED( surf_def_v(l)%ncsws ) )                     &
     1879                      surf_v(l)%ncsws(mm(l))   = surf_def_v(l)%ncsws(m)
    18071880                   IF ( ALLOCATED( surf_def_v(l)%nrsws ) )                     &
    18081881                      surf_v(l)%nrsws(mm(l))   = surf_def_v(l)%nrsws(m)
     
    18291902                   IF ( ALLOCATED( surf_lsm_v(l)%ss ) )                        &
    18301903                      surf_v(l)%ss(mm(l))      = surf_lsm_v(l)%ss(m)
     1904                   IF ( ALLOCATED( surf_lsm_v(l)%qcs ) )                       &
     1905                      surf_v(l)%qcs(mm(l))     = surf_lsm_v(l)%qcs(m)
     1906                   IF ( ALLOCATED( surf_lsm_v(l)%ncs ) )                       &
     1907                      surf_v(l)%ncs(mm(l))     = surf_lsm_v(l)%ncs(m)
    18311908                   IF ( ALLOCATED( surf_lsm_v(l)%qrs ) )                       &
    18321909                      surf_v(l)%qrs(mm(l))     = surf_lsm_v(l)%qrs(m)
     
    18471924                   IF ( ALLOCATED( surf_lsm_v(l)%ssws ) )                      &
    18481925                      surf_v(l)%qsws(mm(l))    = surf_lsm_v(l)%ssws(m)
     1926                   IF ( ALLOCATED( surf_lsm_v(l)%ncsws ) )                     &
     1927                      surf_v(l)%ncsws(mm(l))   = surf_lsm_v(l)%ncsws(m)
    18491928                   IF ( ALLOCATED( surf_lsm_v(l)%nrsws ) )                     &
    18501929                      surf_v(l)%nrsws(mm(l))   = surf_lsm_v(l)%nrsws(m)
     
    18711950                   IF ( ALLOCATED( surf_usm_v(l)%ss ) )                        &
    18721951                      surf_v(l)%ss(mm(l))      = surf_usm_v(l)%ss(m)
     1952                   IF ( ALLOCATED( surf_usm_v(l)%qcs ) )                       &
     1953                      surf_v(l)%qcs(mm(l))     = surf_usm_v(l)%qcs(m)
     1954                   IF ( ALLOCATED( surf_usm_v(l)%ncs ) )                       &
     1955                      surf_v(l)%ncs(mm(l))     = surf_usm_v(l)%ncs(m)
    18731956                   IF ( ALLOCATED( surf_usm_v(l)%qrs ) )                       &
    18741957                      surf_v(l)%qrs(mm(l))     = surf_usm_v(l)%qrs(m)
     
    18891972                   IF ( ALLOCATED( surf_usm_v(l)%ssws ) )                      &
    18901973                      surf_v(l)%qsws(mm(l))    = surf_usm_v(l)%ssws(m)
     1974                   IF ( ALLOCATED( surf_usm_v(l)%ncsws ) )                     &
     1975                      surf_v(l)%ncsws(mm(l))   = surf_usm_v(l)%ncsws(m)
    18911976                   IF ( ALLOCATED( surf_usm_v(l)%nrsws ) )                     &
    18921977                      surf_v(l)%nrsws(mm(l))   = surf_usm_v(l)%nrsws(m)
     
    19462031             WRITE ( 14 )  surf_h(l)%ss
    19472032          ENDIF
     2033          WRITE ( 14 )  'surf_h(' // dum // ')%qcs                 '
     2034          IF ( ALLOCATED ( surf_h(l)%qcs ) )  THEN 
     2035             WRITE ( 14 )  surf_h(l)%qcs
     2036          ENDIF
     2037          WRITE ( 14 )  'surf_h(' // dum // ')%ncs                 ' 
     2038          IF ( ALLOCATED ( surf_h(l)%ncs ) )  THEN
     2039             WRITE ( 14 )  surf_h(l)%ncs
     2040          ENDIF
    19482041          WRITE ( 14 )  'surf_h(' // dum // ')%qrs                 '
    19492042          IF ( ALLOCATED ( surf_h(l)%qrs ) )  THEN 
     
    19812074          IF ( ALLOCATED ( surf_h(l)%ssws ) )  THEN
    19822075             WRITE ( 14 )  surf_h(l)%ssws
     2076          ENDIF
     2077          WRITE ( 14 )  'surf_h(' // dum // ')%qcsws               ' 
     2078          IF ( ALLOCATED ( surf_h(l)%qcsws ) )  THEN
     2079             WRITE ( 14 )  surf_h(l)%qcsws
     2080          ENDIF
     2081          WRITE ( 14 )  'surf_h(' // dum // ')%ncsws               ' 
     2082          IF ( ALLOCATED ( surf_h(l)%ncsws ) )  THEN
     2083             WRITE ( 14 )  surf_h(l)%ncsws
    19832084          ENDIF
    19842085          WRITE ( 14 )  'surf_h(' // dum // ')%qrsws               ' 
     
    20212122             WRITE ( 14 )  surf_v(l)%ss
    20222123          ENDIF
     2124          WRITE ( 14 )  'surf_v(' // dum // ')%qcs                 ' 
     2125          IF ( ALLOCATED ( surf_v(l)%qcs ) )  THEN
     2126             WRITE ( 14 )  surf_v(l)%qcs
     2127          ENDIF
     2128          WRITE ( 14 )  'surf_v(' // dum // ')%ncs                 '
     2129          IF ( ALLOCATED ( surf_v(l)%ncs ) )  THEN
     2130             WRITE ( 14 )  surf_v(l)%ncs
     2131          ENDIF
    20232132          WRITE ( 14 )  'surf_v(' // dum // ')%qrs                 ' 
    20242133          IF ( ALLOCATED ( surf_v(l)%qrs ) )  THEN
     
    20482157          IF ( ALLOCATED ( surf_v(l)%ssws ) )  THEN
    20492158             WRITE ( 14 )  surf_v(l)%ssws
     2159          ENDIF
     2160          WRITE ( 14 )  'surf_v(' // dum // ')%qcsws               ' 
     2161          IF ( ALLOCATED ( surf_v(l)%qcsws ) )  THEN
     2162             WRITE ( 14 )  surf_v(l)%qcsws
     2163          ENDIF
     2164          WRITE ( 14 )  'surf_v(' // dum // ')%ncsws               ' 
     2165          IF ( ALLOCATED ( surf_v(l)%ncsws ) )  THEN
     2166             WRITE ( 14 )  surf_v(l)%ncsws
    20502167          ENDIF
    20512168          WRITE ( 14 )  'surf_v(' // dum // ')%qrsws               ' 
     
    22402357                      IF ( ALLOCATED( surf_h(0)%ss )  .AND.  kk == 1 )         &
    22412358                         READ ( 13 )  surf_h(0)%ss
     2359                   CASE ( 'surf_h(0)%qcs' )         
     2360                      IF ( ALLOCATED( surf_h(0)%qcs )  .AND.  kk == 1 )        &
     2361                         READ ( 13 )  surf_h(0)%qcs
     2362                   CASE ( 'surf_h(0)%ncs' )         
     2363                      IF ( ALLOCATED( surf_h(0)%ncs )  .AND.  kk == 1 )        &
     2364                         READ ( 13 )  surf_h(0)%ncs
    22422365                   CASE ( 'surf_h(0)%qrs' )         
    22432366                      IF ( ALLOCATED( surf_h(0)%qrs )  .AND.  kk == 1 )        &
     
    22672390                      IF ( ALLOCATED( surf_h(0)%ssws )  .AND.  kk == 1 )       &
    22682391                         READ ( 13 )  surf_h(0)%ssws
     2392                   CASE ( 'surf_h(0)%qcsws' )         
     2393                      IF ( ALLOCATED( surf_h(0)%qcsws )  .AND.  kk == 1 )      &
     2394                         READ ( 13 )  surf_h(0)%qcsws
     2395                   CASE ( 'surf_h(0)%ncsws' )         
     2396                      IF ( ALLOCATED( surf_h(0)%ncsws )  .AND.  kk == 1 )      &
     2397                         READ ( 13 )  surf_h(0)%ncsws
    22692398                   CASE ( 'surf_h(0)%qrsws' )         
    22702399                      IF ( ALLOCATED( surf_h(0)%qrsws )  .AND.  kk == 1 )      &
     
    22962425                      IF ( ALLOCATED( surf_h(1)%ss )  .AND.  kk == 1 )         &
    22972426                         READ ( 13 )  surf_h(1)%ss
     2427                   CASE ( 'surf_h(1)%qcs' )         
     2428                      IF ( ALLOCATED( surf_h(1)%qcs )  .AND.  kk == 1 )        &
     2429                         READ ( 13 )  surf_h(1)%qcs
     2430                   CASE ( 'surf_h(1)%ncs' )         
     2431                      IF ( ALLOCATED( surf_h(1)%ncs )  .AND.  kk == 1 )        &
     2432                         READ ( 13 )  surf_h(1)%ncs
    22982433                   CASE ( 'surf_h(1)%qrs' )         
    22992434                      IF ( ALLOCATED( surf_h(1)%qrs )  .AND.  kk == 1 )        &
     
    23232458                      IF ( ALLOCATED( surf_h(1)%ssws )  .AND.  kk == 1 )       &
    23242459                         READ ( 13 )  surf_h(1)%ssws
     2460                   CASE ( 'surf_h(1)%qcsws' )         
     2461                      IF ( ALLOCATED( surf_h(1)%qcsws )  .AND.  kk == 1 )      &
     2462                         READ ( 13 )  surf_h(1)%qcsws
     2463                   CASE ( 'surf_h(1)%ncsws' )         
     2464                      IF ( ALLOCATED( surf_h(1)%ncsws )  .AND.  kk == 1 )      &
     2465                         READ ( 13 )  surf_h(1)%ncsws
    23252466                   CASE ( 'surf_h(1)%qrsws' )         
    23262467                      IF ( ALLOCATED( surf_h(1)%qrsws )  .AND.  kk == 1 )      &
     
    23522493                      IF ( ALLOCATED( surf_h(2)%ss )  .AND.  kk == 1 )         &
    23532494                         READ ( 13 )  surf_h(2)%ss
     2495                   CASE ( 'surf_h(2)%qcs' )         
     2496                      IF ( ALLOCATED( surf_h(2)%qcs )  .AND.  kk == 1 )        &
     2497                         READ ( 13 )  surf_h(2)%qcs
     2498                   CASE ( 'surf_h(2)%ncs' )         
     2499                      IF ( ALLOCATED( surf_h(2)%ncs )  .AND.  kk == 1 )        &
     2500                         READ ( 13 )  surf_h(2)%ncs
    23542501                   CASE ( 'surf_h(2)%qrs' )         
    23552502                      IF ( ALLOCATED( surf_h(2)%qrs )  .AND.  kk == 1 )        &
     
    23792526                      IF ( ALLOCATED( surf_h(2)%ssws )  .AND.  kk == 1 )       &
    23802527                         READ ( 13 )  surf_h(2)%ssws
     2528                   CASE ( 'surf_h(2)%qcsws' )         
     2529                      IF ( ALLOCATED( surf_h(2)%qcsws )  .AND.  kk == 1 )      &
     2530                         READ ( 13 )  surf_h(2)%qcsws
     2531                   CASE ( 'surf_h(2)%ncsws' )         
     2532                      IF ( ALLOCATED( surf_h(2)%ncsws )  .AND.  kk == 1 )      &
     2533                         READ ( 13 )  surf_h(2)%ncsws
    23812534                   CASE ( 'surf_h(2)%qrsws' )         
    23822535                      IF ( ALLOCATED( surf_h(2)%qrsws )  .AND.  kk == 1 )      &
     
    24102563                      IF ( ALLOCATED( surf_v(0)%ss )  .AND.  kk == 1 )         &
    24112564                         READ ( 13 )  surf_v(0)%ss
     2565                   CASE ( 'surf_v(0)%qcs' )         
     2566                      IF ( ALLOCATED( surf_v(0)%qcs )  .AND.  kk == 1 )        &
     2567                         READ ( 13 )  surf_v(0)%qcs
     2568                   CASE ( 'surf_v(0)%ncs' )         
     2569                      IF ( ALLOCATED( surf_v(0)%ncs )  .AND.  kk == 1 )        &
     2570                         READ ( 13 )  surf_v(0)%ncs
    24122571                   CASE ( 'surf_v(0)%qrs' )         
    24132572                      IF ( ALLOCATED( surf_v(0)%qrs )  .AND.  kk == 1 )        &
     
    24312590                      IF ( ALLOCATED( surf_v(0)%ssws )  .AND.  kk == 1 )       &
    24322591                         READ ( 13 )  surf_v(0)%ssws
     2592                   CASE ( 'surf_v(0)%qcsws' )         
     2593                      IF ( ALLOCATED( surf_v(0)%qcsws )  .AND.  kk == 1 )      &
     2594                         READ ( 13 )  surf_v(0)%qcsws
     2595                   CASE ( 'surf_v(0)%ncsws' )         
     2596                      IF ( ALLOCATED( surf_v(0)%ncsws )  .AND.  kk == 1 )      &
     2597                         READ ( 13 )  surf_v(0)%ncsws
    24332598                   CASE ( 'surf_v(0)%qrsws' )         
    24342599                      IF ( ALLOCATED( surf_v(0)%qrsws )  .AND.  kk == 1 )      &
     
    24692634                      IF ( ALLOCATED( surf_v(1)%ss )  .AND.  kk == 1 )         &
    24702635                         READ ( 13 )  surf_v(1)%ss
     2636                   CASE ( 'surf_v(1)%qcs' )         
     2637                      IF ( ALLOCATED( surf_v(1)%qcs )  .AND.  kk == 1 )        &
     2638                         READ ( 13 )  surf_v(1)%qcs
     2639                   CASE ( 'surf_v(1)%ncs' )         
     2640                      IF ( ALLOCATED( surf_v(1)%ncs )  .AND.  kk == 1 )        &
     2641                         READ ( 13 )  surf_v(1)%ncs
    24712642                   CASE ( 'surf_v(1)%qrs' )         
    24722643                      IF ( ALLOCATED( surf_v(1)%qrs )  .AND.  kk == 1 )        &
     
    24902661                      IF ( ALLOCATED( surf_v(1)%ssws )  .AND.  kk == 1 )       &
    24912662                         READ ( 13 )  surf_v(1)%ssws
     2663                   CASE ( 'surf_v(1)%qcsws' )         
     2664                      IF ( ALLOCATED( surf_v(1)%qcsws )  .AND.  kk == 1 )      &
     2665                         READ ( 13 )  surf_v(1)%qcsws
     2666                   CASE ( 'surf_v(1)%ncsws' )         
     2667                      IF ( ALLOCATED( surf_v(1)%ncsws )  .AND.  kk == 1 )      &
     2668                         READ ( 13 )  surf_v(1)%ncsws
    24922669                   CASE ( 'surf_v(1)%qrsws' )         
    24932670                      IF ( ALLOCATED( surf_v(1)%qrsws )  .AND.  kk == 1 )      &
     
    25282705                      IF ( ALLOCATED( surf_v(2)%ss )  .AND.  kk == 1 )         &
    25292706                         READ ( 13 )  surf_v(2)%ss
     2707                   CASE ( 'surf_v(2)%qcs' )         
     2708                      IF ( ALLOCATED( surf_v(2)%qcs )  .AND.  kk == 1 )        &
     2709                         READ ( 13 )  surf_v(2)%qcs
     2710                   CASE ( 'surf_v(2)%ncs' )         
     2711                      IF ( ALLOCATED( surf_v(2)%ncs )  .AND.  kk == 1 )        &
     2712                         READ ( 13 )  surf_v(2)%ncs
    25302713                   CASE ( 'surf_v(2)%qrs' )         
    25312714                      IF ( ALLOCATED( surf_v(2)%qrs )  .AND.  kk == 1 )        &
     
    25492732                      IF ( ALLOCATED( surf_v(2)%ssws )  .AND.  kk == 1 )       &
    25502733                         READ ( 13 )  surf_v(2)%ssws
     2734                   CASE ( 'surf_v(2)%qcsws' )         
     2735                      IF ( ALLOCATED( surf_v(2)%qcsws )  .AND.  kk == 1 )      &
     2736                         READ ( 13 )  surf_v(2)%qcsws
     2737                   CASE ( 'surf_v(2)%ncsws' )         
     2738                      IF ( ALLOCATED( surf_v(2)%ncsws )  .AND.  kk == 1 )      &
     2739                         READ ( 13 )  surf_v(2)%ncsws
    25512740                   CASE ( 'surf_v(2)%qrsws' )         
    25522741                      IF ( ALLOCATED( surf_v(2)%qrsws )  .AND.  kk == 1 )      &
     
    25872776                      IF ( ALLOCATED( surf_v(3)%ss )  .AND.  kk == 1 )         &
    25882777                         READ ( 13 )  surf_v(3)%ss
     2778                   CASE ( 'surf_v(3)%qcs' )         
     2779                      IF ( ALLOCATED( surf_v(3)%qcs )  .AND.  kk == 1 )        &
     2780                         READ ( 13 )  surf_v(3)%qcs
     2781                   CASE ( 'surf_v(3)%ncs' )         
     2782                      IF ( ALLOCATED( surf_v(3)%ncs )  .AND.  kk == 1 )        &
     2783                         READ ( 13 )  surf_v(3)%ncs
    25892784                   CASE ( 'surf_v(3)%qrs' )         
    25902785                      IF ( ALLOCATED( surf_v(3)%qrs )  .AND.  kk == 1 )        &
     
    26082803                      IF ( ALLOCATED( surf_v(3)%ssws )  .AND.  kk == 1 )       &
    26092804                         READ ( 13 )  surf_v(3)%ssws
     2805                   CASE ( 'surf_v(3)%qcsws' )         
     2806                      IF ( ALLOCATED( surf_v(3)%qcsws )  .AND.  kk == 1 )      &
     2807                         READ ( 13 )  surf_v(3)%qcsws
     2808                   CASE ( 'surf_v(3)%ncsws' )         
     2809                      IF ( ALLOCATED( surf_v(3)%ncsws )  .AND.  kk == 1 )      &
     2810                         READ ( 13 )  surf_v(3)%ncsws
    26102811                   CASE ( 'surf_v(3)%qrsws' )         
    26112812                      IF ( ALLOCATED( surf_v(3)%qrsws )  .AND.  kk == 1 )      &
     
    28143015             ENDIF
    28153016
     3017             IF ( SCAN( TRIM( field_chr ), '%qcs' ) /= 0 )  THEN
     3018                IF ( ALLOCATED( surf_target%qcs )  .AND.                       &
     3019                     ALLOCATED( surf_file%qcs   ) )                            &
     3020                  surf_target%qcs(m_target) = surf_file%qcs(m_file)
     3021             ENDIF
     3022
     3023             IF ( SCAN( TRIM( field_chr ), '%qcsws' ) /= 0 )  THEN
     3024                IF ( ALLOCATED( surf_target%qcsws )  .AND.                     &
     3025                     ALLOCATED( surf_file%qcsws   ) )                          &
     3026                   surf_target%qcsws(m_target) = surf_file%qcsws(m_file)
     3027             ENDIF
     3028
     3029             IF ( SCAN( TRIM( field_chr ), '%ncs' ) /= 0 )  THEN
     3030                IF ( ALLOCATED( surf_target%ncs )  .AND.                       &
     3031                     ALLOCATED( surf_file%ncs   ) )                            &
     3032                   surf_target%ncs(m_target) = surf_file%ncs(m_file)
     3033             ENDIF
     3034
     3035             IF ( SCAN( TRIM( field_chr ), '%ncsws' ) /= 0 )  THEN
     3036                IF ( ALLOCATED( surf_target%ncsws )  .AND.                     &
     3037                     ALLOCATED( surf_file%ncsws   ) )                          &
     3038                   surf_target%ncsws(m_target) = surf_file%ncsws(m_file)
     3039             ENDIF
     3040
    28163041             IF ( SCAN( TRIM( field_chr ), '%qrs' ) /= 0 )  THEN
    28173042                IF ( ALLOCATED( surf_target%qrs )  .AND.                       &
  • palm/trunk/SOURCE/swap_timelevel.f90

    r2233 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2233 2017-05-30 18:08:54Z suehring
    2732!
    2833! 2232 2017-05-30 17:47:52Z suehring
     
    103108#if defined( __nopointer )
    104109    USE arrays_3d,                                                             &
    105         ONLY:  e, e_p, nr, nr_p, pt, pt_p, q, q_p, qr, qr_p, s, s_p, sa, sa_p, &
    106                u, u_p, v, v_p, w, w_p
     110        ONLY:  e, e_p, nc, nc_p, nr, nr_p, pt, pt_p, q, q_p, qc, qc_p qr, qr_p,&
     111               s, s_p, sa, sa_p, u, u_p, v, v_p, w, w_p
    107112#else
    108113    USE arrays_3d,                                                             &
    109         ONLY:  e, e_1, e_2, e_p, nr, nr_1, nr_2, nr_p, pt, pt_1, pt_2, pt_p, q,&
    110                q_1, q_2, q_p, qr, qr_1, qr_2, qr_p, s, s_1, s_2, s_p, sa, sa_1,&
    111                sa_2, sa_p, u, u_1, u_2, u_p, v, v_1, v_2, v_p, w, w_1, w_2, w_p
     114        ONLY:  e, e_1, e_2, e_p, nc, nc_1, nc_2, nc_p, nr, nr_1, nr_2, nr_p,   &
     115               pt, pt_1, pt_2, pt_p, q, q_1, q_2, q_p, qc, qc_1, qc_2, qc_p,   &
     116               qr, qr_1, qr_2, qr_p, s, s_1, s_2, s_p, sa, sa_1, sa_2, sa_p,   &
     117               u, u_1, u_2, u_p, v, v_1, v_2, v_p, w, w_1, w_2, w_p
    112118
    113119#endif
     
    121127    USE control_parameters,                                                    &
    122128        ONLY:  cloud_physics, constant_diffusion, humidity, land_surface,      &
    123                microphysics_seifert, neutral, ocean, passive_scalar,           &
    124                timestep_count, urban_surface
     129               microphysics_morrison, microphysics_seifert, neutral, ocean,    &
     130               passive_scalar, timestep_count, urban_surface
    125131
    126132    USE indices,                                                               &
     
    174180
    175181    IF ( humidity )  THEN
    176        q = q_p             
     182       q = q_p
     183       IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     184          qc = qc_p
     185          nc = nc_p
     186       ENDIF             
    177187       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    178188          qr = qr_p
     
    214224          IF ( humidity )  THEN
    215225             q => q_1;    q_p => q_2
     226             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     227                qc => qc_1;    qc_p => qc_2
     228                nc => nc_1;    nc_p => nc_2
     229             ENDIF
    216230             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    217231                qr => qr_1;    qr_p => qr_2
     
    241255          IF ( humidity )  THEN
    242256             q => q_2;    q_p => q_1
     257             IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     258                qc => qc_2;    qc_p => qc_1
     259                nc => nc_2;    nc_p => nc_1
     260             ENDIF
    243261             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    244262                qr => qr_2;    qr_p => qr_1
  • palm/trunk/SOURCE/time_integration.f90

    r2271 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2271 2017-06-09 12:34:55Z sward
    2732! Start timestep message changed
    2833!
     
    267272
    268273    USE arrays_3d,                                                             &
    269         ONLY:  diss, dzu, e, e_p, nr, nr_p, prho, pt, pt_p, pt_init, q_init, q,&
    270                ql, ql_c, ql_v, ql_vp, qr, qr_p, q_p, ref_state, rho_ocean,     &
    271                s, s_p, sa_p, tend, u, u_p, v, vpt, v_p, w, w_p
     274        ONLY:  diss, dzu, e, e_p, nc, nc_p, nr, nr_p, prho, pt, pt_p, pt_init, &
     275               q_init, q, qc, qc_p, ql, ql_c, ql_v, ql_vp, qr, qr_p, q_p,      &
     276               ref_state, rho_ocean, s, s_p, sa_p, tend, u, u_p, v, vpt,       &
     277               v_p, w, w_p
    272278
    273279    USE calc_mean_profile_mod,                                                 &
     
    291297               land_surface, large_scale_forcing,                              &
    292298               loop_optimization, lsf_surf, lsf_vert, masks,                   &
    293                microphysics_seifert, mid, nest_domain,                         &
     299               microphysics_morrison, microphysics_seifert, mid, nest_domain,  &
    294300               neutral, nr_timesteps_this_run, nudging,                        &
    295301               ocean, passive_scalar,                                          &
     
    585591          IF ( humidity )  THEN
    586592             CALL exchange_horiz( q_p, nbgp )
     593             IF ( cloud_physics .AND. microphysics_morrison )  THEN
     594                CALL exchange_horiz( qc_p, nbgp )
     595                CALL exchange_horiz( nc_p, nbgp )
     596             ENDIF
    587597             IF ( cloud_physics .AND. microphysics_seifert )  THEN
    588598                CALL exchange_horiz( qr_p, nbgp )
     
    637647                   CALL exchange_horiz( q, nbgp )
    638648
     649                   IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
     650                       CALL exchange_horiz( qc, nbgp )
     651                       CALL exchange_horiz( nc, nbgp )
     652                   ENDIF
    639653                   IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
    640 !                        CALL exchange_horiz( qc, nbgp )
    641654                       CALL exchange_horiz( qr, nbgp )
    642 !                        CALL exchange_horiz( nc, nbgp )
    643655                       CALL exchange_horiz( nr, nbgp )
    644656                   ENDIF
  • palm/trunk/SOURCE/write_3d_binary.f90

    r2233 r2292  
    2525! -----------------
    2626! $Id$
     27! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
     28! includes two more prognostic equations for cloud drop concentration (nc) 
     29! and cloud water content (qc).
     30!
     31! 2233 2017-05-30 18:08:54Z suehring
    2732!
    2833! 2232 2017-05-30 17:47:52Z suehring
     
    116121
    117122    USE arrays_3d,                                                             &
    118         ONLY:  e, kh, km, p, pt, q, ql, qc, nr, prr, precipitation_amount, qr, &
    119                s, sa, u, u_m_l, u_m_n, u_m_r, u_m_s, v, v_m_l, v_m_n, v_m_r,   &
    120                v_m_s, vpt, w, w_m_l, w_m_n, w_m_r, w_m_s
     123        ONLY:  e, kh, km, p, pt, q, ql, qc, nc, nr, prr, precipitation_amount, &
     124               qr, s, sa, u, u_m_l, u_m_n, u_m_r, u_m_s, v, v_m_l, v_m_n,      &
     125               v_m_r, v_m_s, vpt, w, w_m_l, w_m_n, w_m_r, w_m_s
    121126       
    122127    USE averaging
     
    124129    USE control_parameters,                                                    &
    125130        ONLY:  iran, humidity, passive_scalar, cloud_physics, cloud_droplets,  &
    126                microphysics_seifert, ocean, topography
     131               microphysics_morrison, microphysics_seifert, ocean, topography
    127132               
    128133    USE indices,                                                               &
     
    222227             IF ( ALLOCATED( qc_av ) )  THEN
    223228                WRITE ( 14 )  'qc_av               ';  WRITE ( 14 )  qc_av
     229             ENDIF
     230             IF ( microphysics_morrison )  THEN
     231                WRITE ( 14 )  'nc                  ';  WRITE ( 14 )  nc
     232                IF ( ALLOCATED( nc_av ) )  THEN
     233                   WRITE ( 14 )  'nc_av               ';  WRITE ( 14 )  nc_av
     234                ENDIF
    224235             ENDIF
    225236             IF ( microphysics_seifert )  THEN
Note: See TracChangeset for help on using the changeset viewer.