Changeset 15603 for branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
- Timestamp:
- 2021-12-16T10:39:55+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r12555 r15603 2 2 !!====================================================================== 3 3 !! *** MODULE zdftke *** 4 !! Ocean physics: vertical mixing coefficient computed from the tke 4 !! Ocean physics: vertical mixing coefficient computed from the tke 5 5 !! turbulent closure parameterization 6 6 !!===================================================================== … … 22 22 !! - ! 2008-05 (J.-M. Molines, G. Madec) 2D form of avtb 23 23 !! - ! 2008-06 (G. Madec) style + DOCTOR name for namelist parameters 24 !! - ! 2008-12 (G. Reffray) stable discretization of the production term 25 !! 3.2 ! 2009-06 (G. Madec, S. Masson) TKE restart compatible with key_cpl 24 !! - ! 2008-12 (G. Reffray) stable discretization of the production term 25 !! 3.2 ! 2009-06 (G. Madec, S. Masson) TKE restart compatible with key_cpl 26 26 !! ! + cleaning of the parameters + bugs correction 27 27 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase … … 52 52 USE wrk_nemo ! work arrays 53 53 USE timing ! Timing 54 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 54 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 55 55 #if defined key_agrif 56 56 USE agrif_opa_interp 57 57 USE agrif_opa_update 58 58 #endif 59 60 59 USE stopack 61 60 62 61 IMPLICIT NONE … … 75 74 INTEGER :: nn_pdl ! Prandtl number or not (ratio avt/avm) (=0/1) 76 75 REAL(wp) :: rn_ediff ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 77 REAL(wp) :: rn_ediss ! coefficient of the Kolmogoroff dissipation 76 REAL(wp) :: rn_ediss ! coefficient of the Kolmogoroff dissipation 78 77 REAL(wp) :: rn_ebb ! coefficient of the surface input of tke 79 78 REAL(wp) :: rn_emin ! minimum value of tke [m2/s2] … … 94 93 95 94 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 96 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_niw !: TKE budget- near-inertial waves term97 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: efr ! surface boundary condition for nn_etau = 498 95 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 99 96 #if defined key_c1d … … 102 99 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers 103 100 #endif 101 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 102 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_niw !: TKE budget- near-inertial waves term 103 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: efr ! surface boundary condition for nn_etau = 4 104 104 105 105 !! * Substitutions … … 118 118 !!---------------------------------------------------------------------- 119 119 ALLOCATE( & 120 & efr (jpi,jpj) , e_niw(jpi,jpj,jpk) , & 120 & efr (jpi,jpj) , e_niw(jpi,jpj,jpk) , & 121 121 #if defined key_c1d 122 122 & e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) , & … … 147 147 !! surface: en = max( rn_emin0, rn_ebb * taum ) 148 148 !! bottom : en = rn_emin 149 !! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 150 !! 151 !! The now Turbulent kinetic energy is computed using the following 149 !! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 150 !! 151 !! The now Turbulent kinetic energy is computed using the following 152 152 !! time stepping: implicit for vertical diffusion term, linearized semi 153 !! implicit for kolmogoroff dissipation term, and explicit forward for 154 !! both buoyancy and shear production terms. Therefore a tridiagonal 153 !! implicit for kolmogoroff dissipation term, and explicit forward for 154 !! both buoyancy and shear production terms. Therefore a tridiagonal 155 155 !! linear system is solved. Note that buoyancy and shear terms are 156 156 !! discretized in a energy conserving form (Bruchard 2002). … … 160 160 !! 161 161 !! The now vertical eddy vicosity and diffusivity coefficients are 162 !! given by: 162 !! given by: 163 163 !! avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 164 !! avt = max( avmb, pdl * avm ) 164 !! avt = max( avmb, pdl * avm ) 165 165 !! eav = max( avmb, avm ) 166 166 !! where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and 167 !! given by an empirical funtion of the localRichardson number if nn_pdl=1 167 !! given by an empirical funtion of the localRichardson number if nn_pdl=1 168 168 !! 169 169 !! ** Action : compute en (now turbulent kinetic energy) … … 180 180 ! 181 181 IF( kt /= nit000 ) THEN ! restore before value to compute tke 182 avt (:,:,:) = avt_k (:,:,:) 183 avm (:,:,:) = avm_k (:,:,:) 184 avmu(:,:,:) = avmu_k(:,:,:) 185 avmv(:,:,:) = avmv_k(:,:,:) 186 ENDIF 182 avt (:,:,:) = avt_k (:,:,:) 183 avm (:,:,:) = avm_k (:,:,:) 184 avmu(:,:,:) = avmu_k(:,:,:) 185 avmv(:,:,:) = avmv_k(:,:,:) 186 ENDIF 187 ! 188 #if defined key_traldf_c2d || key_traldf_c3d 189 IF( ln_stopack) THEN 190 IF( nn_spp_tkelc > 0 ) THEN 191 rn_lc0 = rn_lc 192 CALL spp_gen( kt, rn_lc0, nn_spp_tkelc, rn_tkelc_sd, jk_spp_tkelc ) 193 ENDIF 194 IF( nn_spp_tkedf > 0 ) THEN 195 rn_ediff0 = rn_ediff 196 CALL spp_gen( kt, rn_ediff0, nn_spp_tkedf, rn_tkedf_sd, jk_spp_tkedf ) 197 ENDIF 198 IF( nn_spp_tkeds > 0 ) THEN 199 rn_ediss0 = rn_ediss 200 CALL spp_gen( kt, rn_ediss0, nn_spp_tkeds, rn_tkeds_sd, jk_spp_tkeds ) 201 ENDIF 202 IF( nn_spp_tkebb > 0 ) THEN 203 rn_ebb0 = rn_ebb 204 CALL spp_gen( kt, rn_ebb0, nn_spp_tkebb, rn_tkebb_sd, jk_spp_tkebb ) 205 ENDIF 206 IF( nn_spp_tkefr > 0 ) THEN 207 rn_efr0 = rn_efr 208 CALL spp_gen( kt, rn_efr0, nn_spp_tkefr, rn_tkefr_sd, jk_spp_tkefr ) 209 ENDIF 210 ENDIF 211 #else 212 IF ( ln_stopack ) & 213 & CALL ctl_stop( 'zdf_tke: parameter perturbation will only work with '// & 214 'key_traldf_c2d or key_traldf_c3d') 215 #endif 187 216 ! 188 217 CALL tke_tke ! now tke (en) … … 190 219 CALL tke_avn ! now avt, avm, avmu, avmv 191 220 ! 192 avt_k (:,:,:) = avt (:,:,:) 193 avm_k (:,:,:) = avm (:,:,:) 194 avmu_k(:,:,:) = avmu(:,:,:) 195 avmv_k(:,:,:) = avmv(:,:,:) 221 avt_k (:,:,:) = avt (:,:,:) 222 avm_k (:,:,:) = avm (:,:,:) 223 avmu_k(:,:,:) = avmu(:,:,:) 224 avmv_k(:,:,:) = avmv(:,:,:) 196 225 ! 197 226 #if defined key_agrif 198 ! Update child grid f => parent grid 227 ! Update child grid f => parent grid 199 228 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only 200 #endif 201 ! 229 #endif 230 IF ( kt == nitend ) THEN 231 DEALLOCATE ( rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 ) 232 ENDIF 233 ! 234 202 235 END SUBROUTINE zdf_tke 203 236 … … 225 258 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 226 259 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient 227 REAL(wp) :: z bbrau, zesh2! temporary scalars228 REAL(wp) :: zfact1 , zfact2, zfact3! - -260 REAL(wp) :: zesh2 ! temporary scalars 261 REAL(wp) :: zfact1 ! - - 229 262 REAL(wp) :: ztx2 , zty2 , zcof ! - - 230 263 REAL(wp) :: ztau , zdif ! - - … … 233 266 !!bfr REAL(wp) :: zebot ! - - 234 267 INTEGER , POINTER, DIMENSION(:,: ) :: imlc 235 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc 268 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc, zbbrau,zfact2,zfact3 236 269 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 237 270 !!-------------------------------------------------------------------- … … 240 273 ! 241 274 CALL wrk_alloc( jpi,jpj, imlc ) ! integer 242 CALL wrk_alloc( jpi,jpj, zhlc ) 243 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 244 ! 245 zbbrau = rn_ebb / rau0 ! Local constant initialisation 246 zfact1 = -.5_wp * rdt 247 zfact2 = 1.5_wp * rdt * rn_ediss 248 zfact3 = 0.5_wp * rn_ediss 275 CALL wrk_alloc( jpi,jpj, zhlc ) 276 CALL wrk_alloc( jpi,jpj, zbbrau ) 277 CALL wrk_alloc( jpi,jpj, zfact2 ) 278 CALL wrk_alloc( jpi,jpj, zfact3 ) 279 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 280 ! 281 zbbrau = rn_ebb0 / rau0 ! Local constant initialisation 282 zfact1 = -.5_wp * rdt 283 zfact2 = 1.5_wp * rdt * rn_ediss0 284 zfact3 = 0.5_wp * rn_ediss0 249 285 ! 250 286 ! … … 261 297 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 262 298 DO ji = fs_2, fs_jpim1 ! vector opt. 263 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)264 END DO 265 END DO 266 299 en(ji,jj,1) = MAX( rn_emin0, zbbrau(ji,jj) * taum(ji,jj) ) * tmask(ji,jj,1) 300 END DO 301 END DO 302 267 303 !!bfr - start commented area 268 304 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 303 339 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 304 340 DO jk = jpkm1, 2, -1 305 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us 341 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us 306 342 DO ji = 1, jpi ! with us=0.016*wind(starting from jpk-1) 307 343 zus = zcof * taum(ji,jj) … … 311 347 END DO 312 348 ! ! finite LC depth 313 DO jj = 1, jpj 349 DO jj = 1, jpj 314 350 DO ji = 1, jpi 315 351 zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) … … 326 362 ! ! vertical velocity due to LC 327 363 zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) ) 328 zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )364 zwlc = zind * rn_lc0(ji,jj) * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 329 365 ! ! TKE Langmuir circulation source term 330 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( 1._wp - fr_i(ji,jj) ) 366 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( 1._wp - fr_i(ji,jj) )* ( zwlc * zwlc * zwlc ) / & 331 367 & zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 368 332 369 END DO 333 370 END DO … … 347 384 DO ji = 1, jpi 348 385 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 349 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 386 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 350 387 & / ( fse3uw_n(ji,jj,jk) & 351 388 & * fse3uw_b(ji,jj,jk) ) … … 376 413 ! ! shear prod. at w-point weightened by mask 377 414 zesh2 = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 378 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 415 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 379 416 ! 380 417 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 381 418 zd_lw(ji,jj,jk) = zzd_lw 382 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)419 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2(ji,jj) * dissl(ji,jj,jk) * tmask(ji,jj,jk) 383 420 ! 384 421 ! ! right hand side in en 385 422 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 386 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) &423 & + zfact3(ji,jj) * dissl(ji,jj,jk) * en (ji,jj,jk) ) & 387 424 & * wmask(ji,jj,jk) 388 425 END DO … … 433 470 END DO 434 471 435 ! ! Save TKE prior to nn_etau addition 436 e_niw(:,:,:) = en(:,:,:) 437 ! 472 ! ! Save TKE prior to nn_etau addition 473 e_niw(:,:,:) = en(:,:,:) 474 ! 438 475 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 439 476 ! ! TKE due to surface and internal wave breaking … … 455 492 DO jj = 2, jpjm1 456 493 DO ji = fs_2, fs_jpim1 ! vector opt. 457 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) &494 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 458 495 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 459 496 END DO … … 464 501 DO ji = fs_2, fs_jpim1 ! vector opt. 465 502 jk = nmln(ji,jj) 466 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) &503 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 467 504 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 468 505 END DO … … 477 514 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 478 515 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 479 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 480 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 516 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 517 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 481 518 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 482 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) &519 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau(ji,jj) * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 483 520 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 484 521 END DO … … 486 523 END DO 487 524 ELSEIF( nn_etau == 4 ) THEN !* column integral independant of htau (rn_efr must be scaled up) 488 IF( nn_htau == 2 ) THEN ! efr dependant on time-varying htau 525 IF( nn_htau == 2 ) THEN ! efr dependant on time-varying htau 489 526 DO jj = 2, jpjm1 490 527 DO ji = fs_2, fs_jpim1 ! vector opt. … … 504 541 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 505 542 ! 506 DO jk = 2, jpkm1 ! TKE budget: near-inertial waves term 507 DO jj = 2, jpjm1 508 DO ji = fs_2, fs_jpim1 ! vector opt. 509 e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk) 510 END DO 511 END DO 512 END DO 513 ! 514 CALL lbc_lnk( e_niw, 'W', 1. ) 543 DO jk = 2, jpkm1 ! TKE budget: near-inertial waves term 544 DO jj = 2, jpjm1 545 DO ji = fs_2, fs_jpim1 ! vector opt. 546 e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk) 547 END DO 548 END DO 549 END DO 550 ! 551 CALL lbc_lnk( e_niw, 'W', 1. ) 515 552 ! 516 553 CALL wrk_dealloc( jpi,jpj, imlc ) ! integer 517 CALL wrk_dealloc( jpi,jpj, zhlc ) 518 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 554 CALL wrk_dealloc( jpi,jpj, zhlc ) 555 CALL wrk_dealloc( jpi,jpj, zbbrau ) 556 CALL wrk_dealloc( jpi,jpj, zfact2 ) 557 CALL wrk_dealloc( jpi,jpj, zfact3 ) 558 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 519 559 ! 520 560 IF( nn_timing == 1 ) CALL timing_stop('tke_tke') … … 529 569 !! ** Purpose : Compute the vertical eddy viscosity and diffusivity 530 570 !! 531 !! ** Method : At this stage, en, the now TKE, is known (computed in 532 !! the tke_tke routine). First, the now mixing lenth is 571 !! ** Method : At this stage, en, the now TKE, is known (computed in 572 !! the tke_tke routine). First, the now mixing lenth is 533 573 !! computed from en and the strafification (N^2), then the mixings 534 574 !! coefficients are computed. 535 575 !! - Mixing length : a first evaluation of the mixing lengh 536 576 !! scales is: 537 !! mxl = sqrt(2*en) / N 577 !! mxl = sqrt(2*en) / N 538 578 !! where N is the brunt-vaisala frequency, with a minimum value set 539 579 !! to rmxl_min (rn_mxl0) in the interior (surface) ocean. 540 !! The mixing and dissipative length scale are bound as follow : 580 !! The mixing and dissipative length scale are bound as follow : 541 581 !! nn_mxl=0 : mxl bounded by the distance to surface and bottom. 542 582 !! zmxld = zmxlm = mxl 543 583 !! nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl 544 !! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 584 !! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 545 585 !! less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl 546 586 !! nn_mxl=3 : mxl is bounded from the surface to the bottom usings 547 !! |d/dz(xml)|<1 to obtain lup, and from the bottom to 548 !! the surface to obtain ldown. the resulting length 587 !! |d/dz(xml)|<1 to obtain lup, and from the bottom to 588 !! the surface to obtain ldown. the resulting length 549 589 !! scales are: 550 !! zmxld = sqrt( lup * ldown ) 590 !! zmxld = sqrt( lup * ldown ) 551 591 !! zmxlm = min ( lup , ldown ) 552 592 !! - Vertical eddy viscosity and diffusivity: 553 593 !! avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 554 !! avt = max( avmb, pdlr * avm ) 594 !! avt = max( avmb, pdlr * avm ) 555 595 !! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise. 556 596 !! … … 567 607 IF( nn_timing == 1 ) CALL timing_start('tke_avn') 568 608 569 CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 609 CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 570 610 571 611 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 576 616 ! 577 617 ! initialisation of interior minimum value (avoid a 2d loop with mikt) 578 zmxlm(:,:,:) = rmxl_min 618 zmxlm(:,:,:) = rmxl_min 579 619 zmxld(:,:,:) = rmxl_min 580 620 ! … … 586 626 END DO 587 627 END DO 588 ELSE 628 ELSE 589 629 zmxlm(:,:,1) = rn_mxl0 590 630 ENDIF … … 604 644 ! !* Physical limits for the mixing length 605 645 ! 606 zmxld(:,:,1 ) = zmxlm(:,:,1) ! surface set to the minimum value 646 zmxld(:,:,1 ) = zmxlm(:,:,1) ! surface set to the minimum value 607 647 zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value 608 648 ! … … 698 738 DO ji = fs_2, fs_jpim1 ! vector opt. 699 739 zsqen = SQRT( en(ji,jj,jk) ) 700 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen740 zav = rn_ediff0(ji,jj) * zmxlm(ji,jj,jk) * zsqen 701 741 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 702 742 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) … … 704 744 END DO 705 745 END DO 746 #if defined key_traldf_c2d || key_traldf_c3d 747 IF( ln_stopack) THEN 748 IF(nn_spp_avt > 0 ) CALL spp_gen( 1, avt(:,:,jk), nn_spp_avt, rn_avt_sd, jk_spp_avt, jk) 749 IF(nn_spp_avm > 0 ) CALL spp_gen( 1, avm(:,:,jk), nn_spp_avm, rn_avm_sd, jk_spp_avm, jk) 750 ENDIF 751 #else 752 IF ( ln_stopack ) & 753 & CALL ctl_stop( 'tke_avn: parameter perturbation will only work with '// & 754 'key_traldf_c2d or key_traldf_c3d') 755 #endif 706 756 END DO 707 757 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) … … 749 799 ENDIF 750 800 ! 751 CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 801 CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 752 802 ! 753 803 IF( nn_timing == 1 ) CALL timing_stop('tke_avn') … … 759 809 !!---------------------------------------------------------------------- 760 810 !! *** ROUTINE zdf_tke_init *** 761 !! 762 !! ** Purpose : Initialization of the vertical eddy diffivity and 811 !! 812 !! ** Purpose : Initialization of the vertical eddy diffivity and 763 813 !! viscosity when using a tke turbulent closure scheme 764 814 !! … … 776 826 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 777 827 & rn_mxl0 , nn_pdl , ln_lc , rn_lc , & 778 & nn_etau , nn_htau , rn_efr , rn_c 828 & nn_etau , nn_htau , rn_efr , rn_c 779 829 !!---------------------------------------------------------------------- 780 830 … … 843 893 rn_mxl0 = rmxl_min 844 894 ENDIF 845 895 896 ALLOCATE( rn_lc0 (jpi,jpj) ) ; rn_lc0 = rn_lc 897 ALLOCATE( rn_ediff0(jpi,jpj) ) ; rn_ediff0 = rn_ediff 898 ALLOCATE( rn_ediss0(jpi,jpj) ) ; rn_ediss0 = rn_ediss 899 ALLOCATE( rn_ebb0 (jpi,jpj) ) ; rn_ebb0 = rn_ebb 900 ALLOCATE( rn_efr0 (jpi,jpj) ) ; rn_efr0 = rn_efr 901 846 902 IF( nn_etau == 2 ) THEN 847 903 ierr = zdf_mxl_alloc() … … 855 911 856 912 ! !* depth of penetration of surface tke 857 IF( nn_etau /= 0 ) THEN 913 IF( nn_etau /= 0 ) THEN 858 914 htau(:,:) = 0._wp 859 915 SELECT CASE( nn_htau ) ! Choice of the depth of penetration … … 861 917 htau(:,:) = 10._wp 862 918 CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees 863 htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) ) 919 htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) ) 864 920 CASE( 2 ) ! fraction of depth-integrated TKE within mixed-layer 865 921 rhtau = -1._wp / LOG( 1._wp - rn_c ) … … 904 960 END DO 905 961 dissl(:,:,:) = 1.e-12_wp 906 ! 962 ! 907 963 CALL tke_rst( nit000, 'READ' ) !* read or initialize all required files 908 964 ! … … 913 969 !!--------------------------------------------------------------------- 914 970 !! *** ROUTINE tke_rst *** 915 !! 971 !! 916 972 !! ** Purpose : Read or write TKE file (en) in restart file 917 973 !! 918 974 !! ** Method : use of IOM library 919 !! if the restart does not contain TKE, en is either 920 !! set to rn_emin or recomputed 975 !! if the restart does not contain TKE, en is either 976 !! set to rn_emin or recomputed 921 977 !!---------------------------------------------------------------------- 922 978 INTEGER , INTENT(in) :: kt ! ocean time-step … … 927 983 !!---------------------------------------------------------------------- 928 984 ! 929 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 985 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 930 986 ! ! --------------- 931 987 IF( ln_rstart ) THEN !* Read the restart file … … 956 1012 CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation) 957 1013 ! 958 avt_k (:,:,:) = avt (:,:,:)959 avm_k (:,:,:) = avm (:,:,:)960 avmu_k(:,:,:) = avmu(:,:,:)961 avmv_k(:,:,:) = avmv(:,:,:)962 !963 1014 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO 964 1015 ENDIF 965 1016 ELSE !* Start from rest 966 1017 en(:,:,:) = rn_emin * tmask(:,:,:) 967 DO jk = 1, jpk ! set the Kz to the background value968 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)969 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)970 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)971 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)972 END DO973 1018 ENDIF 974 1019 ! 1020 avt_k (:,:,:) = avt (:,:,:) 1021 avm_k (:,:,:) = avm (:,:,:) 1022 avmu_k(:,:,:) = avmu(:,:,:) 1023 avmv_k(:,:,:) = avmv(:,:,:) 1024 975 1025 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 976 1026 ! ! -------------------
Note: See TracChangeset
for help on using the changeset viewer.
