Changeset 14072 for NEMO/trunk/src/OCE/ZDF/zdfgls.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/ZDF/zdfgls.F90 (modified) (38 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ZDF/zdfgls.F90
r13970 r14072 2 2 !!====================================================================== 3 3 !! *** MODULE zdfgls *** 4 !! Ocean physics: vertical mixing coefficient computed from the gls 4 !! Ocean physics: vertical mixing coefficient computed from the gls 5 5 !! turbulent closure parameterization 6 6 !!====================================================================== 7 7 !! History : 3.0 ! 2009-09 (G. Reffray) Original code 8 8 !! 3.3 ! 2010-10 (C. Bricaud) Add in the reference 9 !! 4.0 ! 2017-04 (G. Madec) remove CPP keys & avm at t-point only 9 !! 4.0 ! 2017-04 (G. Madec) remove CPP keys & avm at t-point only 10 10 !! - ! 2017-05 (G. Madec) add top friction as boundary condition 11 11 !!---------------------------------------------------------------------- … … 16 16 !! gls_rst : read/write gls restart in ocean restart file 17 17 !!---------------------------------------------------------------------- 18 USE oce ! ocean dynamics and active tracers 18 USE oce ! ocean dynamics and active tracers 19 19 USE dom_oce ! ocean space and time domain 20 20 USE domvvl ! ocean space and time domain : variable volume layer … … 64 64 REAL(wp) :: rn_hsro ! Minimum surface roughness 65 65 REAL(wp) :: rn_hsri ! Ice ocean roughness 66 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 66 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 67 67 68 68 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters 69 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 70 REAL(wp) :: rl_sf = 0.2_wp ! 0 <rl_sf<vkarmn 69 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 70 REAL(wp) :: rl_sf = 0.2_wp ! 0 <rl_sf<vkarmn 71 71 REAL(wp) :: rghmin = -0.28_wp 72 72 REAL(wp) :: rgh0 = 0.0329_wp … … 75 75 REAL(wp) :: ra2 = 0.74_wp 76 76 REAL(wp) :: rb1 = 16.60_wp 77 REAL(wp) :: rb2 = 10.10_wp 78 REAL(wp) :: re2 = 1.33_wp 77 REAL(wp) :: rb2 = 10.10_wp 78 REAL(wp) :: re2 = 1.33_wp 79 79 REAL(wp) :: rl1 = 0.107_wp 80 80 REAL(wp) :: rl2 = 0.0032_wp … … 146 146 INTEGER :: itop, itopp1 ! - - 147 147 REAL(wp) :: zesh2, zsigpsi, zcoef, zex1 , zex2 ! local scalars 148 REAL(wp) :: ztx2, zty2, zup, zdown, zcof, zdir ! - - 148 REAL(wp) :: ztx2, zty2, zup, zdown, zcof, zdir ! - - 149 149 REAL(wp) :: zratio, zrn2, zflxb, sh , z_en ! - - 150 150 REAL(wp) :: prod, buoy, diss, zdiss, sm ! - - … … 153 153 REAL(wp), DIMENSION(jpi,jpj) :: zdep 154 154 REAL(wp), DIMENSION(jpi,jpj) :: zkar 155 REAL(wp), DIMENSION(jpi,jpj) :: zflxs ! Turbulence fluxed induced by internal waves 155 REAL(wp), DIMENSION(jpi,jpj) :: zflxs ! Turbulence fluxed induced by internal waves 156 156 REAL(wp), DIMENSION(jpi,jpj) :: zhsro ! Surface roughness (surface waves) 157 157 REAL(wp), DIMENSION(jpi,jpj) :: zice_fra ! Tapering of wave breaking under sea ice … … 167 167 ! Preliminary computing 168 168 169 ustar2_surf(:,:) = 0._wp ; psi(:,:,:) = 0._wp 169 ustar2_surf(:,:) = 0._wp ; psi(:,:,:) = 0._wp 170 170 ustar2_top (:,:) = 0._wp ; zwall_psi(:,:,:) = 0._wp 171 171 ustar2_bot (:,:) = 0._wp … … 177 177 CASE( 3 ) ; zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp ) 178 178 END SELECT 179 179 180 180 ! Compute surface, top and bottom friction at T-points 181 181 DO_2D( 0, 0, 0, 0 ) !== surface ocean friction … … 184 184 ! 185 185 !!gm Rq we may add here r_ke0(_top/_bot) ? ==>> think about that... 186 ! 186 ! 187 187 IF( .NOT.ln_drg_OFF ) THEN !== top/bottom friction (explicit before friction) 188 188 DO_2D( 0, 0, 0, 0 ) ! bottom friction (explicit before friction) … … 201 201 ENDIF 202 202 ENDIF 203 203 204 204 SELECT CASE ( nn_z0_met ) !== Set surface roughness length ==! 205 CASE ( 0 ) ! Constant roughness 205 CASE ( 0 ) ! Constant roughness 206 206 zhsro(:,:) = rn_hsro 207 207 CASE ( 1 ) ! Standard Charnock formula … … 271 271 IF( ln_sigpsi ) THEN 272 272 zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1. 273 zwall_psi(ji,jj,jk) = rsc_psi / & 273 zwall_psi(ji,jj,jk) = rsc_psi / & 274 274 & ( zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) ) 275 275 ELSE … … 286 286 & / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk,Kmm) ) 287 287 ! ! diagonal 288 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rn_Dt * zdiss * wmask(ji,jj,jk) 288 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rn_Dt * zdiss * wmask(ji,jj,jk) 289 289 ! ! right hand side in en 290 290 en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zesh2 * wmask(ji,jj,jk) … … 302 302 SELECT CASE ( nn_bc_surf ) 303 303 ! 304 CASE ( 0 ) ! Dirichlet boundary condition (set e at k=1 & 2) 304 CASE ( 0 ) ! Dirichlet boundary condition (set e at k=1 & 2) 305 305 ! First level 306 306 en (:,:,1) = MAX( rn_emin , rc02r * ustar2_surf(:,:) * (1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1)**r2_3 ) … … 308 308 zd_up(:,:,1) = 0._wp 309 309 zdiag(:,:,1) = 1._wp 310 ! 310 ! 311 311 ! One level below 312 312 en (:,:,2) = MAX( rc02r * ustar2_surf(:,:) * ( 1._wp + (1._wp-zice_fra(:,:))*rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm)) & 313 313 & / zhsro(:,:) )**(1.5_wp*ra_sf) )**(2._wp/3._wp) , rn_emin ) 314 zd_lw(:,:,2) = 0._wp 314 zd_lw(:,:,2) = 0._wp 315 315 zd_up(:,:,2) = 0._wp 316 316 zdiag(:,:,2) = 1._wp … … 345 345 SELECT CASE ( nn_bc_bot ) 346 346 ! 347 CASE ( 0 ) ! Dirichlet 347 CASE ( 0 ) ! Dirichlet 348 348 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = rn_lmin 349 349 ! ! Balance between the production and the dissipation terms … … 357 357 z_en = MAX( rc02r * ustar2_bot(ji,jj), rn_emin ) 358 358 ! 359 ! Dirichlet condition applied at: 360 ! Bottom level (ibot) & Just above it (ibotm1) 359 ! Dirichlet condition applied at: 360 ! Bottom level (ibot) & Just above it (ibotm1) 361 361 zd_lw(ji,jj,ibot) = 0._wp ; zd_lw(ji,jj,ibotm1) = 0._wp 362 362 zd_up(ji,jj,ibot) = 0._wp ; zd_up(ji,jj,ibotm1) = 0._wp … … 373 373 ! 374 374 !!gm TO BE VERIFIED !!! 375 ! Dirichlet condition applied at: 376 ! top level (itop) & Just below it (itopp1) 375 ! Dirichlet condition applied at: 376 ! top level (itop) & Just below it (itopp1) 377 377 zd_lw(ji,jj,itop) = 0._wp ; zd_lw(ji,jj,itopp1) = 0._wp 378 378 zd_up(ji,jj,itop) = 0._wp ; zd_up(ji,jj,itopp1) = 0._wp … … 383 383 ! 384 384 CASE ( 1 ) ! Neumman boundary condition 385 ! 385 ! 386 386 DO_2D( 0, 0, 0, 0 ) 387 387 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point … … 391 391 ! 392 392 ! Bottom level Dirichlet condition: 393 ! Bottom level (ibot) & Just above it (ibotm1) 393 ! Bottom level (ibot) & Just above it (ibotm1) 394 394 ! Dirichlet ! Neumann 395 395 zd_lw(ji,jj,ibot) = 0._wp ! ! Remove zd_up from zdiag … … 405 405 ! 406 406 ! Bottom level Dirichlet condition: 407 ! Bottom level (ibot) & Just above it (ibotm1) 407 ! Bottom level (ibot) & Just above it (ibotm1) 408 408 ! Dirichlet ! Neumann 409 409 zd_lw(ji,jj,itop) = 0._wp ! ! Remove zd_up from zdiag … … 427 427 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 428 428 END_3D 429 ! ! set the minimum value of tke 429 ! ! set the minimum value of tke 430 430 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 431 431 … … 455 455 CASE( 3 ) ! generic 456 456 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 457 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 457 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * hmxl_b(ji,jj,jk)**rnn 458 458 END_3D 459 459 ! … … 470 470 ! 471 471 ! psi / k 472 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 472 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 473 473 ! 474 474 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 zdir = 1 (stable) otherwise zdir = 0 (unstable) … … 490 490 zesh2 = zdir * ( prod + buoy ) + (1._wp - zdir ) * prod ! production term 491 491 zdiss = zdir * ( diss / psi(ji,jj,jk) ) + (1._wp - zdir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 492 ! 492 ! 493 493 ! building the matrix 494 494 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) … … 528 528 zd_up(:,:,2) = 0._wp 529 529 zdiag(:,:,2) = 1._wp 530 ! 530 ! 531 531 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 532 532 ! … … 564 564 SELECT CASE ( nn_bc_bot ) ! bottom boundary 565 565 ! 566 CASE ( 0 ) ! Dirichlet 566 CASE ( 0 ) ! Dirichlet 567 567 ! ! en(ibot) = u*^2 / Co2 and hmxl_n(ibot) = vkarmn * r_z0_bot 568 568 ! ! Balance between the production and the dissipation terms … … 585 585 ! 586 586 CASE ( 1 ) ! Neumman boundary condition 587 ! 587 ! 588 588 DO_2D( 0, 0, 0, 0 ) 589 589 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point … … 641 641 CASE( 2 ) ! k-w 642 642 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 643 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 643 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 644 644 END_3D 645 645 ! … … 660 660 eps (ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 661 661 hmxl_n(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 662 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 662 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 663 663 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 664 664 IF( ln_length_lim ) hmxl_n(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), hmxl_n(ji,jj,jk) ) … … 720 720 721 721 ! default value, in case jpk > mbkt(ji,jj)+1. Not needed but avoid a bug when looking for undefined values (-fpe0) 722 zstm(:,:,jpk) = 0. 722 zstm(:,:,jpk) = 0. 723 723 DO_2D( 0, 0, 0, 0 ) ! update bottom with good values 724 724 zstm(ji,jj,mbkt(ji,jj)+1) = zstm(ji,jj,mbkt(ji,jj)) … … 756 756 !!---------------------------------------------------------------------- 757 757 !! *** ROUTINE zdf_gls_init *** 758 !! 759 !! ** Purpose : Initialization of the vertical eddy diffivity and 758 !! 759 !! ** Purpose : Initialization of the vertical eddy diffivity and 760 760 !! viscosity computed using a GLS turbulent closure scheme 761 761 !! … … 983 983 ! 984 984 END SELECT 985 985 986 986 ! !* Set Schmidt number for psi diffusion in the wave breaking case 987 987 ! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009 988 988 ! ! or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001 989 989 IF( ln_sigpsi ) THEN 990 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 990 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 991 991 ! Verification: retrieve Burchard (2001) results by uncomenting the line below: 992 992 ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work … … 996 996 rsc_psi0 = rsc_psi 997 997 ENDIF 998 998 999 999 ! !* Shear free turbulence parameters 1000 1000 ! … … 1039 1039 rc04 = rc03 * rc0 1040 1040 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking 1041 rsbc_tke2 = rn_Dt * rn_crban / rl_sf ! Neumann + Wave breaking 1041 rsbc_tke2 = rn_Dt * rn_crban / rl_sf ! Neumann + Wave breaking 1042 1042 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1043 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1043 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1044 1044 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1045 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1045 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1046 1046 rsbc_psi1 = -0.5_wp * rn_Dt * rc0**(rpp-2._wp*rmm) / rsc_psi 1047 rsbc_psi2 = -0.5_wp * rn_Dt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1047 rsbc_psi2 = -0.5_wp * rn_Dt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1048 1048 ! 1049 1049 rfact_tke = -0.5_wp / rsc_tke * rn_Dt ! Cst used for the Diffusion term of tke … … 1054 1054 zwall(:,:,:) = 1._wp * tmask(:,:,:) 1055 1055 1056 ! !* read or initialize all required files 1056 ! !* read or initialize all required files 1057 1057 CALL gls_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, hmxl_n) 1058 1058 ! … … 1063 1063 !!--------------------------------------------------------------------- 1064 1064 !! *** ROUTINE gls_rst *** 1065 !! 1065 !! 1066 1066 !! ** Purpose : Read or write TKE file (en) in restart file 1067 1067 !! 1068 1068 !! ** Method : use of IOM library 1069 !! if the restart does not contain TKE, en is either 1069 !! if the restart does not contain TKE, en is either 1070 1070 !! set to rn_emin or recomputed (nn_igls/=0) 1071 1071 !!---------------------------------------------------------------------- … … 1081 1081 !!---------------------------------------------------------------------- 1082 1082 ! 1083 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 1083 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 1084 1084 ! ! --------------- 1085 1085 IF( ln_rstart ) THEN !* Read the restart file … … 1094 1094 CALL iom_get( numror, jpdom_auto, 'avm_k' , avm_k ) 1095 1095 CALL iom_get( numror, jpdom_auto, 'hmxl_n', hmxl_n ) 1096 ELSE 1096 ELSE 1097 1097 IF(lwp) WRITE(numout,*) 1098 1098 IF(lwp) WRITE(numout,*) ' ==>> previous run without GLS scheme, set en and hmxl_n to background values'
Note: See TracChangeset
for help on using the changeset viewer.
