[2089] | 1 | subroutine PHY_vdfTKE_RUN |
---|
| 2 | |
---|
| 3 | !------------------------------------------------------------------------------+ |
---|
| 4 | ! Mon 17-Jun-2013 MAR | |
---|
| 5 | ! SubRoutine PHY_vdfTKE_RUN includes Vertical Diffusion of TKE | |
---|
| 6 | ! and epsilon | |
---|
| 7 | ! | |
---|
| 8 | ! | |
---|
| 9 | ! version 3.p.4.1 created by H. Gallee, Tue 19-Mar-2013 | |
---|
| 10 | ! Last Modification by H. Gallee, Mon 17-Jun-2013 | |
---|
| 11 | ! | |
---|
| 12 | !------------------------------------------------------------------------------+ |
---|
| 13 | ! | |
---|
| 14 | ! INPUT: Kzm_AT Vertical Turbulent Coeffic.(momentum) [m2/s2] | |
---|
| 15 | ! ^^^^^^ | |
---|
| 16 | ! | |
---|
| 17 | ! INPUT / OUTPUT: The Vertical Turbulent Fluxes are included for: | |
---|
| 18 | ! ^^^^^^^^^^^^^^ | |
---|
| 19 | ! a) Turbulent Kinetic Energy TKE_AT(_xyz) [m2/s2] | |
---|
| 20 | ! b) Turbulent Kinetic Energy Dissipation eps_AT(_xyz) [m2/s3] | |
---|
| 21 | ! | |
---|
| 22 | ! #OPTIONS: #De: Dirichlet Type Top Boundary Condit. for TKE_AT (TKE ) | |
---|
| 23 | ! #^^^^^^^^ & eps_AT (epsilon) | |
---|
| 24 | ! | |
---|
| 25 | !------------------------------------------------------------------------------+ |
---|
| 26 | |
---|
| 27 | use Mod_Real |
---|
| 28 | use Mod_PHY____dat |
---|
| 29 | use Mod_PHY____grd |
---|
| 30 | use Mod_PHY_AT_grd |
---|
| 31 | use Mod_PHY_AT_kkl |
---|
| 32 | use Mod_PHY_DY_kkl |
---|
| 33 | |
---|
| 34 | |
---|
| 35 | |
---|
| 36 | ! Local Variables |
---|
| 37 | ! ================ |
---|
| 38 | |
---|
| 39 | use Mod_vdfTKE_RUN |
---|
| 40 | |
---|
| 41 | |
---|
| 42 | IMPLICIT NONE |
---|
| 43 | |
---|
| 44 | |
---|
| 45 | real(kind=real8) :: S3DSBC |
---|
| 46 | real(kind=real8) :: sige2k |
---|
| 47 | real(kind=real8) :: psa_sq |
---|
| 48 | ! #De real(kind=real8) :: TKEtop = 0. |
---|
| 49 | |
---|
| 50 | integer :: i ,j ,k |
---|
| 51 | integer :: k1 ,ikl |
---|
| 52 | |
---|
| 53 | |
---|
| 54 | |
---|
| 55 | |
---|
| 56 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 57 | ! ! |
---|
| 58 | ! ALLOCATION |
---|
| 59 | ! ========== |
---|
| 60 | |
---|
| 61 | IF (it_RUN.EQ.1 .OR. FlagDALLOC) THEN ! |
---|
| 62 | allocate ( S3D__1(mzp-1) ) |
---|
| 63 | allocate ( S3D__2(mzp-1) ) |
---|
| 64 | allocate ( S3D__3(mzp-1) ) |
---|
| 65 | allocate ( S3D__4(mzp-1) ) |
---|
| 66 | allocate ( S3D__5(mzp-1) ) |
---|
| 67 | allocate ( S3D__6(mzp-1) ) |
---|
| 68 | allocate ( S3D__7(mzp-1) ) |
---|
| 69 | END IF ! |
---|
| 70 | ! ! |
---|
| 71 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 72 | |
---|
| 73 | |
---|
| 74 | |
---|
| 75 | |
---|
| 76 | ! INITIALIZATION |
---|
| 77 | ! ============== |
---|
| 78 | |
---|
| 79 | sige2k = sige / sigk |
---|
| 80 | |
---|
| 81 | |
---|
| 82 | |
---|
| 83 | ! ===================== |
---|
| 84 | DO ikl = 1,kcolp |
---|
| 85 | ! ===================== |
---|
| 86 | |
---|
| 87 | i = ii__AP(ikl) |
---|
| 88 | j = jj__AP(ikl) |
---|
| 89 | |
---|
| 90 | psa_sq = psa_DY(ikl) *psa_DY(ikl) |
---|
| 91 | |
---|
| 92 | |
---|
| 93 | |
---|
| 94 | |
---|
| 95 | ! Vertical Diffusion of Turbulent Kinetic Energy |
---|
| 96 | ! ============================================== |
---|
| 97 | |
---|
| 98 | |
---|
| 99 | ! Tridiagonal Matrix Coefficients - TKE_AT |
---|
| 100 | ! ---------------------------------------- |
---|
| 101 | |
---|
| 102 | k=mzp-1 |
---|
| 103 | S3DSBC = -GravF2*0.5*(Kzm_AT(ikl,k)+Kzm_AT(ikl,k+1)) &! SBC: TKE and epsilon, Atm SBL |
---|
| 104 | & * sigk *betaAT * roamDY(ikl,k)*roa_DY(ikl,k ) &! |
---|
| 105 | & /(psa_sq * dsigmi(k) *dsigma( k+1)) |
---|
| 106 | S3D__1(mzp-1) = S3DSBC ! SBC: TKE and epsilon, Atm SBL |
---|
| 107 | |
---|
| 108 | DO k=mzp-2,1,-1 |
---|
| 109 | S3D__1(k) =-GravF2*0.5*(Kzm_AT(ikl,k)+Kzm_AT(ikl,k+1)) & |
---|
| 110 | & * sigk *betaAT * roamDY(ikl,k)*roa_DY(ikl,k ) & |
---|
| 111 | & /(psa_sq * dsigmi(k) *dsigma( k+1)) |
---|
| 112 | |
---|
| 113 | S3D__3(k+1) = S3D__1(k) * dsigmi( k)/dsigmi( k+1) & |
---|
| 114 | & / roa_DY(ikl,k)*roa_DY(ikl,k+1) |
---|
| 115 | END DO |
---|
| 116 | |
---|
| 117 | S3D__3(1) = 0.0 ! UBC: Von Neuman , Atm Top |
---|
| 118 | DO k= 1,mzp-1 |
---|
| 119 | S3D__1(k) = S3D__1(k) * dt__AT |
---|
| 120 | S3D__3(k) = S3D__3(k) * dt__AT |
---|
| 121 | S3D__2(k) = 1.0 -S3D__3(k) -S3D__1(k) |
---|
| 122 | END DO |
---|
| 123 | |
---|
| 124 | |
---|
| 125 | ! Second Member of the Tridiagonal System - TKE_AT |
---|
| 126 | ! ------------------------------------------------ |
---|
| 127 | |
---|
| 128 | S3D__4(1) = & |
---|
| 129 | & S3D__1(1) *a_b_AT*(TKE_AT(ikl,1)-TKE_AT(ikl,k1p(1))) |
---|
| 130 | ! #De S3D__1(1) = 0.0 |
---|
| 131 | ! #De S3D__2(1) = 1.0 |
---|
| 132 | ! #De S3D__4(1) = TKEtop |
---|
| 133 | |
---|
| 134 | DO k=k1p(1),mzp-2 |
---|
| 135 | S3D__4(k) = & |
---|
| 136 | & S3D__1(k) *a_b_AT*(TKE_AT(ikl, k )-TKE_AT(ikl,k1p(1))) & |
---|
| 137 | & -S3D__3(k) *a_b_AT*(TKE_AT(ikl,k1m(k))-TKE_AT(ikl, k )) |
---|
| 138 | END DO |
---|
| 139 | |
---|
| 140 | k= mzp-1 |
---|
| 141 | S3D__4(k) = & |
---|
| 142 | & S3D__1(k)*(a_b_AT* TKE_AT(ikl, k )-TKE_AT(ikl,k1p(1)) & |
---|
| 143 | & /betaAT ) & |
---|
| 144 | & -S3D__3(k) *a_b_AT*(TKE_AT(ikl,k1m(k))-TKE_AT(ikl, k )) |
---|
| 145 | |
---|
| 146 | S3D__1(k) = 0.000 |
---|
| 147 | |
---|
| 148 | ! S3D__4(mzp-1)=-(alphAT* TKE_AT(ikl,mzp-1)-TKE_AT(ikl,mzp )) & |
---|
| 149 | ! & * GravF2 *0.5000*(Kzm_AT(ikl,mzp-1)+Kzm_AT(ikl,mzp-2)) & |
---|
| 150 | ! & * roamDY(ikl,mzp-1)*roamDY(ikl,mzp-1) & |
---|
| 151 | ! & / (psa_sq* dsigmi( mzp-1)*dsigma( mzp-1)) & |
---|
| 152 | ! & -S3D__3(mzp-1)* a_b_AT*(TKE_AT(ikl,mzp-2)-TKE_AT(ikl,mzp-1)) |
---|
| 153 | |
---|
| 154 | |
---|
| 155 | ! Tridiagonal Matrix Inversion - TKE_AT |
---|
| 156 | ! ------------------------------------- |
---|
| 157 | |
---|
| 158 | k1= 1 |
---|
| 159 | ! #De k1= 2 |
---|
| 160 | DO k=k1,mzp-1 |
---|
| 161 | S3D__4(k) = S3D__4(k) + TKE_AT(ikl,k) |
---|
| 162 | END DO |
---|
| 163 | |
---|
| 164 | ! Forward Sweep |
---|
| 165 | ! ~~~~~~~~~~~~~~ |
---|
| 166 | S3D__5(1) = S3D__2(1) |
---|
| 167 | S3D__6(1) =-S3D__1(1) /S3D__5(1) |
---|
| 168 | DO k=k1p(1),mzp-1 |
---|
| 169 | S3D__5(k) = S3D__3(k) *S3D__6(k-1)+S3D__2(k) |
---|
| 170 | S3D__6(k) =-S3D__1(k) /S3D__5(k) |
---|
| 171 | END DO |
---|
| 172 | S3D__7(1) = S3D__4(1) /S3D__5(1) |
---|
| 173 | DO k=k1p(1),mzp-1 |
---|
| 174 | S3D__7(k) =(S3D__4(k) -S3D__3(k) & |
---|
| 175 | & *S3D__7(k-1))/S3D__5(k) |
---|
| 176 | END DO |
---|
| 177 | |
---|
| 178 | ! Backward Sweep |
---|
| 179 | ! ~~~~~~~~~~~~~~ |
---|
| 180 | DO k=k1m(mzp-1),1,-1 |
---|
| 181 | S3D__7(k) = S3D__6(k) *S3D__7(k+1)+S3D__7(k) |
---|
| 182 | END DO |
---|
| 183 | |
---|
| 184 | |
---|
| 185 | DO k=1,mzp-1 |
---|
| 186 | TrT_AT(ikl,k) = TrT_AT(ikl,k) & |
---|
| 187 | & +(S3D__7(k) -TKE_AT(ikl,k)) /dt__AT |
---|
| 188 | TKE_AT(ikl,k) = S3D__7(k) |
---|
| 189 | END DO |
---|
| 190 | |
---|
| 191 | |
---|
| 192 | ! Vertical Diffusion of Dissipation |
---|
| 193 | ! ================================= |
---|
| 194 | |
---|
| 195 | |
---|
| 196 | ! Update Tridiagonal Matrix Coefficients - eps_AT |
---|
| 197 | ! ----------------------------------------------- |
---|
| 198 | |
---|
| 199 | S3D__1(mzp-1) = S3DSBC * dt__AT ! SBC: TKE and epsilon, Atm SBL |
---|
| 200 | DO k=1,mzp-1 |
---|
| 201 | S3D__1(k) = S3D__1(k) * sige2k |
---|
| 202 | S3D__3(k) = S3D__3(k) * sige2k |
---|
| 203 | S3D__2(k) = 1.0 - S3D__3(k) - S3D__1(k) |
---|
| 204 | END DO |
---|
| 205 | |
---|
| 206 | |
---|
| 207 | ! Second Member of the Tridiagonal System - eps_AT |
---|
| 208 | ! ------------------------------------------------ |
---|
| 209 | |
---|
| 210 | S3D__4(1) = & |
---|
| 211 | & S3D__1(1) *a_b_AT*(eps_AT(ikl,1)-eps_AT(ikl,k1p(1))) |
---|
| 212 | ! #De S3D__1(1) = 0.0 |
---|
| 213 | ! #De S3D__2(1) = 1.0 |
---|
| 214 | ! #De S3D__4(1) = eps_DI(i,j) |
---|
| 215 | |
---|
| 216 | DO k=k1p(1),mzp-2 |
---|
| 217 | S3D__4(k) = & |
---|
| 218 | & S3D__1(k) *a_b_AT*(eps_AT(ikl,k)-eps_AT(ikl,k1p(k))) & |
---|
| 219 | & -S3D__3(k) *a_b_AT*(eps_AT(ikl,k1m(k))-eps_AT(ikl,k)) |
---|
| 220 | END DO |
---|
| 221 | |
---|
| 222 | k= mzp-1 |
---|
| 223 | S3D__4(k) = & |
---|
| 224 | & S3D__1(k)*(a_b_AT* eps_AT(ikl, k )-eps_AT(ikl,k1p(1)) & |
---|
| 225 | & /betaAT ) & |
---|
| 226 | & -S3D__3(k) *a_b_AT*(eps_AT(ikl,k1m(k))-eps_AT(ikl, k )) |
---|
| 227 | |
---|
| 228 | S3D__1(k) = 0.000 |
---|
| 229 | |
---|
| 230 | ! S3D__4(mzp-1)=-(alphAT* eps_AT(ikl,mzp-1)-eps_AT(ikl,mzp)) & |
---|
| 231 | ! & * GravF2* 0.5000*(Kzm_AT(ikl,mzp-1)+Kzm_AT(ikl,mzp-2)) & |
---|
| 232 | ! & * roamDY(ikl,mzp-1)*roamDY(ikl,mzp-1) & |
---|
| 233 | ! & / (psa_sq* dsigmi( mzp-1)*dsigma( mzp-1)) & |
---|
| 234 | ! & -S3D__3(mzp-1)* a_b_AT*(eps_AT(ikl,mzp-2)-eps_AT(ikl,mzp-1)) |
---|
| 235 | |
---|
| 236 | |
---|
| 237 | ! Tridiagonal Matrix Inversion - eps_AT |
---|
| 238 | ! ------------------------------------- |
---|
| 239 | |
---|
| 240 | k1= 1 |
---|
| 241 | ! #De k1= 2 |
---|
| 242 | DO k=k1,mzp-1 |
---|
| 243 | S3D__4(k) = S3D__4(k) + eps_AT(ikl,k) |
---|
| 244 | END DO |
---|
| 245 | |
---|
| 246 | ! Forward Sweep |
---|
| 247 | ! ~~~~~~~~~~~~~~ |
---|
| 248 | S3D__5(1) = S3D__2(1) |
---|
| 249 | S3D__6(1) =-S3D__1(1) /S3D__5(1) |
---|
| 250 | DO k=k1p(1),mzp-1 |
---|
| 251 | S3D__5(k) = S3D__3(k) *S3D__6(k-1)+S3D__2(k) |
---|
| 252 | S3D__6(k) =-S3D__1(k) /S3D__5(k) |
---|
| 253 | END DO |
---|
| 254 | S3D__7(1) = S3D__4(1) /S3D__5(1) |
---|
| 255 | DO k=k1p(1),mzp-1 |
---|
| 256 | S3D__7(k) =(S3D__4(k) -S3D__3(k) & |
---|
| 257 | & *S3D__7(k-1))/S3D__5(k) |
---|
| 258 | END DO |
---|
| 259 | |
---|
| 260 | ! Backward Sweep |
---|
| 261 | ! ~~~~~~~~~~~~~~ |
---|
| 262 | DO k=k1m(mzp-1),1,-1 |
---|
| 263 | S3D__7(k) = S3D__6(k) *S3D__7(k+1)+S3D__7(k) |
---|
| 264 | END DO |
---|
| 265 | |
---|
| 266 | DO k=1,mzp-1 |
---|
| 267 | eps_AT(ikl,k) = S3D__7(k) |
---|
| 268 | END DO |
---|
| 269 | |
---|
| 270 | |
---|
| 271 | |
---|
| 272 | |
---|
| 273 | ! ===================== |
---|
| 274 | ENDDO ! ikl = 1,kcolp |
---|
| 275 | ! ===================== |
---|
| 276 | |
---|
| 277 | |
---|
| 278 | |
---|
| 279 | |
---|
| 280 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 281 | ! ! |
---|
| 282 | ! DE-ALLOCATION |
---|
| 283 | ! ============= |
---|
| 284 | |
---|
| 285 | IF (FlagDALLOC) THEN ! |
---|
| 286 | deallocate ( S3D__1 ) |
---|
| 287 | deallocate ( S3D__2 ) |
---|
| 288 | deallocate ( S3D__3 ) |
---|
| 289 | deallocate ( S3D__4 ) |
---|
| 290 | deallocate ( S3D__5 ) |
---|
| 291 | deallocate ( S3D__6 ) |
---|
| 292 | deallocate ( S3D__7 ) |
---|
| 293 | ENDIF |
---|
| 294 | ! ! |
---|
| 295 | !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
| 296 | |
---|
| 297 | |
---|
| 298 | |
---|
| 299 | |
---|
| 300 | return |
---|
| 301 | end subroutine PHY_vdfTKE_RUN |
---|