Index: palm/trunk/SOURCE/prandtl_fluxes.f90
===================================================================
--- palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 1362)
+++ palm/trunk/SOURCE/prandtl_fluxes.f90	(revision 1494)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfixes: qs is now calculated before calculation of Rif. Ccalculation of 
+! buoyancy flux in Rif corrected (added missing humidity term), allow use of 
+! topography for coupled runs (not tested)
 ! 
 ! Former revisions:
@@ -129,5 +131,5 @@
    
        IF ( large_scale_forcing .AND. lsf_surf )  THEN
-          pt(0,:,:) = pt_surface
+          pt(nzb_s_inner(j,i),:,:) = pt_surface
        ENDIF
 
@@ -153,6 +155,6 @@
                 b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p )
 
-                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (          &
-                          LOG( z_p / z0h(j,i) ) -                           &
+                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (             &
+                          LOG( z_p / z0h(j,i) ) -                              &
                           2.0_wp * LOG( ( 1.0_wp + a ) / ( 1.0_wp + b ) ) )
              ENDIF
@@ -160,4 +162,65 @@
           ENDDO
        ENDDO
+    ENDIF
+
+!
+!-- If required compute q*
+    IF ( humidity  .OR.  passive_scalar )  THEN
+       IF ( constant_waterflux )  THEN
+!
+!--       For a given water flux in the Prandtl layer:
+          !$OMP PARALLEL DO
+          !$acc kernels loop
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
+             ENDDO
+          ENDDO
+          
+       ELSE
+          coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. run_coupled )
+
+           IF ( large_scale_forcing .AND. lsf_surf )  THEN
+              q(nzb_s_inner(j,i),:,:) = q_surface
+           ENDIF
+
+          !$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)
+
+!
+!--             Assume saturation for atmosphere coupled to ocean (but not
+!--             in case of precursor runs)
+                IF ( coupled_run )  THEN
+                   e_q = 6.1_wp * &
+                           EXP( 0.07_wp * ( MIN(pt(k,j,i),pt(k+1,j,i)) - 273.15_wp ) )
+                   q(k,j,i) = 0.622_wp * e_q / ( surface_pressure - e_q )
+                ENDIF
+                IF ( rif(j,i) >= 0.0_wp )  THEN
+!
+!--                Stable stratification
+                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (             &
+                                  LOG( z_p / z0h(j,i) ) +                      &
+                                  5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
+                                                                 )
+                ELSE
+!
+!--                Unstable stratification
+                   a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) 
+                   b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p ) 
+ 
+                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (           &
+                             LOG( z_p / z0h(j,i) ) -                           &
+                              2.0_wp * LOG( (1.0_wp + a ) / ( 1.0_wp + b ) ) )
+                ENDIF
+
+             ENDDO
+          ENDDO
+       ENDIF
     ENDIF
 
@@ -190,6 +253,7 @@
              k   = nzb_s_inner(j,i)
              z_p = zu(k+1) - zw(k)
-             rif(j,i) = z_p * kappa * g *                            &
-                        ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i) ) / &
+             rif(j,i) = z_p * kappa * g *                                      &
+                        ( ts(j,i) + 0.61_wp * pt(k+1,j,i) * qs(j,i) + 0.61_wp  &
+                          * q(k+1,j,i) * ts(j,i)) /                            &
                         ( vpt(k+1,j,i) * ( us(j,i)**2 + 1E-30_wp ) )
 !
@@ -218,7 +282,7 @@
 !--       Compute the absolute value of the horizontal velocity 
 !--       (relative to the surface)
-          uv_total = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)        &
-                                      - u(k,j,i)   - u(k,j,i+1) ) )**2 + &
-                           ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)        &
+          uv_total = SQRT( ( 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)              &
+                                      - u(k,j,i)   - u(k,j,i+1) ) )**2 +       &
+                           ( 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)              &
                                       - v(k,j,i)   - v(k,j+1,i) ) )**2 )    
 
@@ -227,7 +291,7 @@
 !
 !--          Stable stratification
