Changeset 19 for palm/trunk/SOURCE/production_e.f90

Timestamp:

Feb 23, 2007 4:53:48 AM (16 years ago)

Author:

raasch

Message:

preliminary version of modified boundary conditions at top

File:

• palm/trunk/SOURCE/production_e.f90 (modified) (12 diffs)

palm/trunk/SOURCE/production_e.f90

@@ -4 +4 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation extended for gridpoint nzt, extended for given temperature /
+! humidity fluxes at the top
 !
 ! Former revisions:
@@ -67 +68 @@
 
           DO  j = nys, nyn
-             DO  k = nzb_diff_s_outer(j,i), nzt-1
+             DO  k = nzb_diff_s_outer(j,i), nzt
 
                 dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
@@ -327 +328 @@
              DO  j = nys, nyn
 
-                DO  k = nzb_diff_s_inner(j,i), nzt-1
+                DO  k = nzb_diff_s_inner(j,i), nzt_diff
                    tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
                                     ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
                 ENDDO
+
                 IF ( use_surface_fluxes )  THEN
                    k = nzb_diff_s_inner(j,i)-1
@@ -336 +338 @@
                 ENDIF
 
+                IF ( use_top_fluxes )  THEN
+                   k = nzt
+                   tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
+                ENDIF
+
              ENDDO
 
@@ -342 +349 @@
              DO  j = nys, nyn
 
-                DO  k = nzb_diff_s_inner(j,i), nzt-1
+                DO  k = nzb_diff_s_inner(j,i), nzt_diff
 
                    IF ( .NOT. cloud_physics )  THEN
@@ -375 +382 @@
                 DO  j = nys, nyn
 
-                   k = nzb_diff_s_inner(j,i)-1
+                   k = nzb_diff_s_inner(j,i)
 
                    IF ( .NOT. cloud_physics )  THEN
@@ -402 +409 @@
              ENDIF
 
+             IF ( use_top_fluxes )  THEN
+
+                DO  j = nys, nyn
+
+                   k = nzt
+
+                   IF ( .NOT. cloud_physics )  THEN
+                      k1 = 1.0 + 0.61 * q(k,j,i)
+                      k2 = 0.61 * pt(k,j,i)
+                   ELSE
+                      IF ( ql(k,j,i) == 0.0 )  THEN
+                         k1 = 1.0 + 0.61 * q(k,j,i)
+                         k2 = 0.61 * pt(k,j,i)
+                      ELSE
+                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                         temp  = theta * t_d_pt(k)
+                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+                                    ( q(k,j,i) - ql(k,j,i) ) *          &
+                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                      ENDIF
+                   ENDIF
+
+                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
+                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
+                ENDDO
+
+             ENDIF
+
           ENDIF
 
@@ -430 +468 @@
 !
 !--    Calculate TKE production by shear
-       DO  k = nzb_diff_s_outer(j,i), nzt-1
+       DO  k = nzb_diff_s_outer(j,i), nzt
 
           dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
@@ -657 +695 @@
        IF ( .NOT. moisture )  THEN
 
-          DO  k = nzb_diff_s_inner(j,i), nzt-1
+          DO  k = nzb_diff_s_inner(j,i), nzt_diff
              tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
                                     ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
           ENDDO
+
           IF ( use_surface_fluxes )  THEN
              k = nzb_diff_s_inner(j,i)-1
@@ -666 +705 @@
           ENDIF
 
+          IF ( use_top_fluxes )  THEN
+             k = nzt
+             tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
+          ENDIF
+
        ELSE
 
-          DO  k = nzb_diff_s_inner(j,i), nzt-1
+          DO  k = nzb_diff_s_inner(j,i), nzt_diff
 
              IF ( .NOT. cloud_physics )  THEN
@@ -694 +738 @@
                                       ) * dd2zu(k)
           ENDDO
+
           IF ( use_surface_fluxes )  THEN
              k = nzb_diff_s_inner(j,i)-1
@@ -720 +765 @@
           ENDIF
 
+          IF ( use_top_fluxes )  THEN
+             k = nzt
+
+             IF ( .NOT. cloud_physics ) THEN
+                k1 = 1.0 + 0.61 * q(k,j,i)
+                k2 = 0.61 * pt(k,j,i)
+             ELSE
+                IF ( ql(k,j,i) == 0.0 )  THEN
+                   k1 = 1.0 + 0.61 * q(k,j,i)
+                   k2 = 0.61 * pt(k,j,i)
+                ELSE
+                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
+                   temp  = theta * t_d_pt(k)
+                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
+                              ( q(k,j,i) - ql(k,j,i) ) *          &
+                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
+                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
+                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
+                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
+                ENDIF
+             ENDIF
+
+             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
+                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
+          ENDIF
+
        ENDIF