Changeset 65 for palm/trunk/SOURCE/diffusion_e.f90
 Timestamp:
 Mar 13, 2007 12:11:43 PM (17 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/diffusion_e.f90
r57 r65 64 64 65 65 66 DO i = nxl, nxr 67 DO j = nys, nyn 68 ! 69 ! First, calculate phifunction for eventually adjusting the & 70 ! mixing length to the prandtl mixing length 71 IF ( adjust_mixing_length .AND. prandtl_layer ) THEN 72 IF ( rif(j,i) >= 0.0 ) THEN 73 phi_m = 1.0 + 5.0 * rif(j,i) 74 ELSE 75 phi_m = 1.0 / SQRT( SQRT( 1.0  16.0 * rif(j,i) ) ) 66 ! 67 ! This if clause must be outside the kloop because otherwise 68 ! runtime errors occur with C hopt on NEC 69 IF ( use_pt_reference ) THEN 70 71 DO i = nxl, nxr 72 DO j = nys, nyn 73 ! 74 ! First, calculate phifunction for eventually adjusting the & 75 ! mixing length to the prandtl mixing length 76 IF ( adjust_mixing_length .AND. prandtl_layer ) THEN 77 IF ( rif(j,i) >= 0.0 ) THEN 78 phi_m = 1.0 + 5.0 * rif(j,i) 79 ELSE 80 phi_m = 1.0 / SQRT( SQRT( 1.0  16.0 * rif(j,i) ) ) 81 ENDIF 76 82 ENDIF 77 ENDIF 78 79 DO k = nzb_s_inner(j,i)+1, nzt 80 ! 81 ! Calculate the mixing length (for dissipation) 82 dpt_dz = ( theta(k+1,j,i)  theta(k1,j,i) ) * dd2zu(k) 83 IF ( dpt_dz > 0.0 ) THEN 84 IF ( use_pt_reference ) THEN 83 84 DO k = nzb_s_inner(j,i)+1, nzt 85 ! 86 ! Calculate the mixing length (for dissipation) 87 dpt_dz = ( theta(k+1,j,i)  theta(k1,j,i) ) * dd2zu(k) 88 IF ( dpt_dz > 0.0 ) THEN 85 89 l_stable = 0.76 * SQRT( e(k,j,i) ) / & 86 90 SQRT( g / pt_reference * dpt_dz ) + 1E5 87 91 ELSE 88 l_stable = 0.76 * SQRT( e(k,j,i) ) / & 89 SQRT( g / theta(k,j,i) * dpt_dz ) + 1E5 90 ENDIF 91 ELSE 92 l_stable = l_grid(k) 93 ENDIF 94 ! 95 ! Adjustment of the mixing length 96 IF ( wall_adjustment ) THEN 97 l(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k), & 98 l_stable ) 99 ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) ) 100 ELSE 101 l(k,j) = MIN( l_grid(k), l_stable ) 102 ll(k,j) = l_grid(k) 103 ENDIF 104 IF ( adjust_mixing_length .AND. prandtl_layer ) THEN 105 l(k,j) = MIN( l(k,j), kappa * zu(k) / phi_m ) 106 ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m ) 107 ENDIF 108 109 ENDDO 110 ENDDO 111 ! 112 ! Calculate the tendency terms 113 DO j = nys, nyn 114 DO k = nzb_s_inner(j,i)+1, nzt 115 116 dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * & 117 e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j) 118 119 tend(k,j,i) = tend(k,j,i) & 92 l_stable = l_grid(k) 93 ENDIF 94 ! 95 ! Adjustment of the mixing length 96 IF ( wall_adjustment ) THEN 97 l(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k),& 98 l_stable ) 99 ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) ) 100 ELSE 101 l(k,j) = MIN( l_grid(k), l_stable ) 102 ll(k,j) = l_grid(k) 103 ENDIF 104 IF ( adjust_mixing_length .AND. prandtl_layer ) THEN 105 l(k,j) = MIN( l(k,j), kappa * zu(k) / phi_m ) 106 ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m ) 107 ENDIF 108 109 ENDDO 110 ENDDO 111 112 ! 113 ! Calculate the tendency terms 114 DO j = nys, nyn 115 DO k = nzb_s_inner(j,i)+1, nzt 116 117 dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * & 118 e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j) 119 120 tend(k,j,i) = tend(k,j,i) & 120 121 + ( & 121 122 ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)e(k,j,i) ) & … … 132 133  dissipation(k,j) 133 134 134 ENDDO 135 ENDDO 136 ENDDO 137 138 ! 139 ! Store dissipation if needed for calculating the sgs particle 140 ! velocities 141 IF ( use_sgs_for_particles ) THEN 142 DO j = nys, nyn 143 DO k = nzb_s_inner(j,i)+1, nzt 144 diss(k,j,i) = dissipation(k,j) 145 ENDDO 146 ENDDO 147 ENDIF 148 135 149 ENDDO 136 150 137 ! 138 ! Store dissipation if needed for calculating the sgs particle 139 ! velocities 140 IF ( use_sgs_for_particles ) THEN 141 DO j = nys, nyn 151 ELSE 152 153 DO i = nxl, nxr 154 DO j = nys, nyn 155 ! 156 ! First, calculate phifunction for eventually adjusting the & 157 ! mixing length to the prandtl mixing length 158 IF ( adjust_mixing_length .AND. prandtl_layer ) THEN 159 IF ( rif(j,i) >= 0.0 ) THEN 160 phi_m = 1.0 + 5.0 * rif(j,i) 161 ELSE 162 phi_m = 1.0 / SQRT( SQRT( 1.0  16.0 * rif(j,i) ) ) 163 ENDIF 164 ENDIF 165 142 166 DO k = nzb_s_inner(j,i)+1, nzt 143 diss(k,j,i) = dissipation(k,j) 144 ENDDO 145 ENDDO 146 ENDIF 147 148 ENDDO 167 ! 168 ! Calculate the mixing length (for dissipation) 169 dpt_dz = ( theta(k+1,j,i)  theta(k1,j,i) ) * dd2zu(k) 170 IF ( dpt_dz > 0.0 ) THEN 171 l_stable = 0.76 * SQRT( e(k,j,i) ) / & 172 SQRT( g / theta(k,j,i) * dpt_dz ) + 1E5 173 ELSE 174 l_stable = l_grid(k) 175 ENDIF 176 ! 177 ! Adjustment of the mixing length 178 IF ( wall_adjustment ) THEN 179 l(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k),& 180 l_stable ) 181 ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) ) 182 ELSE 183 l(k,j) = MIN( l_grid(k), l_stable ) 184 ll(k,j) = l_grid(k) 185 ENDIF 186 IF ( adjust_mixing_length .AND. prandtl_layer ) THEN 187 l(k,j) = MIN( l(k,j), kappa * zu(k) / phi_m ) 188 ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m ) 189 ENDIF 190 191 ENDDO 192 ENDDO 193 194 ! 195 ! Calculate the tendency terms 196 DO j = nys, nyn 197 DO k = nzb_s_inner(j,i)+1, nzt 198 199 dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * & 200 e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j) 201 202 tend(k,j,i) = tend(k,j,i) & 203 + ( & 204 ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)e(k,j,i) ) & 205  ( km(k,j,i)+km(k,j,i1) ) * ( e(k,j,i)e(k,j,i1) ) & 206 ) * ddx2 & 207 + ( & 208 ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)e(k,j,i) ) & 209  ( km(k,j,i)+km(k,j1,i) ) * ( e(k,j,i)e(k,j1,i) ) & 210 ) * ddy2 & 211 + ( & 212 ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)e(k,j,i) ) * ddzu(k+1) & 213  ( km(k,j,i)+km(k1,j,i) ) * ( e(k,j,i)e(k1,j,i) ) * ddzu(k) & 214 ) * ddzw(k) & 215  dissipation(k,j) 216 217 ENDDO 218 ENDDO 219 220 ! 221 ! Store dissipation if needed for calculating the sgs particle 222 ! velocities 223 IF ( use_sgs_for_particles ) THEN 224 DO j = nys, nyn 225 DO k = nzb_s_inner(j,i)+1, nzt 226 diss(k,j,i) = dissipation(k,j) 227 ENDDO 228 ENDDO 229 ENDIF 230 231 ENDDO 232 233 ENDIF 149 234 150 235 !
Note: See TracChangeset
for help on using the changeset viewer.