-             us(j,i) = kappa * uv_total / (                                &
-                                  LOG( z_p / z0(j,i) ) +                   &
-                                  5.0_wp * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
+             us(j,i) = kappa * uv_total / (                                    &
+                                  LOG( z_p / z0(j,i) ) +                       &
+                                  5.0_wp * rif(j,i) * ( z_p - z0(j,i) ) / z_p  &
                                           )
           ELSE
@@ -237,9 +301,9 @@
              b = SQRT( SQRT( 1.0_wp - 16.0_wp * rif(j,i) / z_p * z0(j,i) ) )
 
-             us(j,i) = kappa * uv_total / (                                  &
-                       LOG( z_p / z0(j,i) ) -                                &
-                       LOG( ( 1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (              &
-                            ( 1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +          &
-                               2.0_wp * ( ATAN( a ) - ATAN( b ) )                  &
+             us(j,i) = kappa * uv_total / (                                    &
+                       LOG( z_p / z0(j,i) ) -                                  &
+                       LOG( ( 1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (          &
+                            ( 1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +      &
+                               2.0_wp * ( ATAN( a ) - ATAN( b ) )              &
                                            )
           ENDIF
@@ -271,7 +335,7 @@
 !
 !--          Stable stratification
-             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )/ (              &
-                                     LOG( z_p / z0(j,i) ) +               &
-                                     5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p &
+             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) )/ (                  &
+                                     LOG( z_p / z0(j,i) ) +                    &
+                                     5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p   &
                                                             )
           ELSE
@@ -281,9 +345,9 @@
              b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) )
 
-             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) ) / (            &
-                         LOG( z_p / z0(j,i) ) -                           &
-                         LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (          &
-                              (1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +      &
-                                 2.0_wp * ( ATAN( a ) - ATAN( b ) )             &
+             usws(j,i) = kappa * ( u(k+1,j,i) - u(k,j,i) ) / (                 &
+                         LOG( z_p / z0(j,i) ) -                                &
+                         LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (         &
+                              (1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +     &
+                                 2.0_wp * ( ATAN( a ) - ATAN( b ) )            &
                                                  )
           ENDIF
@@ -309,7 +373,7 @@
 !
 !--          Stable stratification
-             vsws(j,i) = kappa * ( v(k+1,j,i) -  v(k,j,i) ) / (           &
-                                     LOG( z_p / z0(j,i) ) +               &
-                                     5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p &
+             vsws(j,i) = kappa * ( v(k+1,j,i) -  v(k,j,i) ) / (                &
+                                     LOG( z_p / z0(j,i) ) +                    &
+                                     5.0_wp * rifm * ( z_p - z0(j,i) ) / z_p   &
                                                               )
           ELSE
@@ -319,9 +383,9 @@
              b = SQRT( SQRT( 1.0_wp - 16.0_wp * rifm / z_p * z0(j,i) ) )
 
-             vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / (            &
-                         LOG( z_p / z0(j,i) ) -                           &
-                         LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (          &
-                              (1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +      &
-                                 2.0_wp * ( ATAN( a ) - ATAN( b ) )             &
+             vsws(j,i) = kappa * ( v(k+1,j,i) - v(k,j,i) ) / (                 &
+                         LOG( z_p / z0(j,i) ) -                                &
+                         LOG( (1.0_wp + a )**2 * ( 1.0_wp + a**2 ) / (         &
+                              (1.0_wp + b )**2 * ( 1.0_wp + b**2 )   ) ) +     &
+                                 2.0_wp * ( ATAN( a ) - ATAN( b ) )            &
                                                  )
           ENDIF
@@ -331,104 +395,43 @@
 
 !
-!-- If required compute q*
-    IF ( humidity  .OR.  passive_scalar )  THEN
-       IF ( constant_waterflux )  THEN
-!
-!--       For a given water flux in the Prandtl layer:
-          !$OMP PARALLEL DO
-          !$acc kernels loop
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-                qs(j,i) = -qsws(j,i) / ( us(j,i) + 1E-30_wp )
-             ENDDO
-          ENDDO
-          
-       ELSE
-          coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND. run_coupled )
-
-           IF ( large_scale_forcing .AND. lsf_surf )  THEN
-              q(0,:,:) = q_surface
-           ENDIF
-
-          !$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)
-
-!
-!--             Assume saturation for atmosphere coupled to ocean (but not
-!--             in case of precursor runs)
-                IF ( coupled_run )  THEN
-                   e_q = 6.1_wp * &
-                           EXP( 0.07_wp * ( MIN(pt(0,j,i),pt(1,j,i)) - 273.15_wp ) )
-                   q(k,j,i) = 0.622_wp * e_q / ( surface_pressure - e_q )
-                ENDIF
-                IF ( rif(j,i) >= 0.0_wp )  THEN
-!
-!--                Stable stratification
-                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (          &
-                                  LOG( z_p / z0h(j,i) ) +                   &
-                                  5.0_wp * rif(j,i) * ( z_p - z0h(j,i) ) / z_p &
-                                                                 )
-                ELSE
-!
-!--                Unstable stratification
-                   a = SQRT( 1.0_wp - 16.0_wp * rif(j,i) ) 
-                   b = SQRT( 1.0_wp - 16.0_wp * rif(j,i) * z0h(j,i) / z_p ) 
+!-- If required compute qr* and nr*
+    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 ) 
  
-                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (        &
-                             LOG( z_p / z0h(j,i) ) -                        &
-                              2.0_wp * LOG( (1.0_wp + a ) / ( 1.0_wp + b ) ) )
-                ENDIF
-
-             ENDDO
-          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
+                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
