Changeset 14072 for NEMO/trunk/src/OCE/ZDF/zdftke.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/ZDF/zdftke.F90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ZDF/zdftke.F90
r14057 r14072 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 28 28 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 29 !! 4.0 ! 2017-04 (G. Madec) remove CPP ddm key & avm at t-point only 29 !! 4.0 ! 2017-04 (G. Madec) remove CPP ddm key & avm at t-point only 30 30 !! - ! 2017-05 (G. Madec) add top/bottom friction as boundary condition 31 31 !! 4.2 ! 2020-12 (G. Madec, E. Clementi) add wave coupling … … 59 59 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 60 60 USE prtctl ! Print control 61 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 61 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 62 62 USE sbcwave ! Surface boundary waves 63 63 … … 78 78 INTEGER :: nn_pdl ! Prandtl number or not (ratio avt/avm) (=0/1) 79 79 REAL(wp) :: rn_ediff ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 80 REAL(wp) :: rn_ediss ! coefficient of the Kolmogoroff dissipation 80 REAL(wp) :: rn_ediss ! coefficient of the Kolmogoroff dissipation 81 81 REAL(wp) :: rn_ebb ! coefficient of the surface input of tke 82 82 REAL(wp) :: rn_emin ! minimum value of tke [m2/s2] … … 90 90 LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not 91 91 REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells 92 INTEGER :: nn_eice ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3) 92 INTEGER :: nn_eice ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3) 93 93 94 94 REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values) … … 139 139 !! surface: en = max( rn_emin0, rn_ebb * taum ) 140 140 !! bottom : en = rn_emin 141 !! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 142 !! 143 !! The now Turbulent kinetic energy is computed using the following 141 !! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 142 !! 143 !! The now Turbulent kinetic energy is computed using the following 144 144 !! time stepping: implicit for vertical diffusion term, linearized semi 145 !! implicit for kolmogoroff dissipation term, and explicit forward for 146 !! both buoyancy and shear production terms. Therefore a tridiagonal 145 !! implicit for kolmogoroff dissipation term, and explicit forward for 146 !! both buoyancy and shear production terms. Therefore a tridiagonal 147 147 !! linear system is solved. Note that buoyancy and shear terms are 148 148 !! discretized in a energy conserving form (Bruchard 2002). … … 152 152 !! 153 153 !! The now vertical eddy vicosity and diffusivity coefficients are 154 !! given by: 154 !! given by: 155 155 !! avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 156 !! avt = max( avmb, pdl * avm ) 156 !! avt = max( avmb, pdl * avm ) 157 157 !! eav = max( avmb, avm ) 158 158 !! where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and 159 !! given by an empirical funtion of the localRichardson number if nn_pdl=1 159 !! given by an empirical funtion of the localRichardson number if nn_pdl=1 160 160 !! 161 161 !! ** Action : compute en (now turbulent kinetic energy) … … 193 193 !! a tridiagonal linear system by a "methode de chasse" 194 194 !! - increase TKE due to surface and internal wave breaking 195 !! NB: when sea-ice is present, both LC parameterization 196 !! and TKE penetration are turned off when the ice fraction 197 !! is smaller than 0.25 195 !! NB: when sea-ice is present, both LC parameterization 196 !! and TKE penetration are turned off when the ice fraction 197 !! is smaller than 0.25 198 198 !! 199 199 !! ** Action : - en : now turbulent kinetic energy) … … 223 223 zbbrau = rn_ebb / rho0 ! Local constant initialisation 224 224 zbbirau = 3.75_wp / rho0 225 zfact1 = -.5_wp * rn_Dt 225 zfact1 = -.5_wp * rn_Dt 226 226 zfact2 = 1.5_wp * rn_Dt * rn_ediss 227 227 zfact3 = 0.5_wp * rn_ediss … … 244 244 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) 245 245 zdiag(ji,jj,1) = 1._wp/en(ji,jj,1) 246 zd_lw(ji,jj,1) = 1._wp 246 zd_lw(ji,jj,1) = 1._wp 247 247 zd_up(ji,jj,1) = 0._wp 248 248 END_2D … … 345 345 END_2D 346 346 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !* TKE Langmuir circulation source term added to en 347 IF ( zus3(ji,jj) /= 0._wp ) THEN 347 IF ( zus3(ji,jj) /= 0._wp ) THEN 348 348 IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 349 349 ! ! vertical velocity due to LC … … 376 376 END_3D 377 377 ENDIF 378 ! 378 ! 379 379 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !* Matrix and right hand side in en 380 380 zcof = zfact1 * tmask(ji,jj,jk) … … 403 403 ! 404 404 IF ( cpl_phioc .and. ln_phioc ) THEN 405 SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves 405 SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves 406 406 407 407 CASE ( 0 ) ! Dirichlet BC … … 456 456 ! 457 457 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 458 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 458 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 459 459 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 460 460 & * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) … … 470 470 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 471 471 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 472 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 473 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 472 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 473 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 474 474 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 475 475 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & … … 487 487 !! ** Purpose : Compute the vertical eddy viscosity and diffusivity 488 488 !! 489 !! ** Method : At this stage, en, the now TKE, is known (computed in 490 !! the tke_tke routine). First, the now mixing lenth is 489 !! ** Method : At this stage, en, the now TKE, is known (computed in 490 !! the tke_tke routine). First, the now mixing lenth is 491 491 !! computed from en and the strafification (N^2), then the mixings 492 492 !! coefficients are computed. 493 493 !! - Mixing length : a first evaluation of the mixing lengh 494 494 !! scales is: 495 !! mxl = sqrt(2*en) / N 495 !! mxl = sqrt(2*en) / N 496 496 !! where N is the brunt-vaisala frequency, with a minimum value set 497 497 !! to rmxl_min (rn_mxl0) in the interior (surface) ocean. 498 !! The mixing and dissipative length scale are bound as follow : 498 !! The mixing and dissipative length scale are bound as follow : 499 499 !! nn_mxl=0 : mxl bounded by the distance to surface and bottom. 500 500 !! zmxld = zmxlm = mxl 501 501 !! nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl 502 !! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 502 !! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 503 503 !! less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl 504 504 !! nn_mxl=3 : mxl is bounded from the surface to the bottom usings 505 !! |d/dz(xml)|<1 to obtain lup, and from the bottom to 506 !! the surface to obtain ldown. the resulting length 505 !! |d/dz(xml)|<1 to obtain lup, and from the bottom to 506 !! the surface to obtain ldown. the resulting length 507 507 !! scales are: 508 !! zmxld = sqrt( lup * ldown ) 508 !! zmxld = sqrt( lup * ldown ) 509 509 !! zmxlm = min ( lup , ldown ) 510 510 !! - Vertical eddy viscosity and diffusivity: 511 511 !! avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 512 !! avt = max( avmb, pdlr * avm ) 512 !! avt = max( avmb, pdlr * avm ) 513 513 !! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise. 514 514 !! … … 534 534 ! 535 535 ! initialisation of interior minimum value (avoid a 2d loop with mikt) 536 zmxlm(:,:,:) = rmxl_min 536 zmxlm(:,:,:) = rmxl_min 537 537 zmxld(:,:,:) = rmxl_min 538 538 ! … … 543 543 zmxlm(:,:,1)= zcoef * MAX ( 1.6 * hsw(:,:) , 0.02 ) ! surface mixing length = F(wave height) 544 544 ELSE 545 ! 545 ! 546 546 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g) 547 547 ! … … 603 603 ! !* Physical limits for the mixing length 604 604 ! 605 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the minimum value 605 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the minimum value 606 606 zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value 607 607 ! … … 686 686 !!---------------------------------------------------------------------- 687 687 !! *** ROUTINE zdf_tke_init *** 688 !! 689 !! ** Purpose : Initialization of the vertical eddy diffivity and 688 !! 689 !! ** Purpose : Initialization of the vertical eddy diffivity and 690 690 !! viscosity when using a tke turbulent closure scheme 691 691 !! … … 707 707 & rn_mxl0 , nn_mxlice, rn_mxlice, & 708 708 & nn_pdl , ln_lc , rn_lc , & 709 & nn_etau , nn_htau , rn_efr , nn_eice , & 709 & nn_etau , nn_htau , rn_efr , nn_eice , & 710 710 & nn_bc_surf, nn_bc_bot, ln_mxhsw 711 711 !!---------------------------------------------------------------------- … … 760 760 WRITE(numout,*) ' fraction of TKE that penetrates rn_efr = ', rn_efr 761 761 WRITE(numout,*) ' langmuir & surface wave breaking under ice nn_eice = ', nn_eice 762 SELECT CASE( nn_eice ) 762 SELECT CASE( nn_eice ) 763 763 CASE( 0 ) ; WRITE(numout,*) ' ==>>> no impact of ice cover on langmuir & surface wave breaking' 764 764 CASE( 1 ) ; WRITE(numout,*) ' ==>>> weigthed by 1-TANH( fr_i(:,:) * 10 )' … … 767 767 CASE DEFAULT 768 768 CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3') 769 END SELECT 769 END SELECT 770 770 WRITE(numout,*) 771 771 WRITE(numout,*) ' ==>>> critical Richardson nb with your parameters ri_cri = ', ri_cri … … 796 796 rn_mxl0 = rmxl_min 797 797 ENDIF 798 799 IF( nn_etau == 2 ) CALL zdf_mxl( nit000, Kmm ) ! Initialization of nmln 798 799 IF( nn_etau == 2 ) CALL zdf_mxl( nit000, Kmm ) ! Initialization of nmln 800 800 801 801 ! !* depth of penetration of surface tke 802 IF( nn_etau /= 0 ) THEN 802 IF( nn_etau /= 0 ) THEN 803 803 SELECT CASE( nn_htau ) ! Choice of the depth of penetration 804 804 CASE( 0 ) ! constant depth penetration (here 10 meters) 805 805 htau(:,:) = 10._wp 806 806 CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees 807 htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) ) 807 htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) ) 808 808 END SELECT 809 809 ENDIF 810 810 ! !* read or initialize all required files 811 CALL tke_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, dissl) 811 CALL tke_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, dissl) 812 812 ! 813 813 END SUBROUTINE zdf_tke_init … … 817 817 !!--------------------------------------------------------------------- 818 818 !! *** ROUTINE tke_rst *** 819 !! 819 !! 820 820 !! ** Purpose : Read or write TKE file (en) in restart file 821 821 !! 822 822 !! ** Method : use of IOM library 823 !! if the restart does not contain TKE, en is either 824 !! set to rn_emin or recomputed 823 !! if the restart does not contain TKE, en is either 824 !! set to rn_emin or recomputed 825 825 !!---------------------------------------------------------------------- 826 826 USE zdf_oce , ONLY : en, avt_k, avm_k ! ocean vertical physics … … 833 833 !!---------------------------------------------------------------------- 834 834 ! 835 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 835 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 836 836 ! ! --------------- 837 837 IF( ln_rstart ) THEN !* Read the restart file
Note: See TracChangeset
for help on using the changeset viewer.
