- Timestamp:
- 2018-07-13T09:28:50+02:00 (6 years ago)
- Location:
- NEMO/branches/2018/dev_r9838_ENHANCE04_RK3
- Files:
-
- 1 edited
- 1 copied
-
. (copied) (copied from NEMO/trunk)
-
src/OCE/DYN/dynzdf.F90 (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DYN/dynzdf.F90
r9598 r9939 11 11 !!---------------------------------------------------------------------- 12 12 !! dyn_zdf : compute the after velocity through implicit calculation of vertical mixing 13 !! zdf_trd : diagnose the zdf velocity trends and the KE dissipation trend 14 !!gm ==>>> zdf_trd currently not used 13 15 !!---------------------------------------------------------------------- 14 16 USE oce ! ocean dynamics and tracers variables … … 26 28 USE in_out_manager ! I/O manager 27 29 USE lib_mpp ! MPP library 30 USE iom ! IOM library 28 31 USE prtctl ! Print control 29 32 USE timing ! Timing … … 67 70 INTEGER, INTENT(in) :: kt ! ocean time-step index 68 71 ! 69 INTEGER :: ji, jj, jk ! dummy loop indices70 INTEGER :: iku, ikv ! local integers71 REAL(wp) :: zzwi, ze3ua, z dt! local scalars72 REAL(wp) :: zzws, ze3va ! - -73 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwd, zws ! 3D workspace74 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - -72 INTEGER :: ji, jj, jk ! dummy loop indices 73 INTEGER :: iku, ikv ! local integers 74 REAL(wp) :: zzwi, ze3ua, z2dt_2 ! local scalars 75 REAL(wp) :: zzws, ze3va ! - - 76 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwd, zws ! 3D workspace 77 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - - 75 78 !!--------------------------------------------------------------------- 76 79 ! … … 86 89 ENDIF 87 90 ENDIF 88 ! !* set time step 89 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! = rdt (restart with Euler time stepping) 90 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! = 2 rdt (leapfrog) 91 ENDIF 91 ! 92 z2dt_2 = rDt * 0.5_wp !* =rn_Dt except in 1st Euler time step where it is equal to rn_Dt/2 93 ! 92 94 ! 93 95 ! !* explicit top/bottom drag case … … 106 108 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 107 109 DO jk = 1, jpkm1 108 ua(:,:,jk) = ( ub(:,:,jk) + r 2dt * ua(:,:,jk) ) * umask(:,:,jk)109 va(:,:,jk) = ( vb(:,:,jk) + r 2dt * va(:,:,jk) ) * vmask(:,:,jk)110 ua(:,:,jk) = ( ub(:,:,jk) + rDt * ua(:,:,jk) ) * umask(:,:,jk) 111 va(:,:,jk) = ( vb(:,:,jk) + rDt * va(:,:,jk) ) * vmask(:,:,jk) 110 112 END DO 111 113 ELSE ! applied on thickness weighted velocity 112 114 DO jk = 1, jpkm1 113 ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) & 114 & + r2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 115 va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) & 116 & + r2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 115 ua(:,:,jk) = ( e3u_b(:,:,jk)*ub(:,:,jk) + rDt * e3u_n(:,:,jk)*ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 116 va(:,:,jk) = ( e3v_b(:,:,jk)*vb(:,:,jk) + rDt * e3v_n(:,:,jk)*va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 117 117 END DO 118 118 ENDIF … … 133 133 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 134 134 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 135 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua136 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va135 ua(ji,jj,iku) = ua(ji,jj,iku) + z2dt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 136 va(ji,jj,ikv) = va(ji,jj,ikv) + z2dt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 137 137 END DO 138 138 END DO … … 144 144 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 145 145 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 146 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua147 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va146 ua(ji,jj,iku) = ua(ji,jj,iku) + z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 147 va(ji,jj,ikv) = va(ji,jj,ikv) + z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 148 148 END DO 149 149 END DO … … 153 153 ! !== Vertical diffusion on u ==! 154 154 ! 155 ! !* Matrix construction 156 zdt = r2dt * 0.5 157 SELECT CASE( nldf_dyn ) 158 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 155 SELECT CASE( nldf_dyn ) !* Matrix construction 156 ! 157 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 159 158 DO jk = 1, jpkm1 160 159 DO jj = 2, jpjm1 161 160 DO ji = fs_2, fs_jpim1 ! vector opt. 162 161 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 163 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk) + akzu(ji,jj,jk ) ) &164 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk )165 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &166 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)162 zzwi = - rDt * ( 0.5 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) + akzu(ji,jj,jk ) ) & 163 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 164 zzws = - rDt * ( 0.5 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) + akzu(ji,jj,jk+1) ) & 165 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 167 166 zwi(ji,jj,jk) = zzwi 168 167 zws(ji,jj,jk) = zzws … … 171 170 END DO 172 171 END DO 173 CASE DEFAULT ! iso-level lateral mixing172 CASE DEFAULT ! iso-level lateral mixing 174 173 DO jk = 1, jpkm1 175 174 DO jj = 2, jpjm1 176 175 DO ji = fs_2, fs_jpim1 ! vector opt. 177 176 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 178 zzwi = - z dt* ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk )179 zzws = - z dt* ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)177 zzwi = - z2dt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 178 zzws = - z2dt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 180 179 zwi(ji,jj,jk) = zzwi 181 180 zws(ji,jj,jk) = zzws … … 186 185 END SELECT 187 186 ! 188 DO jj = 2, jpjm1 !* Surface boundary conditions189 DO ji = fs_2, fs_jpim1 ! vector opt.187 DO jj = 2, jpjm1 !* Surface boundary conditions 188 DO ji = fs_2, fs_jpim1 ! vector opt. 190 189 zwi(ji,jj,1) = 0._wp 191 190 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) … … 204 203 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 205 204 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 206 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua205 zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 207 206 END DO 208 207 END DO … … 213 212 iku = miku(ji,jj) ! ocean top level at u- and v-points 214 213 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 215 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua214 zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 216 215 END DO 217 216 END DO … … 245 244 DO ji = fs_2, fs_jpim1 ! vector opt. 246 245 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 247 ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 248 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 246 ua(ji,jj,1) = ua(ji,jj,1) + z2dt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( ze3ua * rho0 ) * umask(ji,jj,1) 249 247 END DO 250 248 END DO … … 272 270 ! !== Vertical diffusion on v ==! 273 271 ! 274 ! !* Matrix construction 275 zdt = r2dt * 0.5 276 SELECT CASE( nldf_dyn ) 277 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 272 ! 273 SELECT CASE( nldf_dyn ) !* Matrix construction 274 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 278 275 DO jk = 1, jpkm1 279 276 DO jj = 2, jpjm1 280 277 DO ji = fs_2, fs_jpim1 ! vector opt. 281 278 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 282 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk) + akzv(ji,jj,jk ) ) &283 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk )284 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &285 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)279 zzwi = - rDt * ( 0.5 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) + akzv(ji,jj,jk ) ) & 280 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 281 zzws = - rDt * ( 0.5 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) + akzv(ji,jj,jk+1) ) & 282 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 286 283 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 287 284 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) … … 290 287 END DO 291 288 END DO 292 CASE DEFAULT ! iso-level lateral mixing289 CASE DEFAULT ! iso-level lateral mixing 293 290 DO jk = 1, jpkm1 294 291 DO jj = 2, jpjm1 295 292 DO ji = fs_2, fs_jpim1 ! vector opt. 296 293 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 297 zzwi = - z dt* ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk )298 zzws = - z dt* ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)294 zzwi = - z2dt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 295 zzws = - z2dt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 299 296 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 300 297 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) … … 305 302 END SELECT 306 303 ! 307 DO jj = 2, jpjm1 !* Surface boundary conditions308 DO ji = fs_2, fs_jpim1 ! vector opt.304 DO jj = 2, jpjm1 !* Surface boundary conditions 305 DO ji = fs_2, fs_jpim1 ! vector opt. 309 306 zwi(ji,jj,1) = 0._wp 310 307 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) … … 322 319 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 323 320 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 324 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va321 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - z2dt_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 325 322 END DO 326 323 END DO … … 330 327 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 331 328 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 332 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va329 zwd(ji,jj,iku) = zwd(ji,jj,iku) - z2dt_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 333 330 END DO 334 331 END DO … … 362 359 DO ji = fs_2, fs_jpim1 ! vector opt. 363 360 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 364 va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 365 & / ( ze3va * rau0 ) * vmask(ji,jj,1) 361 va(ji,jj,1) = va(ji,jj,1) + z2dt_2 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( ze3va * rho0 ) * vmask(ji,jj,1) 366 362 END DO 367 363 END DO … … 387 383 END DO 388 384 ! 389 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 390 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 391 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 385 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 386 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 387 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_Dt - ztrdu(:,:,:) 388 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_Dt - ztrdv(:,:,:) 389 ELSE ! applied on thickness weighted velocity 390 ztrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_Dt - ztrdu(:,:,:) 391 ztrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_Dt - ztrdv(:,:,:) 392 ENDIF 392 393 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 393 394 DEALLOCATE( ztrdu, ztrdv ) … … 401 402 END SUBROUTINE dyn_zdf 402 403 404 !!gm currently not used : just for memory to be able to add dissipation trend through vertical mixing 405 406 SUBROUTINE zdf_trd( ptrdu, ptrdv ,kt ) 407 !!---------------------------------------------------------------------- 408 !! *** ROUTINE zdf_trd *** 409 !! 410 !! ** Purpose : compute the trend due to the vert. momentum diffusion 411 !! together with the Leap-Frog time stepping using an 412 !! implicit scheme. 413 !! 414 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 415 !! ua = ub + 2*dt * ua vector form or linear free surf. 416 !! ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a otherwise 417 !! - update the after velocity with the implicit vertical mixing. 418 !! This requires to solver the following system: 419 !! ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ] 420 !! with the following surface/top/bottom boundary condition: 421 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) 422 !! top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 423 !! 424 !! ** Action : (ua,va) after velocity 425 !!--------------------------------------------------------------------- 426 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: ptrdu, ptrdv ! 3D work arrays use for zdf trends diag 427 INTEGER , INTENT(in ) :: kt ! ocean time-step index 428 ! 429 INTEGER :: ji, jj, jk ! dummy loop indices 430 REAL(wp) :: zzz ! local scalar 431 REAL(wp) :: zavmu, zavmum1 ! - - 432 REAL(wp) :: zavmv, zavmvm1 ! - - 433 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z2d ! - - 434 !!--------------------------------------------------------------------- 435 ! 436 CALL lbc_lnk_multi( ua, 'U', -1., va, 'V', -1. ) ! apply lateral boundary condition on (ua,va) 437 ! 438 ! 439 ! !== momentum trend due to vertical diffusion ==! 440 ! 441 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 442 ptrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * r1_Dt - ptrdu(:,:,:) 443 ptrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * r1_Dt - ptrdv(:,:,:) 444 ELSE ! applied on thickness weighted velocity 445 ptrdu(:,:,:) = ( e3u_a(:,:,:)*ua(:,:,:) - e3u_b(:,:,:)*ub(:,:,:) ) / e3u_n(:,:,:) * r1_Dt - ptrdu(:,:,:) 446 ptrdv(:,:,:) = ( e3v_a(:,:,:)*va(:,:,:) - e3v_b(:,:,:)*vb(:,:,:) ) / e3v_n(:,:,:) * r1_Dt - ptrdv(:,:,:) 447 ENDIF 448 CALL trd_dyn( ptrdu, ptrdv, jpdyn_zdf, kt ) 449 ! 450 ! 451 ! !== KE dissipation trend due to vertical diffusion ==! 452 ! 453 IF( iom_use( 'dispkevfo' ) ) THEN ! ocean kinetic energy dissipation per unit area 454 ! ! due to v friction (v=vertical) 455 ! ! see NEMO_book appendix C, §C.8 (N.B. here averaged at t-points) 456 ! ! Note that formally, in a Leap-Frog environment, the shear**2 should be the product of 457 ! ! now by before shears, i.e. the source term of TKE (local positivity is not ensured). 458 ! ! Note also that now e3 scale factors are used as after one are not computed ! 459 ! 460 CALL wrk_alloc(jpi,jpj, z2d ) 461 z2d(:,:) = 0._wp 462 DO jk = 1, jpkm1 463 DO jj = 2, jpjm1 464 DO ji = 2, jpim1 465 zavmu = 0.5 * ( avm(ji+1,jj,jk) + avm(ji ,jj,jk) ) 466 zavmum1 = 0.5 * ( avm(ji ,jj,jk) + avm(ji-1,jj,jk) ) 467 zavmv = 0.5 * ( avm(ji,jj+1,jk) + avm(ji,jj ,jk) ) 468 zavmvm1 = 0.5 * ( avm(ji,jj ,jk) + avm(ji,jj-1,jk) ) 469 470 z2d(ji,jj) = z2d(ji,jj) + ( & 471 & zavmu * ( ua(ji ,jj,jk-1) - ua(ji ,jj,jk) )**2 / e3uw_n(ji ,jj,jk) * wumask(ji ,jj,jk) & 472 & + zavmum1 * ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) )**2 / e3uw_n(ji-1,jj,jk) * wumask(ji-1,jj,jk) & 473 & + zavmv * ( va(ji,jj ,jk-1) - va(ji,jj ,jk) )**2 / e3vw_n(ji,jj ,jk) * wvmask(ji,jj ,jk) & 474 & + zavmvm1 * ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) )**2 / e3vw_n(ji,jj-1,jk) * wvmask(ji,jj-1,jk) & 475 & ) 476 !!gm --- This trends is in fact properly computed in zdf_sh2 but with a backward shift of one time-step ===>>> use it ? 477 !! No since in zdfshé only kz tke (or gls) is used 478 !! 479 !!gm --- formally, as done below, in a Leap-Frog environment, the shear**2 should be the product of 480 !!gm now by before shears, i.e. the source term of TKE (local positivity is not ensured). 481 !! CAUTION: requires to compute e3uw_a and e3vw_a !!! 482 ! z2d(ji,jj) = z2d(ji,jj) + ( & 483 ! & avmu(ji ,jj,jk) * ( un(ji ,jj,jk-1) - un(ji ,jj,jk) ) / e3uw_n(ji ,jj,jk) & 484 ! & * ( ua(ji ,jj,jk-1) - ua(ji ,jj,jk) ) / e3uw_a(ji ,jj,jk) * wumask(ji ,jj,jk) & 485 ! & + avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) / e3uw_n(ji-1,jj,jk) & 486 ! & ( ua(ji-1,jj,jk-1) - ua(ji-1,jj,jk) ) / e3uw_a(ji-1,jj,jk) * wumask(ji-1,jj,jk) & 487 ! & + avmv(ji,jj ,jk) * ( vn(ji,jj ,jk-1) - vn(ji,jj ,jk) ) / e3vw_n(ji,jj ,jk) & 488 ! & ( va(ji,jj ,jk-1) - va(ji,jj ,jk) ) / e3vw_a(ji,jj ,jk) * wvmask(ji,jj ,jk) & 489 ! & + avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) / e3vw_n(ji,jj-1,jk) & 490 ! & ( va(ji,jj-1,jk-1) - va(ji,jj-1,jk) ) / e3vw_a(ji,jj-1,jk) * wvmask(ji,jj-1,jk) & 491 ! & ) 492 !!gm end 493 END DO 494 END DO 495 END DO 496 zzz= - 0.5_wp* rho0 ! caution sign minus here 497 z2d(:,:) = zzz * z2d(:,:) 498 CALL lbc_lnk( z2d,'T', 1. ) 499 CALL iom_put( 'dispkevfo', z2d ) 500 ENDIF 501 ! 502 END SUBROUTINE zdf_trd 503 504 !!gm end not used 505 403 506 !!============================================================================== 404 507 END MODULE dynzdf
Note: See TracChangeset
for help on using the changeset viewer.
