Changeset 2155 for palm/trunk


Ignore:
Timestamp:
Feb 21, 2017 9:57:40 AM (8 years ago)
Author:
hoffmann
Message:

bugfix in cloud microphysics

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

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

    r2101 r2155  
    2020! Current revisions:
    2121! ------------------
    22 !
    23 ! 
     22! Bugfix in the calculation of microphysical quantities on ghost points.
     23!
    2424! Former revisions:
    2525! -----------------
     
    2828! 2031 2016-10-21 15:11:58Z knoop
    2929! renamed variable rho to rho_ocean
    30 ! 
     30!
    3131! 2000 2016-08-20 18:09:15Z knoop
    3232! Forced header and separation lines into 80 columns
    33 ! 
     33!
    3434! 1850 2016-04-08 13:29:27Z maronga
    3535! Module renamed
     
    5252! Added new routine calc_precipitation_amount. The routine now allows to account
    5353! for precipitation due to sedimenation of cloud (fog) droplets
    54 ! 
     54!
    5555! 1682 2015-10-07 23:56:08Z knoop
    56 ! Code annotations made doxygen readable 
     56! Code annotations made doxygen readable
    5757!
    5858! 1646 2015-09-02 16:00:10Z hoffmann
     
    6363! Vectorized version of adjust_cloud added.
    6464! Little reformatting of the code.
    65 ! 
     65!
    6666! 1353 2014-04-08 15:21:23Z heinze
    6767! REAL constants provided with KIND-attribute
    68 ! 
     68!
    6969! 1346 2014-03-27 13:18:20Z heinze
    70 ! Bugfix: REAL constants provided with KIND-attribute especially in call of 
     70! Bugfix: REAL constants provided with KIND-attribute especially in call of
    7171! intrinsic function like MAX, MIN, SIGN
    72 ! 
     72!
    7373! 1334 2014-03-25 12:21:40Z heinze
    7474! Bugfix: REAL constants provided with KIND-attribute
     
    106106! 1053 2012-11-13 17:11:03Z hoffmann
    107107! initial revision
    108 ! 
     108!
    109109! Description:
    110110! ------------
     
    140140    REAL(wp) ::  c_term = 600.0_wp         !< coef. for terminal velocity (m-1)
    141141    REAL(wp) ::  diff_coeff_l = 0.23E-4_wp  !< diffusivity of water vapor (m2 s-1)
    142     REAL(wp) ::  eps_sb = 1.0E-20_wp       !< threshold in two-moments scheme
     142    REAL(wp) ::  eps_sb = 1.0E-10_wp       !< threshold in two-moments scheme
    143143    REAL(wp) ::  k_cc = 9.44E09_wp         !< const. cloud-cloud kernel (m3 kg-2 s-1)
    144144    REAL(wp) ::  k_cr0 = 4.33_wp           !< const. cloud-rain kernel (m3 kg-1 s-1)
     
    221221       MODULE PROCEDURE sedimentation_cloud_ij
    222222    END INTERFACE sedimentation_cloud
    223  
     223
    224224    INTERFACE sedimentation_rain
    225225       MODULE PROCEDURE sedimentation_rain
     
    337337             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
    338338             t_d_pt(k) = 1.0_wp / pt_d_t(k)
    339              hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
     339             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
    340340          ENDDO
    341341
     
    346346
    347347!
    348 !--    Compute length of time step 
     348!--    Compute length of time step
    349349       IF ( call_microphysics_at_all_substeps )  THEN
    350350          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
     
    383383! Description:
    384384! ------------
    385 !> Adjust number of raindrops to avoid nonlinear effects in sedimentation and 
    386 !> evaporation of rain drops due to too small or too big weights 
     385!> Adjust number of raindrops to avoid nonlinear effects in sedimentation and
     386!> evaporation of rain drops due to too small or too big weights
    387387!> of rain drops (Stevens and Seifert, 2008).
    388388!------------------------------------------------------------------------------!
     
    399399
    400400       USE indices,                                                            &
    401            ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
     401           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
    402402
    403403       USE kinds
     
    411411       CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' )
    412412
    413        DO  i = nxl, nxr
    414           DO  j = nys, nyn
     413       DO  i = nxlg, nxrg
     414          DO  j = nysg, nyng
    415415             DO  k = nzb_s_inner(j,i)+1, nzt
    416416                IF ( qr(k,j,i) <= eps_sb )  THEN
     
    456456
    457457       USE indices,                                                            &
    458            ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
     458           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
    459459
    460460       USE kinds
     
    466466       INTEGER(iwp) ::  k                 !<
    467467
    468        REAL(wp)     ::  alpha_cc          !<                   
     468       REAL(wp)     ::  alpha_cc          !<
    469469       REAL(wp)     ::  autocon           !<
    470470       REAL(wp)     ::  dissipation       !<
     
    482482       CALL cpu_log( log_point_s(55), 'autoconversion', 'start' )
    483483
    484        DO  i = nxl, nxr
    485           DO  j = nys, nyn
     484       DO  i = nxlg, nxrg
     485          DO  j = nysg, nyng
    486486             DO  k = nzb_s_inner(j,i)+1, nzt
    487487
     
    494494                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) )
    495495!
    496 !--                Universal function for autoconversion process 
     496!--                Universal function for autoconversion process
    497497!--                (Seifert and Beheng, 2006):
    498498                   phi_au = 600.0_wp * tau_cloud**0.68_wp *                    &
     
    506506                   xc = hyrho(k) * qc(k,j,i) / nc_const
    507507!
    508 !--                Parameterized turbulence effects on autoconversion (Seifert, 
     508!--                Parameterized turbulence effects on autoconversion (Seifert,
    509509!--                Nuijens and Stevens, 2010)
    510510                   IF ( collision_turbulence )  THEN
     
    517517                      sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )
    518518!
    519 !--                   Mixing length (neglecting distance to ground and 
     519!--                   Mixing length (neglecting distance to ground and
    520520!--                   stratification)
    521521                      l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
    522522!
    523 !--                   Limit dissipation rate according to Seifert, Nuijens and 
     523!--                   Limit dissipation rate according to Seifert, Nuijens and
    524524!--                   Stevens (2010)
    525525                      dissipation = MIN( 0.06_wp, diss(k,j,i) )
     
    531531                                  dissipation**( 1.0_wp / 6.0_wp )
    532532!
    533 !--                   The factor of 1.0E4 is needed to convert the dissipation 
     533!--                   The factor of 1.0E4 is needed to convert the dissipation
    534534!--                   rate from m2 s-3 to cm2 s-3.
    535535                      k_au = k_au * ( 1.0_wp +                                 &
    536536                                      dissipation * 1.0E4_wp *                 &
    537537                                      ( re_lambda * 1.0E-3_wp )**0.25_wp *     &
    538                                       ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &     
     538                                      ( alpha_cc * EXP( -1.0_wp * ( ( rc -     &
    539539                                                                      r_cc ) / &
    540540                                                        sigma_cc )**2          &
     
    552552
    553553                   qr(k,j,i) = qr(k,j,i) + autocon * dt_micro
    554                    qc(k,j,i) = qc(k,j,i) - autocon * dt_micro 
     554                   qc(k,j,i) = qc(k,j,i) - autocon * dt_micro
    555555                   nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro
    556556
     
    580580
    581581       USE indices,                                                            &
    582            ONLY:  nxl, nxr, nyn, nys, nzb_s_inner, nzt
     582           ONLY:  nxlg, nxrg, nyng, nysg, nzb_s_inner, nzt
    583583
    584584       USE kinds
     
    593593       REAL(wp)    ::  dqdt_precip !<
    594594
    595        DO  i = nxl, nxr
    596           DO  j = nys, nyn
     595       DO  i = nxlg, nxrg
     596          DO  j = nysg, nyng
    597597             DO  k = nzb_s_inner(j,i)+1, nzt
    598598
     
    640640
    641641       USE indices,                                                            &
    642            ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
     642           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
    643643
    644644       USE kinds
     
    657657       CALL cpu_log( log_point_s(56), 'accretion', 'start' )
    658658
    659        DO  i = nxl, nxr
    660           DO  j = nys, nyn
     659       DO  i = nxlg, nxrg
     660          DO  j = nysg, nyng
    661661             DO  k = nzb_s_inner(j,i)+1, nzt
    662662
     
    664664!
    665665!--                Intern time scale of coagulation (Seifert and Beheng, 2006):
    666                    tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ) 
    667 !
    668 !--                Universal function for accretion process (Seifert and 
     666                   tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )
     667!
     668!--                Universal function for accretion process (Seifert and
    669669!--                Beheng, 2001):
    670670                   phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
    671671!
    672 !--                Parameterized turbulence effects on autoconversion (Seifert, 
    673 !--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 
     672!--                Parameterized turbulence effects on autoconversion (Seifert,
     673!--                Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
    674674!--                convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
    675675                   IF ( collision_turbulence )  THEN
     
    679679                                     )
    680680                   ELSE
    681                       k_cr = k_cr0                       
     681                      k_cr = k_cr0
    682682                   ENDIF
    683683!
     
    687687                   accr = MIN( accr, qc(k,j,i) / dt_micro )
    688688
    689                    qr(k,j,i) = qr(k,j,i) + accr * dt_micro 
     689                   qr(k,j,i) = qr(k,j,i) + accr * dt_micro
    690690                   qc(k,j,i) = qc(k,j,i) - accr * dt_micro
    691691
     
    721721
    722722       USE indices,                                                            &
    723            ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
     723           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
    724724
    725725       USE kinds
    726    
     726
    727727       IMPLICIT NONE
    728728
     
    738738       CALL cpu_log( log_point_s(57), 'selfcollection', 'start' )
    739739
    740        DO  i = nxl, nxr
    741           DO  j = nys, nyn
     740       DO  i = nxlg, nxrg
     741          DO  j = nysg, nyng
    742742             DO  k = nzb_s_inner(j,i)+1, nzt
    743743                IF ( qr(k,j,i) > eps_sb )  THEN
     
    762762                   nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro
    763763
    764                 ENDIF         
     764                ENDIF
    765765             ENDDO
    766766          ENDDO
     
    775775! Description:
    776776! ------------
    777 !> Evaporation of precipitable water. Condensation is neglected for 
     777!> Evaporation of precipitable water. Condensation is neglected for
    778778!> precipitable water.
    779779!------------------------------------------------------------------------------!
     
    793793
    794794       USE indices,                                                            &
    795            ONLY:  nxl, nxr, nys, nyn, nzb_s_inner, nzt
     795           ONLY:  nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt
    796796
    797797       USE kinds
     
    823823       CALL cpu_log( log_point_s(58), 'evaporation', 'start' )
    824824
    825        DO  i = nxl, nxr
    826           DO  j = nys, nyn
     825       DO  i = nxlg, nxrg
     826          DO  j = nysg, nyng
    827827             DO  k = nzb_s_inner(j,i)+1, nzt
    828828                IF ( qr(k,j,i) > eps_sb )  THEN
     
    850850!--                   Actual temperature:
    851851                      temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) )
    852              
     852
    853853                      g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) *     &
    854854                                          l_v / ( thermal_conductivity_l * temp ) &
     
    862862                      dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
    863863!
    864 !--                   Compute ventilation factor and intercept parameter 
     864!--                   Compute ventilation factor and intercept parameter
    865865!--                   (Seifert and Beheng, 2006; Seifert, 2008):
    866866                      IF ( ventilation_effect )  THEN
    867867!
    868 !--                      Shape parameter of gamma distribution (Milbrandt and Yau, 
     868!--                      Shape parameter of gamma distribution (Milbrandt and Yau,
    869869!--                      2005; Stevens and Seifert, 2008):
    870870                         mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp *          &
     
    877877
    878878                         mu_r_2   = mu_r + 2.0_wp
    879                          mu_r_5d2 = mu_r + 2.5_wp 
     879                         mu_r_5d2 = mu_r + 2.5_wp
    880880
    881881                         f_vent = a_vent * gamm( mu_r_2 ) *                     &
     
    893893                                    ( lambda_r / ( 3.0_wp * c_term + lambda_r ) &
    894894                                    )**mu_r_5d2 -                               &
    895                                     0.0390625_wp * ( b_term / a_term )**4 *     & 
     895                                    0.0390625_wp * ( b_term / a_term )**4 *     &
    896896                                    ( lambda_r / ( 4.0_wp * c_term + lambda_r ) &
    897897                                    )**mu_r_5d2                                 &
     
    899899
    900900                         nr_0   = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) /    &
    901                                   gamm( mu_r + 1.0_wp ) 
     901                                  gamm( mu_r + 1.0_wp )
    902902                      ELSE
    903903                         f_vent = 1.0_wp
     
    905905                      ENDIF
    906906   !
    907    !--                Evaporation rate of rain water content (Seifert and 
     907   !--                Evaporation rate of rain water content (Seifert and
    908908   !--                Beheng, 2006):
    909909                      evap    = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat /   &
     
    917917
    918918                   ENDIF
    919                 ENDIF         
     919                ENDIF
    920920
    921921             ENDDO
     
    948948
    949949       USE indices,                                                            &
    950            ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
     950           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzb_s_inner, nzt
    951951
    952952       USE kinds
     
    954954       USE statistics,                                                         &
    955955           ONLY:  weight_substep
    956    
     956
    957957
    958958       IMPLICIT NONE
     
    968968       sed_qc(nzt+1) = 0.0_wp
    969969
    970        DO  i = nxl, nxr
    971           DO  j = nys, nyn
     970       DO  i = nxlg, nxrg
     971          DO  j = nysg, nyng
    972972             DO  k = nzt, nzb_s_inner(j,i)+1, -1
    973973
     
    10291029
    10301030       USE indices,                                                            &
    1031            ONLY:  nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt
     1031           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzb_s_inner, nzt
    10321032
    10331033       USE kinds
     
    10351035       USE statistics,                                                         &
    10361036           ONLY:  weight_substep
    1037        
     1037
    10381038       IMPLICIT NONE
    10391039
     
    10651065
    10661066!
    1067 !--    Compute velocities 
    1068        DO  i = nxl, nxr
    1069           DO  j = nys, nyn
     1067!--    Compute velocities
     1068       DO  i = nxlg, nxrg
     1069          DO  j = nysg, nyng
    10701070             DO  k = nzb_s_inner(j,i)+1, nzt
    10711071                IF ( qr(k,j,i) > eps_sb )  THEN
     
    11191119                                      2.0_wp * w_qr(k) + w_qr(k+1) ) *         &
    11201120                          dt_micro * ddzu(k)
    1121              ENDDO     
     1121             ENDDO
    11221122!
    11231123!--          Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
     
    11551155             DO  k = nzt, nzb_s_inner(j,i)+1, -1
    11561156!
    1157 !--             Sum up all rain drop number densities which contribute to the flux 
     1157!--             Sum up all rain drop number densities which contribute to the flux
    11581158!--             through k-1/2
    11591159                flux  = 0.0_wp
     
    11701170                ENDDO
    11711171!
    1172 !--             It is not allowed to sediment more rain drop number density than 
     1172!--             It is not allowed to sediment more rain drop number density than
    11731173!--             available
    11741174                flux = MIN( flux,                                              &
     
    11811181                                        ddzu(k+1) / hyrho(k) * dt_micro
    11821182!
    1183 !--             Sum up all rain water content which contributes to the flux 
     1183!--             Sum up all rain water content which contributes to the flux
    11841184!--             through k-1/2
    11851185                flux  = 0.0_wp
     
    11991199                ENDDO
    12001200!
    1201 !--             It is not allowed to sediment more rain water content than 
     1201!--             It is not allowed to sediment more rain water content than
    12021202!--             available
    12031203                flux = MIN( flux,                                              &
     
    12111211                                        ddzu(k+1) / hyrho(k) * dt_micro
    12121212                q(k,j,i)  = q(k,j,i)  + ( sed_qr(k+1) - sed_qr(k) ) *          &
    1213                                         ddzu(k+1) / hyrho(k) * dt_micro 
     1213                                        ddzu(k+1) / hyrho(k) * dt_micro
    12141214                pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) *          &
    12151215                                        ddzu(k+1) / hyrho(k) * l_d_cp *        &
     
    13301330             pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp
    13311331             t_d_pt(k) = 1.0_wp / pt_d_t(k)
    1332              hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
     1332             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )
    13331333          ENDDO
    13341334!
     
    13381338
    13391339!
    1340 !--    Compute length of time step 
     1340!--    Compute length of time step
    13411341       IF ( call_microphysics_at_all_substeps )  THEN
    13421342          dt_micro = dt_3d * weight_pres(intermediate_timestep_count)
     
    13951395! Description:
    13961396! ------------
    1397 !> Adjust number of raindrops to avoid nonlinear effects in 
     1397!> Adjust number of raindrops to avoid nonlinear effects in
    13981398!> sedimentation and evaporation of rain drops due to too small or
    13991399!> too big weights of rain drops (Stevens and Seifert, 2008).
     
    14241424          ELSE
    14251425!
    1426 !--          Adjust number of raindrops to avoid nonlinear effects in 
     1426!--          Adjust number of raindrops to avoid nonlinear effects in
    14271427!--          sedimentation and evaporation of rain drops due to too small or
    14281428!--          too big weights of rain drops (Stevens and Seifert, 2008).
     
    14701470       INTEGER(iwp) ::  k                 !<
    14711471
    1472        REAL(wp)     ::  alpha_cc          !<                   
     1472       REAL(wp)     ::  alpha_cc          !<
    14731473       REAL(wp)     ::  autocon           !<
    14741474       REAL(wp)     ::  dissipation       !<
     
    14941494             tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )
    14951495!
    1496 !--          Universal function for autoconversion process 
     1496!--          Universal function for autoconversion process
    14971497!--          (Seifert and Beheng, 2006):
    14981498             phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3
     
    15051505             xc = hyrho(k) * qc_1d(k) / nc_1d(k)
    15061506!
    1507 !--          Parameterized turbulence effects on autoconversion (Seifert, 
     1507!--          Parameterized turbulence effects on autoconversion (Seifert,
    15081508!--          Nuijens and Stevens, 2010)
    15091509             IF ( collision_turbulence )  THEN
     
    15191519                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
    15201520!
    1521 !--             Limit dissipation rate according to Seifert, Nuijens and 
     1521!--             Limit dissipation rate according to Seifert, Nuijens and
    15221522!--             Stevens (2010)
    15231523                dissipation = MIN( 0.06_wp, diss(k,j,i) )
     
    15491549
    15501550             qr_1d(k) = qr_1d(k) + autocon * dt_micro
    1551              qc_1d(k) = qc_1d(k) - autocon * dt_micro 
     1551             qc_1d(k) = qc_1d(k) - autocon * dt_micro
    15521552             nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro
    15531553
     
    16421642!
    16431643!--          Intern time scale of coagulation (Seifert and Beheng, 2006):
    1644              tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 
    1645 !
    1646 !--          Universal function for accretion process 
     1644             tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )
     1645!
     1646!--          Universal function for accretion process
    16471647!--          (Seifert and Beheng, 2001):
    16481648             phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4
    16491649!
    1650 !--          Parameterized turbulence effects on autoconversion (Seifert, 
    1651 !--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 
     1650!--          Parameterized turbulence effects on autoconversion (Seifert,
     1651!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
    16521652!--          convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.
    16531653             IF ( collision_turbulence )  THEN
     
    16571657                               )
    16581658             ELSE
    1659                 k_cr = k_cr0                       
     1659                k_cr = k_cr0
    16601660             ENDIF
    16611661!
     
    16641664             accr = MIN( accr, qc_1d(k) / dt_micro )
    16651665
    1666              qr_1d(k) = qr_1d(k) + accr * dt_micro 
     1666             qr_1d(k) = qr_1d(k) + accr * dt_micro
    16671667             qc_1d(k) = qc_1d(k) - accr * dt_micro
    16681668
     
    16911691
    16921692       USE kinds
    1693    
     1693
    16941694       IMPLICIT NONE
    16951695
     
    17231723             nr_1d(k) = nr_1d(k) + selfcoll * dt_micro
    17241724
    1725           ENDIF         
     1725          ENDIF
    17261726       ENDDO
    17271727
     
    17321732! Description:
    17331733! ------------
    1734 !> Evaporation of precipitable water. Condensation is neglected for 
     1734!> Evaporation of precipitable water. Condensation is neglected for
    17351735!> precipitable water. Call for grid point i,j
    17361736!------------------------------------------------------------------------------!
     
    17991799!--             Actual temperature:
    18001800                temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) )
    1801        
     1801
    18021802                g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v /  &
    18031803                                    ( thermal_conductivity_l * temp ) +        &
     
    18111811                dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
    18121812!
    1813 !--             Compute ventilation factor and intercept parameter 
     1813!--             Compute ventilation factor and intercept parameter
    18141814!--             (Seifert and Beheng, 2006; Seifert, 2008):
    18151815                IF ( ventilation_effect )  THEN
     
    18251825
    18261826                   mu_r_2   = mu_r + 2.0_wp
    1827                    mu_r_5d2 = mu_r + 2.5_wp 
    1828 
    1829                    f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  & 
     1827                   mu_r_5d2 = mu_r + 2.5_wp
     1828
     1829                   f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) +  &
    18301830                            b_vent * schmidt_p_1d3 *                           &
    18311831                            SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) *  &
     
    18411841                              ( lambda_r / ( 3.0_wp * c_term + lambda_r )      &
    18421842                              )**mu_r_5d2 -                                    &
    1843                               0.0390625_wp * ( b_term / a_term )**4 *          & 
     1843                              0.0390625_wp * ( b_term / a_term )**4 *          &
    18441844                              ( lambda_r / ( 4.0_wp * c_term + lambda_r )      &
    18451845                              )**mu_r_5d2                                      &
     
    18471847
    18481848                   nr_0   = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) /           &
    1849                             gamm( mu_r + 1.0_wp ) 
     1849                            gamm( mu_r + 1.0_wp )
    18501850                ELSE
    18511851                   f_vent = 1.0_wp
     
    18631863
    18641864             ENDIF
    1865           ENDIF         
     1865          ENDIF
    18661866
    18671867       ENDDO
     
    18731873! Description:
    18741874! ------------
    1875 !> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR). 
     1875!> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).
    18761876!> Call for grid point i,j
    18771877!------------------------------------------------------------------------------!
     
    18911891
    18921892       USE kinds
    1893        
     1893
    18941894       USE statistics,                                                         &
    18951895           ONLY:  weight_substep
     
    19191919          q_1d(k)  = q_1d(k)  + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
    19201920                                hyrho(k) * dt_micro
    1921           qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      & 
     1921          qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
    19221922                                hyrho(k) * dt_micro
    19231923          pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) /      &
     
    19621962       USE statistics,                                                         &
    19631963           ONLY:  weight_substep
    1964        
     1964
    19651965       IMPLICIT NONE
    19661966
     
    19901990
    19911991!
    1992 !--    Compute velocities 
     1992!--    Compute velocities
    19931993       DO  k = nzb_s_inner(j,i)+1, nzt
    19941994          IF ( qr_1d(k) > eps_sb )  THEN
     
    20352035          c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) *   &
    20362036                    dt_micro * ddzu(k)
    2037        ENDDO     
     2037       ENDDO
    20382038!
    20392039!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
     
    20712071       DO  k = nzt, nzb_s_inner(j,i)+1, -1
    20722072!
    2073 !--       Sum up all rain drop number densities which contribute to the flux 
     2073!--       Sum up all rain drop number densities which contribute to the flux
    20742074!--       through k-1/2
    20752075          flux  = 0.0_wp
     
    20862086          ENDDO
    20872087!
    2088 !--       It is not allowed to sediment more rain drop number density than 
     2088!--       It is not allowed to sediment more rain drop number density than
    20892089!--       available
    20902090          flux = MIN( flux,                                                    &
     
    20952095                                    hyrho(k) * dt_micro
    20962096!
    2097 !--       Sum up all rain water content which contributes to the flux 
     2097!--       Sum up all rain water content which contributes to the flux
    20982098!--       through k-1/2
    20992099          flux  = 0.0_wp
     
    21222122                                hyrho(k) * dt_micro
    21232123          q_1d(k)  = q_1d(k)  + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
    2124                                 hyrho(k) * dt_micro 
     2124                                hyrho(k) * dt_micro
    21252125          pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /      &
    21262126                                hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro
     
    21852185! ------------
    21862186!> This function computes the gamma function (Press et al., 1992).
    2187 !> The gamma function is needed for the calculation of the evaporation 
     2187!> The gamma function is needed for the calculation of the evaporation
    21882188!> of rain drops.
    21892189!------------------------------------------------------------------------------!
    2190     FUNCTION gamm( xx ) 
     2190    FUNCTION gamm( xx )
    21912191
    21922192       USE kinds
    21932193
    2194        IMPLICIT NONE 
     2194       IMPLICIT NONE
    21952195
    21962196       INTEGER(iwp) ::  j            !<
     
    22182218       ser = 1.000000000190015_wp
    22192219
    2220        DO  j = 1, 6 
    2221           y_gamm = y_gamm + 1.0_wp 
    2222           ser    = ser + cof( j ) / y_gamm 
     2220       DO  j = 1, 6
     2221          y_gamm = y_gamm + 1.0_wp
     2222          ser    = ser + cof( j ) / y_gamm
    22232223       ENDDO
    22242224
    2225 ! 
    2226 !--    Until this point the algorithm computes the logarithm of the gamma 
    2227 !--    function. Hence, the exponential function is used. 
    2228 !       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) ) 
    2229        gamm = EXP( tmp ) * stp * ser / x_gamm 
    2230 
    2231        RETURN 
    2232  
    2233     END FUNCTION gamm 
     2225!
     2226!--    Until this point the algorithm computes the logarithm of the gamma
     2227!--    function. Hence, the exponential function is used.
     2228!       gamm = EXP( tmp + LOG( stp * ser / x_gamm ) )
     2229       gamm = EXP( tmp ) * stp * ser / x_gamm
     2230
     2231       RETURN
     2232
     2233    END FUNCTION gamm
    22342234
    22352235 END MODULE microphysics_mod
  • palm/trunk/SOURCE/prognostic_equations.f90

    r2119 r2155  
    2020! Current revisions:
    2121! ------------------
    22 !
    23 ! 
     22! Bugfix in the calculation of microphysical quantities on ghost points.
     23!
    2424! Former revisions:
    2525! -----------------
     
    2828! 2118 2017-01-17 16:38:49Z raasch
    2929! OpenACC version of subroutine removed
    30 ! 
     30!
    3131! 2031 2016-10-21 15:11:58Z knoop
    3232! renamed variable rho to rho_ocean
    33 ! 
     33!
    3434! 2011 2016-09-19 17:29:57Z kanani
    3535! Flag urban_surface is now defined in module control_parameters.
    36 ! 
     36!
    3737! 2007 2016-08-24 15:47:17Z kanani
    3838! Added pt tendency calculation based on energy balance at urban surfaces
    3939! (new urban surface model)
    40 ! 
     40!
    4141! 2000 2016-08-20 18:09:15Z knoop
    4242! Forced header and separation lines into 80 columns
    43 ! 
     43!
    4444! 1976 2016-07-27 13:28:04Z maronga
    4545! Simplied calls to radiation model
    46 ! 
     46!
    4747! 1960 2016-07-12 16:34:24Z suehring
    4848! Separate humidity and passive scalar
    49 ! 
     49!
    5050! 1914 2016-05-26 14:44:07Z witha
    5151! Added calls for wind turbine model
     
    5353! 1873 2016-04-18 14:50:06Z maronga
    5454! Module renamed (removed _mod)
    55 ! 
     55!
    5656! 1850 2016-04-08 13:29:27Z maronga
    5757! Module renamed
    58 ! 
     58!
    5959! 1826 2016-04-07 12:01:39Z maronga
    6060! Renamed canopy model calls.
    61 ! 
     61!
    6262! 1822 2016-04-07 07:49:42Z hoffmann
    6363! Kessler microphysics scheme moved to microphysics.
    6464!
    6565! 1757 2016-02-22 15:49:32Z maronga
    66 ! 
     66!
    6767! 1691 2015-10-26 16:17:44Z maronga
    6868! Added optional model spin-up without radiation / land surface model calls.
    6969! Formatting corrections.
    70 ! 
     70!
    7171! 1682 2015-10-07 23:56:08Z knoop
    72 ! Code annotations made doxygen readable 
    73 ! 
     72! Code annotations made doxygen readable
     73!
    7474! 1585 2015-04-30 07:05:52Z maronga
    7575! Added call for temperature tendency calculation due to radiative flux divergence
    76 ! 
     76!
    7777! 1517 2015-01-07 19:12:25Z hoffmann
    7878! advec_s_bc_mod addded, since advec_s_bc is now a module
     
    8080! 1496 2014-12-02 17:25:50Z maronga
    8181! Renamed "radiation" -> "cloud_top_radiation"
    82 ! 
     82!
    8383! 1484 2014-10-21 10:53:05Z kanani
    8484! Changes due to new module structure of the plant canopy model:
     
    8686! Removed double-listing of use_upstream_for_tke in ONLY-list of module
    8787! control_parameters
    88 ! 
     88!
    8989! 1409 2014-05-23 12:11:32Z suehring
    90 ! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary. 
     90! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary.
    9191! This ensures that left-hand side fluxes are also calculated for nxl in that
    92 ! case, even though the solution at nxl is overwritten in boundary_conds() 
    93 ! 
     92! case, even though the solution at nxl is overwritten in boundary_conds()
     93!
    9494! 1398 2014-05-07 11:15:00Z heinze
    9595! Rayleigh-damping for horizontal velocity components changed: instead of damping
    96 ! against ug and vg, damping against u_init and v_init is used to allow for a 
     96! against ug and vg, damping against u_init and v_init is used to allow for a
    9797! homogenized treatment in case of nudging
    98 ! 
     98!
    9999! 1380 2014-04-28 12:40:45Z heinze
    100 ! Change order of calls for scalar prognostic quantities: 
    101 ! ls_advec -> nudging -> subsidence since initial profiles 
    102 ! 
     100! Change order of calls for scalar prognostic quantities:
     101! ls_advec -> nudging -> subsidence since initial profiles
     102!
    103103! 1374 2014-04-25 12:55:07Z raasch
    104104! missing variables added to ONLY lists
    105 ! 
     105!
    106106! 1365 2014-04-22 15:03:56Z boeske
    107 ! Calls of ls_advec for large scale advection added, 
     107! Calls of ls_advec for large scale advection added,
    108108! subroutine subsidence is only called if use_subsidence_tendencies = .F.,
    109109! new argument ls_index added to the calls of subsidence
    110110! +ls_index
    111 ! 
     111!
    112112! 1361 2014-04-16 15:17:48Z hoffmann
    113113! Two-moment microphysics moved to the start of prognostic equations. This makes
     
    117117!
    118118! Two-moment cloud physics added for vector and accelerator optimization.
    119 ! 
     119!
    120120! 1353 2014-04-08 15:21:23Z heinze
    121121! REAL constants provided with KIND-attribute
    122 ! 
     122!
    123123! 1337 2014-03-25 15:11:48Z heinze
    124124! Bugfix: REAL constants provided with KIND-attribute
    125 ! 
     125!
    126126! 1332 2014-03-25 11:59:43Z suehring
    127 ! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke 
    128 ! 
     127! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke
     128!
    129129! 1330 2014-03-24 17:29:32Z suehring
    130 ! In case of SGS-particle velocity advection of TKE is also allowed with 
     130! In case of SGS-particle velocity advection of TKE is also allowed with
    131131! dissipative 5th-order scheme.
    132132!
     
    161161!
    162162! 1115 2013-03-26 18:16:16Z hoffmann
    163 ! optimized cloud physics: calculation of microphysical tendencies transfered 
     163! optimized cloud physics: calculation of microphysical tendencies transfered
    164164! to microphysics.f90; qr and nr are only calculated if precipitation is required
    165165!
     
    211211 MODULE prognostic_equations_mod
    212212
    213  
     213
    214214
    215215    USE arrays_3d,                                                             &
     
    225225               te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m, tsa_m, tswst, tu_m, tv_m,&
    226226               tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p
    227        
     227
    228228    USE control_parameters,                                                    &
    229229        ONLY:  call_microphysics_at_all_substeps, cloud_physics,               &
     
    249249
    250250    USE indices,                                                               &
    251         ONLY:  nxl, nxlu, nxr, nyn, nys, nysv, nzb_s_inner, nzb_u_inner,       &
    252                nzb_v_inner, nzb_w_inner, nzt
     251        ONLY:  nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,         &
     252               nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
    253253
    254254    USE advec_ws,                                                              &
     
    287287    USE calc_radiation_mod,                                                    &
    288288        ONLY:  calc_radiation
    289  
     289
    290290    USE coriolis_mod,                                                          &
    291291        ONLY:  coriolis
     
    370370!> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.)
    371371!------------------------------------------------------------------------------!
    372  
     372
    373373 SUBROUTINE prognostic_equations_cache
    374374
     
    382382    INTEGER(iwp) ::  omp_get_thread_num  !<
    383383    INTEGER(iwp) ::  tn = 0              !<
    384    
     384
    385385    LOGICAL      ::  loop_start          !<
    386386
     
    394394!$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn)
    395395
    396 !$  tn = omp_get_thread_num() 
     396!$  tn = omp_get_thread_num()
    397397    loop_start = .TRUE.
    398398!$OMP DO
     399
     400!
     401!-- If required, calculate cloud microphysics
     402    IF ( cloud_physics  .AND.  .NOT. microphysics_sat_adjust  .AND.            &
     403         ( intermediate_timestep_count == 1  .OR.                              &
     404           call_microphysics_at_all_substeps )                                 &
     405       )  THEN
     406       DO  i = nxlg, nxrg
     407          DO  j = nysg, nyng
     408             CALL microphysics_control( i, j )
     409           END DO
     410       END DO
     411    ENDIF
     412
    399413    DO  i = nxl, nxr
    400414
     
    404418       IF ( loop_start )  THEN
    405419          loop_start  = .FALSE.
    406           i_omp_start = i 
     420          i_omp_start = i
    407421       ENDIF
    408422
    409423       DO  j = nys, nyn
    410 !
    411 !--       If required, calculate cloud microphysics
    412           IF ( cloud_physics  .AND.  .NOT. microphysics_sat_adjust  .AND.      &
    413                ( intermediate_timestep_count == 1  .OR.                        &
    414                  call_microphysics_at_all_substeps )                           &
    415              )  THEN
    416              CALL microphysics_control( i, j )
    417           ENDIF
    418424!
    419425!--       Tendency terms for u-velocity component
     
    424430                IF ( ws_scheme_mom )  THEN
    425431                   CALL advec_u_ws( i, j, i_omp_start, tn )
    426                 ELSE 
     432                ELSE
    427433                   CALL advec_u_pw( i, j )
    428                 ENDIF 
     434                ENDIF
    429435             ELSE
    430436                CALL advec_u_up( i, j )
     
    490496                IF ( ws_scheme_mom )  THEN
    491497                    CALL advec_v_ws( i, j, i_omp_start, tn )
    492                 ELSE 
     498                ELSE
    493499                    CALL advec_v_pw( i, j )
    494500                ENDIF
     
    501507!
    502508!--          Drag by plant canopy
    503              IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )       
     509             IF ( plant_canopy )  CALL pcm_tendency( i, j, 2 )
    504510
    505511!
     
    551557             IF ( ws_scheme_mom )  THEN
    552558                CALL advec_w_ws( i, j, i_omp_start, tn )
    553              ELSE 
     559             ELSE
    554560                CALL advec_w_pw( i, j )
    555561             END IF
     
    646652             IF ( large_scale_forcing )  THEN
    647653                CALL ls_advec( i, j, simulated_time, 'pt' )
    648              ENDIF     
     654             ENDIF
    649655
    650656!
    651657!--          Nudging
    652              IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' ) 
     658             IF ( nudging )  CALL nudge( i, j, simulated_time, 'pt' )
    653659
    654660!
     
    708714                    CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa,  &
    709715                                diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn  )
    710                 ELSE 
     716                ELSE
    711717                    CALL advec_s_pw( i, j, sa )
    712718                ENDIF
     
    760766             THEN
    761767                IF ( ws_scheme_sca )  THEN
    762                    CALL advec_s_ws( i, j, q, 'q', flux_s_q, & 
     768                   CALL advec_s_ws( i, j, q, 'q', flux_s_q, &
    763769                                diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn )
    764                 ELSE 
     770                ELSE
    765771                   CALL advec_s_pw( i, j, q )
    766772                ENDIF
     
    782788!
    783789!--          Nudging
    784              IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' ) 
     790             IF ( nudging )  CALL nudge( i, j, simulated_time, 'q' )
    785791
    786792!
     
    820826
    821827!
    822 !--          If required, calculate prognostic equations for rain water content 
     828!--          If required, calculate prognostic equations for rain water content
    823829!--          and rain drop concentration
    824830             IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    829835                THEN
    830836                   IF ( ws_scheme_sca )  THEN
    831                       CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       & 
     837                      CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr,       &
    832838                                       diss_s_qr, flux_l_qr, diss_l_qr, &
    833839                                       i_omp_start, tn )
    834                    ELSE 
     840                   ELSE
    835841                      CALL advec_s_pw( i, j, qr )
    836842                   ENDIF
     
    869875                IF ( timestep_scheme(1:5) == 'runge' )  THEN
    870876                   IF ( ws_scheme_sca )  THEN
    871                       CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    & 
     877                      CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr,    &
    872878                                    diss_s_nr, flux_l_nr, diss_l_nr, &
    873879                                    i_omp_start, tn )
    874                    ELSE 
     880                   ELSE
    875881                      CALL advec_s_pw( i, j, nr )
    876882                   ENDIF
     
    907913
    908914          ENDIF
    909          
     915
    910916!
    911917!--       If required, compute prognostic equation for scalar
     
    917923             THEN
    918924                IF ( ws_scheme_sca )  THEN
    919                    CALL advec_s_ws( i, j, s, 's', flux_s_s, & 
     925                   CALL advec_s_ws( i, j, s, 's', flux_s_s, &
    920926                                diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn )
    921                 ELSE 
     927                ELSE
    922928                   CALL advec_s_pw( i, j, s )
    923929                ENDIF
     
    939945!
    940946!--          Nudging, still need to be extended for scalars
    941 !              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' ) 
     947!              IF ( nudging )  CALL nudge( i, j, simulated_time, 's' )
    942948
    943949!
    944950!--          If required compute influence of large-scale subsidence/ascent.
    945 !--          Note, the last argument is of no meaning in this case, as it is 
    946 !--          only used in conjunction with large_scale_forcing, which is to 
     951!--          Note, the last argument is of no meaning in this case, as it is
     952!--          only used in conjunction with large_scale_forcing, which is to
    947953!--          date not implemented for scalars.
    948954             IF ( large_scale_subsidence  .AND.                                &
     
    980986             ENDIF
    981987
    982           ENDIF         
    983 !
    984 !--       If required, compute prognostic equation for turbulent kinetic 
     988          ENDIF
     989!
     990!--       If required, compute prognostic equation for turbulent kinetic
    985991!--       energy (TKE)
    986992          IF ( .NOT. constant_diffusion )  THEN
     
    990996             tend(:,j,i) = 0.0_wp
    991997             IF ( timestep_scheme(1:5) == 'runge'  &
    992                  .AND.  .NOT. use_upstream_for_tke )  THEN 
     998                 .AND.  .NOT. use_upstream_for_tke )  THEN
    993999                 IF ( ws_scheme_sca )  THEN
    9941000                     CALL advec_s_ws( i, j, e, 'e', flux_s_e, diss_s_e, &
     
    10131019!
    10141020!--          Additional sink term for flows through plant canopies
    1015              IF ( plant_canopy )  CALL pcm_tendency( i, j, 6 ) 
     1021             IF ( plant_canopy )  CALL pcm_tendency( i, j, 6 )
    10161022
    10171023             CALL user_actions( i, j, 'e-tendency' )
     
    10611067!> Version for vector machines
    10621068!------------------------------------------------------------------------------!
    1063  
     1069
    10641070 SUBROUTINE prognostic_equations_vector
    10651071
     
    11761182       IF ( ws_scheme_mom )  THEN
    11771183          CALL advec_v_ws
    1178        ELSE 
     1184       ELSE
    11791185          CALL advec_v_pw
    11801186       END IF
     
    13891395!
    13901396!--    Nudging
    1391        IF ( nudging )  CALL nudge( simulated_time, 'pt' ) 
     1397       IF ( nudging )  CALL nudge( simulated_time, 'pt' )
    13921398
    13931399!
     
    14851491
    14861492       CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux )
    1487        
     1493
    14881494       CALL user_actions( 'sa-tendency' )
    14891495
     
    15731579
    15741580       CALL diffusion_s( q, qsws, qswst, wall_qflux )
    1575        
     1581
    15761582!
    15771583!--    Sink or source of humidity due to canopy elements
     
    15861592!
    15871593!--    Nudging
    1588        IF ( nudging )  CALL nudge( simulated_time, 'q' ) 
     1594       IF ( nudging )  CALL nudge( simulated_time, 'q' )
    15891595
    15901596!
     
    16371643
    16381644!
    1639 !--    If required, calculate prognostic equations for rain water content 
     1645!--    If required, calculate prognostic equations for rain water content
    16401646!--    and rain drop concentration
    16411647       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
     
    16751681          CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux )
    16761682
    1677           CALL user_actions( 'qr-tendency' )
    1678 
    16791683!
    16801684!--       Prognostic equation for rain water content
     
    18281832
    18291833       CALL diffusion_s( s, ssws, sswst, wall_sflux )
    1830        
     1834
    18311835!
    18321836!--    Sink or source of humidity due to canopy elements
     
    18411845!
    18421846!--    Nudging. Not implemented for scalars so far.
    1843 !        IF ( nudging )  CALL nudge( simulated_time, 'q' ) 
     1847!        IF ( nudging )  CALL nudge( simulated_time, 'q' )
    18441848
    18451849!
     
    18951899    ENDIF
    18961900!
    1897 !-- If required, compute prognostic equation for turbulent kinetic 
     1901!-- If required, compute prognostic equation for turbulent kinetic
    18981902!-- energy (TKE)
    18991903    IF ( .NOT. constant_diffusion )  THEN
     
    19271931                IF ( ws_scheme_sca )  THEN
    19281932                   CALL advec_s_ws( e, 'e' )
    1929                 ELSE 
     1933                ELSE
    19301934                   CALL advec_s_pw( e )
    19311935                ENDIF
Note: See TracChangeset for help on using the changeset viewer.