Changeset 1481 for trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90
- Timestamp:
- 2009-07-03T17:07:08+02:00 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90
r1268 r1481 50 50 51 51 PUBLIC zdf_tke2 ! routine called in step module 52 PUBLIC tke2_rst ! routine called in step module 52 53 53 54 LOGICAL , PUBLIC, PARAMETER :: lk_zdftke2 = .TRUE. !: TKE vertical mixing flag 54 55 REAL(wp), PUBLIC :: eboost !: multiplicative coeff of the shear product. 55 56 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: en !: now turbulent kinetic energy 56 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: avm t!: now mixing coefficient of diffusion for TKE57 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: avm !: now mixing coefficient of diffusion for TKE 57 58 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: dissl !: now mixing lenght of dissipation 58 59 … … 166 167 INTEGER, INTENT(in) :: kt ! ocean time step 167 168 168 IF( kt == nit000 ) THEN 169 CALL zdf_tke2_init ! Initialization (first time-step only) 170 IF( .NOT. ln_rstart ) CALL tke_tke( kt ) 171 CALL tke_mlmc 172 ELSE 173 CALL tke_tke( kt ) 174 IF( kt /= nitend+1 ) CALL tke_mlmc 175 IF( lrst_oce ) CALL tke2_rst( kt, 'WRITE' ) ! write en in restart file 176 ENDIF 169 IF( kt == nit000 ) CALL zdf_tke2_init ! Initialization (first time-step only) 170 ! 171 CALL tke_tke ! now tke (en) 172 CALL tke_mlmc ! now avmu, avmv, avt 177 173 178 174 END SUBROUTINE zdf_tke2 179 175 180 SUBROUTINE tke_tke ( kt )176 SUBROUTINE tke_tke 181 177 !!---------------------------------------------------------------------- 182 178 !! Now Turbulent kinetic energy (output in en) … … 185 181 !! Resolution of a tridiagonal linear system by a "methode de chasse" 186 182 !! computation from level 2 to jpkm1 (e(1) computed and e(jpk)=0 ). 187 !! avm trepresents the turbulent coefficient of the TKE183 !! avm represents the turbulent coefficient of the TKE 188 184 !!---------------------------------------------------------------------- 189 185 USE oce, zwd => ua ! use ua as workspace … … 191 187 USE oce, zdiag2 => ta ! use ta as workspace 192 188 USE oce, ztkelc => sa ! use sa as workspace 193 !194 INTEGER, INTENT(in) :: kt ! ocean time step195 189 ! 196 190 INTEGER :: ji, jj, jk ! dummy loop arguments … … 209 203 zfact1 = -.5 * rdt * rn_efave 210 204 zfact2 = 1.5 * rdt * rn_ediss 211 zfact3 = 0.5 * rdt * rn_ediss 205 zfact3 = 0.5 * rn_ediss 206 207 ! Surface boundary condition on tke 208 ! ------------------------------------------------- 209 ! en(1) = rn_ebb sqrt(utau^2+vtau^2) / rau0 (min value rn_emin0) 210 !CDIR NOVERRCHK 211 DO jj = 2, jpjm1 212 !CDIR NOVERRCHK 213 DO ji = fs_2, fs_jpim1 ! vector opt. 214 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 215 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 216 zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 217 en(ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1) 218 END DO 219 END DO 212 220 213 221 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 217 225 ! 218 226 ! Computation of total energy produce by LC : cumulative sum over jk 219 zpelc(:,:,1) = MAX( rn2 (:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)227 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1) 220 228 DO jk = 2, jpk 221 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2 (:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)229 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk) 222 230 END DO 223 231 ! … … 246 254 # endif 247 255 ! 248 ! TKE Langmuir circulation source term 256 ! TKE Langmuir circulation source term added to en 249 257 !CDIR NOVERRCHK 250 258 DO jk = 2, jpkm1 … … 259 267 zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 260 268 ! TKE Langmuir circulation source term 261 ztkelc(ji,jj,jk) = ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)269 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) ) * tmask(ji,jj,jk) 262 270 END DO 263 271 END DO … … 266 274 ENDIF 267 275 268 ! Surface boundary condition on tke 269 ! ------------------------------------------------- 270 ! en(1) = rn_ebb sqrt(utau^2+vtau^2) / rau0 (min value rn_emin0) 271 !CDIR NOVERRCHK 272 DO jj = 2, jpjm1 273 !CDIR NOVERRCHK 274 DO ji = fs_2, fs_jpim1 ! vector opt. 275 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 276 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 277 zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 278 en (ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1) 279 END DO 280 END DO 281 282 ! Preliminary computing 276 277 ! Shear production at uw- and vw-points (energy conserving form) 283 278 DO jk = 2, jpkm1 284 279 DO jj = 2, jpjm1 285 280 DO ji = fs_2, fs_jpim1 ! vector opt. 286 281 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) * & 287 & ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) 282 & ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) / & 283 & ( fse3uw_n(ji,jj,jk) * fse3uw_b(ji,jj,jk) ) 288 284 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) * & 289 & ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) 285 & ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) / & 286 & ( fse3vw_n(ji,jj,jk) * fse3vw_b(ji,jj,jk) ) 290 287 ENDDO 291 288 ENDDO … … 304 301 DO jj = 2, jpjm1 305 302 DO ji = fs_2, fs_jpim1 ! vector opt. 306 ! ! shear prod. - stratif. destruction 307 zcoef = 0.5 / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) ) 308 ! shear 309 zdku = zcoef * ( avmu(ji-1,jj ,jk) + avmu(ji,jj,jk) ) ! shear 310 zdkv = zcoef * ( avmv(ji ,jj-1,jk) + avmv(ji,jj,jk) ) 311 zesh2 = eboost * ( zdku + zdkv ) - avt(ji,jj,jk) * rn2(ji,jj,jk) ! coefficient (zesh2) 312 ! 303 ! ! shear prod. at w-point weightened by mask 304 zesh2 = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 305 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 313 306 ! ! Matrix 314 307 zcof = zfact1 * tmask(ji,jj,jk) 315 308 ! ! lower diagonal 316 zdiag2(ji,jj,jk) = zcof * ( avm t (ji,jj,jk ) + avmt(ji,jj,jk-1) ) &309 zdiag2(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & 317 310 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) ) 318 311 ! ! upper diagonal 319 zdiag1(ji,jj,jk) = zcof * ( avm t (ji,jj,jk+1) + avmt(ji,jj,jk ) ) &312 zdiag1(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & 320 313 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) ) 321 314 ! ! diagonal … … 323 316 ! 324 317 ! ! right hand side in en 325 IF( .NOT. ln_lc ) THEN ! No Langmuir cells 326 en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * dissl (ji,jj,jk) * en(ji,jj,jk) * tmask(ji,jj,jk) & 327 & + rdt * zesh2 * tmask(ji,jj,jk) 328 ELSE ! Langmuir cell source term 329 en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * dissl (ji,jj,jk) * en(ji,jj,jk) * tmask(ji,jj,jk) & 330 & + rdt * zesh2 * tmask(ji,jj,jk) & 331 & + rdt * ztkelc(ji,jj,jk) * tmask(ji,jj,jk) 332 ENDIF 318 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 319 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) * tmask(ji,jj,jk) 333 320 END DO 334 321 END DO … … 434 421 USE oce, zmxlm => va ! use va as workspace 435 422 USE oce, zmxld => ta ! use ta as workspace 436 USE oce, visco => sa ! use sa as workspace437 423 ! 438 424 INTEGER :: ji, jj, jk ! dummy loop arguments 439 REAL(wp) :: z bbrau, zrn2, zraug! temporary scalars440 REAL(wp) :: z fact1, ztx2, zdku! - -441 REAL(wp) :: z fact2, zty2, zdkv! - -442 REAL(wp) :: z fact3, zcoef, zav! - -425 REAL(wp) :: zrn2, zraug ! temporary scalars 426 REAL(wp) :: ztx2, zdku ! - - 427 REAL(wp) :: zty2, zdkv ! - - 428 REAL(wp) :: zcoef, zav ! - - 443 429 REAL(wp) :: zpdl, zri, zsqen ! - - 444 430 REAL(wp) :: zemxl, zemlm, zemlp ! - - … … 569 555 ! II Vertical eddy viscosity on tke (put in zmxlm) and first estimate of avt ! 570 556 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<! 571 ! Surface boundary condition on avt and avm t: jk = 1557 ! Surface boundary condition on avt and avm : jk = 1 572 558 ! ------------------------------------------------- 573 ! avt (1) = avmb(1) and avt(jpk) = 0.574 ! avm t(1) = avmb(1) and avmt(jpk) = 0.559 ! avt(1) = avmb(1) and avt(jpk) = 0. 560 ! avm(1) = avmb(1) and avm(jpk) = 0. 575 561 ! 576 562 !CDIR NOVERRCHK … … 582 568 zsqen = SQRT( en(ji,jj,jk) ) 583 569 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 584 avmt (ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk) 585 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 586 visco(ji,jj,jk) = avmt (ji,jj,jk) 570 avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk) 571 avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 587 572 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 588 573 END DO … … 591 576 592 577 ! Lateral boundary conditions (sign unchanged) 593 CALL lbc_lnk( avm t, 'W', 1. ) ; CALL lbc_lnk( visco, 'W', 1. )578 CALL lbc_lnk( avm, 'W', 1. ) 594 579 595 580 DO jk = 2, jpkm1 ! only vertical eddy viscosity 596 581 DO jj = 2, jpjm1 597 582 DO ji = fs_2, fs_jpim1 ! vector opt. 598 avmu(ji,jj,jk) = ( visco(ji,jj,jk) + visco(ji+1,jj ,jk) ) * umask(ji,jj,jk) & 599 & / MAX( 1., tmask(ji,jj,jk) + tmask(ji+1,jj ,jk) ) 600 avmv(ji,jj,jk) = ( visco(ji,jj,jk) + visco(ji ,jj+1,jk) ) * vmask(ji,jj,jk) & 601 & / MAX( 1., tmask(ji,jj,jk) + tmask(ji ,jj+1,jk) ) 583 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 584 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) 602 585 END DO 603 586 END DO 604 587 END DO 605 588 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 606 ! 607 ! 608 DO jk = 2, jpkm1 ! Minimum value on the eddy viscosity 609 avmu(:,:,jk) = MAX( avmu(:,:,jk), avmb(jk) ) * umask(:,:,jk) 610 avmv(:,:,jk) = MAX( avmv(:,:,jk), avmb(jk) ) * vmask(:,:,jk) 611 END DO 589 612 590 613 591 ! Prandtl number … … 623 601 zdkv = avmv(ji,jj-1,jk) * (vn(ji,jj-1,jk-1)-vn(ji,jj-1,jk)) * (vb(ji,jj-1,jk-1)-vb(ji,jj-1,jk)) + & 624 602 & avmv(ji,jj ,jk) * (vn(ji,jj ,jk-1)-vn(ji,jj ,jk)) * (vb(ji,jj ,jk-1)-vb(ji,jj ,jk)) 625 zri = MAX( rn2 (ji,jj,jk), 0. ) / (zcoef * (zdku + zdkv + 1.e-20) / avmt(ji,jj,jk)) ! local Richardson number603 zri = MAX( rn2b(ji,jj,jk), 0. ) * avm(ji,jj,jk)/ (zcoef * (zdku + zdkv + 1.e-20) ) ! local Richardson number 626 604 # if defined key_cfg_1d 627 605 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri … … 770 748 DO jk = 1, jpk 771 749 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 750 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 772 751 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 773 752 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 774 avmt(:,:,jk) = avmb(jk) * tmask(:,:,jk)775 753 END DO 776 754 dissl(:,:,:) = 1.e-12 … … 796 774 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 797 775 ! 798 INTEGER :: jit ! dummy loop indices 776 INTEGER :: jit, jk ! dummy loop indices 777 INTEGER :: id1, id2, id3, id4, id5, id6 799 778 !!---------------------------------------------------------------------- 800 779 ! 801 IF( TRIM(cdrw) == 'READ' ) THEN 802 IF( ln_rstart ) THEN 803 IF( iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 .AND. .NOT.(ln_rstke) ) THEN 780 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 781 ! ! --------------- 782 IF( ln_rstart ) THEN !* Read the restart file 783 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 784 id2 = iom_varid( numror, 'avt' , ldstop = .FALSE. ) 785 id3 = iom_varid( numror, 'avm' , ldstop = .FALSE. ) 786 id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. ) 787 id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. ) 788 id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 789 ! 790 IF( id1 > 0 ) THEN ! 'en' exists 804 791 CALL iom_get( numror, jpdom_autoglo, 'en', en ) 805 ELSE 806 IF( lwp .AND. iom_varid( numror, 'en', ldstop = .FALSE. ) > 0 ) & 807 & WRITE(numout,*) ' ===>>>> : previous run without tke scheme' 808 IF( lwp .AND. ln_rstke ) WRITE(numout,*) ' ===>>>> : We do not use en from the restart file' 809 IF( lwp ) WRITE(numout,*) ' ===>>>> : en is set by iterative loop' 792 IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN ! all required arrays exist 793 CALL iom_get( numror, jpdom_autoglo, 'avt' , avt ) 794 CALL iom_get( numror, jpdom_autoglo, 'avm' , avm ) 795 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu ) 796 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv ) 797 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 798 ELSE ! one at least array is missing 799 CALL tke_mlmc ! recompute avt, avm, avmu, avmv and dissl (approximation) 800 ENDIF 801 ELSE ! No TKE array found: initialisation 802 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop' 810 803 en (:,:,:) = rn_emin * tmask(:,:,:) 811 DO jit = 2, nn_itke + 1 812 CALL zdf_tke2( jit ) 813 END DO 804 CALL tke_mlmc ! recompute avt, avm, avmu, avmv and dissl (approximation) 805 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke2( jit ) ; END DO 814 806 ENDIF 815 ELSE 816 en(:,:,:) = rn_emin * tmask(:,:,:) ! no restart: en set to its minimum value 807 ELSE !* Start from rest 808 en(:,:,:) = rn_emin * tmask(:,:,:) 809 DO jk = 1, jpk ! set the Kz to the background value 810 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk) 811 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk) 812 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk) 813 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 814 END DO 817 815 ENDIF 818 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 816 ! 817 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 818 ! ! ------------------- 819 819 IF(lwp) WRITE(numout,*) '---- tke2-rst ----' 820 CALL iom_rstput( kt, nitrst_tke2, numrow, 'en', en ) 820 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 821 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt ) 822 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm ) 823 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu ) 824 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv ) 825 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 826 ! 821 827 ENDIF 822 828 !
Note: See TracChangeset
for help on using the changeset viewer.