Index: palm/trunk/SOURCE/surface_layer_fluxes.f90
===================================================================
--- palm/trunk/SOURCE/surface_layer_fluxes.f90	(revision 1758)
+++ palm/trunk/SOURCE/surface_layer_fluxes.f90	(revision 1788)
@@ -19,5 +19,7 @@
 ! Current revisions:
 ! ------------------
-! 
+! Added parameter z0q which replaces z0h in the similarity functions for
+! humidity.
+! Syntax layout improved.
 ! 
 ! Former revisions:
@@ -117,5 +119,5 @@
     USE arrays_3d,                                                             &
         ONLY:  e, kh, nr, nrs, nrsws, ol, pt, q, ql, qr, qrs, qrsws, qs, qsws, &
-               shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h
+               shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h, z0q
 
     USE cloud_parameters,                                                      &
@@ -234,5 +236,5 @@
 !--    Use either Newton iteration or a lookup table for the bulk Richardson
 !--    number to calculate the Obukhov length 
-       ELSEIF ( most_method == 'newton' .OR. most_method == 'lookup' )  THEN
+       ELSEIF ( most_method == 'newton'  .OR.  most_method == 'lookup' )  THEN
 
           CALL calc_uv_total
@@ -302,5 +304,5 @@
 !--       Check for roughness heterogeneity. In that case terminate run and
 !--       inform user
-          IF ( MINVAL( z0h ) /= MAXVAL( z0h ) .OR.                             &
+          IF ( MINVAL( z0h ) /= MAXVAL( z0h )  .OR.                            &
                MINVAL( z0  ) /= MAXVAL( z0  ) )  THEN
              terminate_run_l = .TRUE.
@@ -347,5 +349,5 @@
 !
 !--       Calculate stretching factor for the free convection range
-          DO WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
+          DO  WHILE ( ABS( (regr-regr_old) / regr_old ) > 1.0E-10_wp )
              regr_old = regr
              regr = ( 1.0_wp - ( -zeta_min / zeta_step ) * ( 1.0_wp - regr )   &
@@ -533,5 +535,5 @@
 !--             Ensure that the bulk Richardson number and the Obukhov 
 !--             lengtH have the same sign
-                IF ( rib(j,i) * ol(j,i) < 0.0_wp .OR.                        &
+                IF ( rib(j,i) * ol(j,i) < 0.0_wp  .OR.                         &
                      ABS( ol(j,i) ) == ol_max )  THEN
                    IF ( rib(j,i) > 0.0_wp ) ol(j,i) =  0.01_wp
@@ -638,5 +640,5 @@
 
           !$OMP PARALLEL DO PRIVATE( k, z_mo )
-          !# WARNING: does not work on GPU so far because of DO WHILE construct
+          !# WARNING: does not work on GPU so far because of DO  WHILE construct
           !!!!!!$acc kernels loop
           DO  i = nxl, nxr
@@ -656,10 +658,10 @@
                 l = l_bnd
                 IF ( rib_tab(l) - rib(j,i) > 0.0_wp )  THEN
-                   DO WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp .AND. l > 0 )
+                   DO  WHILE ( rib_tab(l-1) - rib(j,i) > 0.0_wp  .AND.  l > 0 )
                       l = l-1
                    ENDDO
                 ELSE
-                   DO WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp                &
-                              .AND. l < num_steps-1 )
+                   DO  WHILE ( rib_tab(l) - rib(j,i) < 0.0_wp                  &
+                              .AND.  l < num_steps-1 )
                       l = l+1
                    ENDDO
@@ -669,6 +671,6 @@
 !
 !--             Linear interpolation to find the correct value of z/L
-                ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) )       &
-                            / (  rib_tab(l) - rib_tab(l-1) )                &
+                ol(j,i) = ( ol_tab(l-1) + ( ol_tab(l) - ol_tab(l-1) )          &
+                            / (  rib_tab(l) - rib_tab(l-1) )                   &
                             * ( rib(j,i) - rib_tab(l-1) ) )
 
@@ -797,5 +799,5 @@
 !
 !--       For a given surface temperature:
-          IF ( large_scale_forcing .AND. lsf_surf )  THEN
+          IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
              !$OMP PARALLEL DO
              !$acc kernels loop private( j, k )
@@ -837,5 +839,5 @@
           IF ( constant_waterflux )  THEN
 !
-!--          For a given water flux in the Prandtl layer:
+!--          For a given water flux in the surface layer
              !$OMP PARALLEL DO
              !$acc kernels loop private( j )
@@ -847,8 +849,8 @@
 
           ELSE
