Changeset 1361 for palm/trunk/SOURCE/prandtl_fluxes.f90

Timestamp:

Apr 16, 2014 3:17:48 PM (10 years ago)

Author:

hoffmann

Message:

improved version of two-moment cloud physics

File:

• palm/trunk/SOURCE/prandtl_fluxes.f90 (modified) (6 diffs)

palm/trunk/SOURCE/prandtl_fluxes.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain 
+! drop concentration (nrsws) added
 ! 
 ! Former revisions:
@@ -64 +65 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  e, pt, q, qs, qsws, rif, shf, ts, u, us, usws, v, vpt, vsws,    &
-               zu, zw, z0, z0h
+        ONLY:  e, nr, nrs, nrsws, pt, q, qr, qrs, qrsws, qs, qsws, rif, shf,   &
+               ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h
 
     USE control_parameters,                                                    &
-        ONLY:  constant_heatflux, constant_waterflux, coupling_mode, g,        &
-               humidity, ibc_e_b, kappa, large_scale_forcing, lsf_surf,        &
-               passive_scalar, pt_surface, q_surface, rif_max, rif_min,        &
-               run_coupled, surface_pressure
+        ONLY:  cloud_physics, constant_heatflux, constant_waterflux,           &
+               coupling_mode, g, humidity, ibc_e_b, icloud_scheme, kappa,      &
+               large_scale_forcing, lsf_surf, passive_scalar, precipitation,   &
+               pt_surface, q_surface, rif_max, rif_min, run_coupled,           &
+               surface_pressure
 
     USE indices,                                                               &
@@ -96 +98 @@
 !
 !-- Data information for accelerators
-    !$acc data present( e, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt, q, qs ) &
-    !$acc      present( qsws, rif, shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h )
+    !$acc data present( e, nrsws, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt )  &
+    !$acc      present( q, qs, qsws, qrsws, rif, shf, ts, u, us, usws, v )     &
+    !$acc      present( vpt, vsws, zu, zw, z0, z0h )
 ! 
 !-- Compute theta*
@@ -383 +386 @@
           ENDDO
        ENDIF
+
+       IF ( cloud_physics  .AND.  icloud_scheme == 0  &
+            .AND.  precipitation )  THEN
+
+          !$OMP PARALLEL DO PRIVATE( a, b, k, z_p )
+          !$acc kernels loop independent
+          DO  i = nxlg, nxrg
+             !$acc loop independent
+             DO  j = nysg, nyng
+
+                k   = nzb_s_inner(j,i)
+                z_p = zu(k+1) - zw(k)
+
+                IF ( rif(j,i) >= 0.0 )  THEN
+!
+!--                Stable stratification
+                   qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) / (          &
+                                  LOG( z_p / z0h(j,i) ) +                      &
+                                  5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p )
+                   nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) / (          &
+                                  LOG( z_p / z0h(j,i) ) +                      &
+                                  5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p )
+
+                ELSE
+!
+!--                Unstable stratification
+                   a = SQRT( 1.0 - 16.0 * rif(j,i) ) 
+                   b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p ) 
+ 
+                   qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) /   (        &
+                             LOG( z_p / z0h(j,i) ) -                           &
+                             2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
+                   nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) /   (        &
+                             LOG( z_p / z0h(j,i) ) -                           &
+                             2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
+
+                ENDIF
+
+             ENDDO
+          ENDDO
+
+       ENDIF
+
     ENDIF
 
@@ -396 +442 @@
        CALL exchange_horiz_2d( qsws )
        !$acc update device( qsws )
+       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  &
+            precipitation )  THEN
+          !$acc update host( qrsws, nrsws )
+          CALL exchange_horiz_2d( qrsws )
+          CALL exchange_horiz_2d( nrsws )
+          !$acc update device( qrsws, nrsws )
+       ENDIF
     ENDIF
 
@@ -425 +478 @@
 
 !
+!-- Compute (turbulent) fluxes of rain water content and rain drop concentartion
+    IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
+       !$OMP PARALLEL DO
+       !$acc kernels loop independent
+       DO  i = nxlg, nxrg
+          !$acc loop independent
+          DO  j = nysg, nyng
+             qrsws(j,i) = -qrs(j,i) * us(j,i)
+             nrsws(j,i) = -nrs(j,i) * us(j,i)
+          ENDDO
+       ENDDO
+    ENDIF
+
+!
 !-- Bottom boundary condition for the TKE
     IF ( ibc_e_b == 2 )  THEN