-             coupled_run = ( coupling_mode == 'atmosphere_to_ocean' .AND.      &
+             coupled_run = ( coupling_mode == 'atmosphere_to_ocean'  .AND.     &
                              run_coupled )
 
-             IF ( large_scale_forcing .AND. lsf_surf )  THEN
+             IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
                 !$OMP PARALLEL DO
                 !$acc kernels loop private( j, k )
@@ -862,5 +864,5 @@
 
              !$OMP PARALLEL DO PRIVATE( e_s, k, z_mo )
-             !$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0h ) private( e_s, j, k, z_mo )
+             !$acc kernels loop independent present( nzb_s_inner, ol, pt, q, qs, qv1, zu, zw, z0q ) private( e_s, j, k, z_mo )
              DO  i = nxlg, nxrg
                 !$acc loop independent
@@ -882,13 +884,13 @@
                    IF ( cloud_physics )  THEN
                       qs(j,i) = kappa * ( qv1(j,i) - q(k,j,i) )                &
-                                        / ( LOG( z_mo / z0h(j,i) )             &
+                                        / ( LOG( z_mo / z0q(j,i) )             &
                                             - psi_h( z_mo / ol(j,i) )          &
-                                            + psi_h( z0h(j,i) / ol(j,i) ) )
+                                            + psi_h( z0q(j,i) / ol(j,i) ) )
 
                    ELSE
                       qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) )              &
-                                        / ( LOG( z_mo / z0h(j,i) )             &
+                                        / ( LOG( z_mo / z0q(j,i) )             &
                                             - psi_h( z_mo / ol(j,i) )          &
-                                            + psi_h( z0h(j,i) / ol(j,i) ) )
+                                            + psi_h( z0q(j,i) / ol(j,i) ) )
                    ENDIF
 
@@ -901,8 +903,9 @@
 !
 !--    If required compute qr* and nr*
-       IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation )  THEN
+       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )   &
+       THEN
 
           !$OMP PARALLEL DO PRIVATE( k, z_mo )
-          !$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0h ) private( j, k, z_mo )
+          !$acc kernels loop independent present( nr, nrs, nzb_s_inner, ol, qr, qrs, zu, zw, z0q ) private( j, k, z_mo )
           DO  i = nxlg, nxrg
              !$acc loop independent
@@ -913,12 +916,12 @@
 
                 qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) )              &
-                                 / ( LOG( z_mo / z0h(j,i) )                 &
+                                 / ( LOG( z_mo / z0q(j,i) )                 &
                                      - psi_h( z_mo / ol(j,i) )              &
-                                     + psi_h( z0h(j,i) / ol(j,i) ) )
+                                     + psi_h( z0q(j,i) / ol(j,i) ) )
 
                 nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) )              &
-                                 / ( LOG( z_mo / z0h(j,i) )                 &
+                                 / ( LOG( z_mo / z0q(j,i) )                 &
                                      - psi_h( z_mo / ol(j,i) )              &
-                                     + psi_h( z0h(j,i) / ol(j,i) ) )
+                                     + psi_h( z0q(j,i) / ol(j,i) ) )
              ENDDO
           ENDDO
@@ -1003,6 +1006,6 @@
 !
 !--    Compute the vertical kinematic heat flux
-       IF ( .NOT. constant_heatflux .AND. ( simulated_time <=                &
-            skip_time_do_lsm .OR. .NOT. land_surface ) )  THEN
+       IF (  .NOT.  constant_heatflux .AND.  ( simulated_time <=               &
+            skip_time_do_lsm  .OR.  .NOT.  land_surface ) )  THEN
           !$OMP PARALLEL DO
           !$acc kernels loop independent present( shf, ts, us )
@@ -1018,7 +1021,7 @@
 !
 !-- Compute the vertical water/scalar flux
-       IF ( .NOT. constant_waterflux .AND. ( humidity .OR. passive_scalar )    &
-            .AND. ( simulated_time <= skip_time_do_lsm .OR. .NOT.            &
-            land_surface ) )  THEN
+       IF (  .NOT.  constant_waterflux  .AND.  ( humidity  .OR.                &
+             passive_scalar )  .AND.  ( simulated_time <= skip_time_do_lsm     &
+            .OR.  .NOT.  land_surface ) )  THEN
           !$OMP PARALLEL DO
           !$acc kernels loop independent present( qs, qsws, us )
@@ -1089,5 +1092,5 @@
 
        IF ( zeta < 0.0_wp )  THEN
-          x = SQRT( SQRT(1.0_wp  - 16.0_wp * zeta ) )
+          x = SQRT( SQRT( 1.0_wp  - 16.0_wp * zeta ) )
           psi_m = pi * 0.5_wp - 2.0_wp * ATAN( x ) + LOG( ( 1.0_wp + x )**2    &
                   * ( 1.0_wp + x**2 ) * 0.125_wp )
@@ -1126,5 +1129,5 @@
 
        IF ( zeta < 0.0_wp )  THEN
-          x = SQRT(1.0_wp  - 16.0_wp * zeta )
+          x = SQRT( 1.0_wp  - 16.0_wp * zeta )
           psi_h = 2.0_wp * LOG( (1.0_wp + x ) / 2.0_wp )
        ELSE
