Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
- Property svn:keywords set to Id
r1756 r2528 25 25 !! 3.2 ! 2009-06 (G. Madec, S. Masson) TKE restart compatible with key_cpl 26 26 !! ! + cleaning of the parameters + bugs correction 27 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 27 28 !!---------------------------------------------------------------------- 28 29 #if defined key_zdftke || defined key_esopa … … 30 31 !! 'key_zdftke' TKE vertical physics 31 32 !!---------------------------------------------------------------------- 32 !! zdf_tke 33 !! tke_tke 34 !! tke_avn 35 !! tke_init: initialization, namelist read, and parameters control36 !! tke_rst 33 !! zdf_tke : update momentum and tracer Kz from a tke scheme 34 !! tke_tke : tke time stepping: update tke at now time step (en) 35 !! tke_avn : compute mixing length scale and deduce avm and avt 36 !! zdf_tke_init : initialization, namelist read, and parameters control 37 !! tke_rst : read/write tke restart in ocean restart file 37 38 !!---------------------------------------------------------------------- 38 USE oce ! ocean dynamics and active tracers39 USE dom_oce ! ocean space and time domain40 USE dom vvl ! ocean space and time domain : variable volume layer41 USE zdf_oce ! ocean vertical physics39 USE oce ! ocean: dynamics and active tracers variables 40 USE phycst ! physical constants 41 USE dom_oce ! domain: ocean 42 USE domvvl ! domain: variable volume layer 42 43 USE sbc_oce ! surface boundary condition: ocean 43 USE phycst ! physical constants44 USE zdfmxl ! mixed layer45 USE restart ! o nly for lrst_oce44 USE zdf_oce ! vertical physics: ocean variables 45 USE zdfmxl ! vertical physics: mixed layer 46 USE restart ! ocean restart 46 47 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 47 48 USE prtctl ! Print control 48 49 USE in_out_manager ! I/O manager 49 50 USE iom ! I/O manager library 50 USE zdfbfr ! bottom friction51 51 52 52 IMPLICIT NONE 53 53 PRIVATE 54 54 55 PUBLIC zdf_tke ! routine called in step module 56 PUBLIC tke_rst ! routine called in step module 55 PUBLIC zdf_tke ! routine called in step module 56 PUBLIC zdf_tke_init ! routine called in opa module 57 PUBLIC tke_rst ! routine called in step module 57 58 58 59 LOGICAL , PUBLIC, PARAMETER :: lk_zdftke = .TRUE. !: TKE vertical mixing flag 59 60 60 61 #if defined key_c1d 61 ! !!* 1D cfg only *62 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e_dis, e_mix 63 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e_pdl, e_ric 62 ! !!** 1D cfg only ** ('key_c1d') 63 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e_dis, e_mix !: dissipation and mixing turbulent lengh scales 64 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: e_pdl, e_ric !: prandl and local Richardson numbers 64 65 #endif 65 66 66 ! !!! ** Namelist namzdf_tke ** 67 LOGICAL :: ln_mxl0 = .FALSE. ! mixing length scale surface value as function of wind stress or not 68 INTEGER :: nn_mxl = 2 ! type of mixing length (=0/1/2/3) 69 REAL(wp) :: rn_lmin0 = 0.4_wp ! surface min value of mixing length [m] 70 REAL(wp) :: rn_lmin = 0.1_wp ! interior min value of mixing length [m] 71 INTEGER :: nn_pdl = 1 ! Prandtl number or not (ratio avt/avm) (=0/1) 72 REAL(wp) :: rn_ediff = 0.1_wp ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 73 REAL(wp) :: rn_ediss = 0.7_wp ! coefficient of the Kolmogoroff dissipation 74 REAL(wp) :: rn_ebb = 3.75_wp ! coefficient of the surface input of tke 75 REAL(wp) :: rn_emin = 0.7071e-6_wp ! minimum value of tke [m2/s2] 76 REAL(wp) :: rn_emin0 = 1.e-4_wp ! surface minimum value of tke [m2/s2] 77 REAL(wp) :: rn_bshear = 1.e-20 ! background shear (>0) 78 INTEGER :: nn_etau = 0 ! type of depth penetration of surface tke (=0/1/2) 79 INTEGER :: nn_htau = 0 ! type of tke profile of penetration (=0/1) 80 REAL(wp) :: rn_efr = 1.0_wp ! fraction of TKE surface value which penetrates in the ocean 81 REAL(wp) :: rn_addhft = 0.0_wp ! add offset applied to HF tau 82 REAL(wp) :: rn_sclhft = 1.0_wp ! scale factor applied to HF tau 83 LOGICAL :: ln_lc = .FALSE. ! Langmuir cells (LC) as a source term of TKE or not 84 REAL(wp) :: rn_lc = 0.15_wp ! coef to compute vertical velocity of Langmuir cells 85 86 REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 87 88 REAL(wp), DIMENSION(jpi,jpj) :: htau ! depth of tke penetration (nn_htau) 89 REAL(wp), DIMENSION(jpi,jpj,jpk) :: en ! now turbulent kinetic energy [m2/s2] 90 REAL(wp), DIMENSION(jpi,jpj,jpk) :: dissl ! now mixing lenght of dissipation 67 ! !!** Namelist namzdf_tke ** 68 LOGICAL :: ln_mxl0 = .FALSE. ! mixing length scale surface value as function of wind stress or not 69 INTEGER :: nn_mxl = 2 ! type of mixing length (=0/1/2/3) 70 REAL(wp) :: rn_mxl0 = 0.04_wp ! surface min value of mixing length (kappa*z_o=0.4*0.1 m) [m] 71 INTEGER :: nn_pdl = 1 ! Prandtl number or not (ratio avt/avm) (=0/1) 72 REAL(wp) :: rn_ediff = 0.1_wp ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 73 REAL(wp) :: rn_ediss = 0.7_wp ! coefficient of the Kolmogoroff dissipation 74 REAL(wp) :: rn_ebb = 3.75_wp ! coefficient of the surface input of tke 75 REAL(wp) :: rn_emin = 0.7071e-6_wp ! minimum value of tke [m2/s2] 76 REAL(wp) :: rn_emin0 = 1.e-4_wp ! surface minimum value of tke [m2/s2] 77 REAL(wp) :: rn_bshear = 1.e-20_wp ! background shear (>0) currently a numerical threshold (do not change it) 78 INTEGER :: nn_etau = 0 ! type of depth penetration of surface tke (=0/1/2/3) 79 INTEGER :: nn_htau = 0 ! type of tke profile of penetration (=0/1) 80 REAL(wp) :: rn_efr = 1.0_wp ! fraction of TKE surface value which penetrates in the ocean 81 LOGICAL :: ln_lc = .FALSE. ! Langmuir cells (LC) as a source term of TKE or not 82 REAL(wp) :: rn_lc = 0.15_wp ! coef to compute vertical velocity of Langmuir cells 83 84 REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 85 REAL(wp) :: rmxl_min ! minimum mixing length value (deduced from rn_ediff and rn_emin values) [m] 86 REAL(wp) :: rhftau_add = 1.e-3_wp ! add offset applied to HF part of taum (nn_etau=3) 87 REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3) 88 89 REAL(wp), DIMENSION(jpi,jpj,jpk), PUBLIC :: en ! now turbulent kinetic energy [m2/s2] 90 91 REAL(wp), DIMENSION(jpi,jpj) :: htau ! depth of tke penetration (nn_htau) 92 REAL(wp), DIMENSION(jpi,jpj,jpk) :: dissl ! now mixing lenght of dissipation 91 93 92 94 !! * Substitutions … … 94 96 # include "vectopt_loop_substitute.h90" 95 97 !!---------------------------------------------------------------------- 96 !! NEMO/OPA 3. 2 , LOCEAN-IPSL (2009)97 !! $Id : zdftke2.F90 1201 2008-09-24 13:24:21Z rblod$98 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)98 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 99 !! $Id$ 100 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 99 101 !!---------------------------------------------------------------------- 100 101 102 CONTAINS 102 103 … … 149 150 !!---------------------------------------------------------------------- 150 151 ! 151 IF( kt == nit000 ) CALL tke_init ! initialisation 152 ! 153 CALL tke_tke ! now tke (en) 154 ! 155 CALL tke_avn ! now avt, avm, avmu, avmv 152 CALL tke_tke ! now tke (en) 153 ! 154 CALL tke_avn ! now avt, avm, avmu, avmv 156 155 ! 157 156 END SUBROUTINE zdf_tke … … 165 164 !! 166 165 !! ** Method : - TKE surface boundary condition 167 !! - source term due to Langmuir cells ( ln_lc=T)166 !! - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T) 168 167 !! - source term due to shear (saved in avmu, avmv arrays) 169 168 !! - Now TKE : resolution of the TKE equation by inverting … … 180 179 !! 181 180 INTEGER :: ji, jj, jk ! dummy loop arguments 182 !! $INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! temporary scalar183 !! $INTEGER :: ikbt, ikbumm1, ikbvmm1 ! temporary scalar181 !!bfr INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! temporary scalar 182 !!bfr INTEGER :: ikbt, ikbumm1, ikbvmm1 ! temporary scalar 184 183 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 185 184 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient … … 190 189 REAL(wp) :: zus , zwlc , zind ! - - 191 190 REAL(wp) :: zzd_up, zzd_lw ! - - 192 !! $REAL(wp) :: zebot ! - -191 !!bfr REAL(wp) :: zebot ! - - 193 192 INTEGER , DIMENSION(jpi,jpj) :: imlc ! 2D workspace 194 193 REAL(wp), DIMENSION(jpi,jpj) :: zhlc ! - - … … 197 196 ! 198 197 zbbrau = rn_ebb / rau0 ! Local constant initialisation 199 zfact1 = -.5 * rdt200 zfact2 = 1.5 * rdt * rn_ediss201 zfact3 = 0.5 * rn_ediss198 zfact1 = -.5_wp * rdt 199 zfact2 = 1.5_wp * rdt * rn_ediss 200 zfact3 = 0.5_wp * rn_ediss 202 201 ! 203 202 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 209 208 END DO 210 209 END DO 211 ! 210 211 !!bfr - start commented area 212 212 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 213 213 ! ! Bottom boundary condition on tke … … 219 219 ! computational cost is justified 220 220 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 221 !222 221 ! en(bot) = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 223 222 !CDIR NOVERRCHK … … 225 224 !CDIR NOVERRCHK 226 225 !! DO ji = fs_2, fs_jpim1 ! vector opt. 227 !! ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) ) 228 !! ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 229 !! ikbum1 = MAX( ikbu-1, 1 ) 230 !! ikbvm1 = MAX( ikbv-1, 1 ) 231 !! ikbu = MIN( mbathy(ji,jj), mbathy(ji-1,jj) ) 232 !! ikbv = MIN( mbathy(ji,jj), mbathy(ji,jj-1) ) 233 !! ikbumm1 = MAX( ikbu-1, 1 ) 234 !! ikbvmm1 = MAX( ikbv-1, 1 ) 235 !! ikbt = MAX( mbathy(ji,jj), 1 ) 236 !! ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,ikbumm1) + & 237 !! bfrua(ji ,jj) * ub(ji ,jj ,ikbum1 ) 238 !! zty2 = bfrva(ji,jj ) * vb(ji ,jj ,ikbvm1) + & 239 !! bfrva(ji,jj-1) * vb(ji ,jj-1,ikbvmm1 ) 226 !! ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + & 227 !! bfrua(ji ,jj) * ub(ji ,jj,mbku(ji ,jj) ) 228 !! zty2 = bfrva(ji,jj ) * vb(ji,jj ,mbkv(ji,jj )) + & 229 !! bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) ) 240 230 !! zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. 241 !! en (ji,jj, ikbt) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)231 !! en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1) 242 232 !! END DO 243 233 !! END DO 244 ! 245 ! 246 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 247 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke 234 !!bfr - end commented area 235 ! 236 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 237 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke (Axell JGR 2002) 248 238 ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 249 239 ! 250 240 ! !* total energy produce by LC : cumulative sum over jk 251 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)241 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * fsdepw(:,:,1) * fse3w(:,:,1) 252 242 DO jk = 2, jpk 253 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)243 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * fsdepw(:,:,jk) * fse3w(:,:,jk) 254 244 END DO 255 245 ! !* finite Langmuir Circulation depth 256 246 zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag ) 257 imlc(:,:) = mb athy(:,:) ! Initialization to the number of w ocean point mbathy247 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 258 248 DO jk = jpkm1, 2, -1 259 249 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us … … 275 265 END DO 276 266 END DO 277 # if defined key_c1d278 hlc(:,:) = zhlc(:,:) * tmask(:,:,1) ! c1d configuration: save finite Langmuir Circulation depth279 # endif280 267 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 281 268 !CDIR NOVERRCHK … … 328 315 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) ) 329 316 ! ! shear prod. at w-point weightened by mask 330 zesh2 = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1. e0, umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &331 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1. e0, vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )317 zesh2 = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 318 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 332 319 ! 333 320 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 334 321 zd_lw(ji,jj,jk) = zzd_lw 335 zdiag(ji,jj,jk) = 1. e0- zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)322 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk) 336 323 ! 337 324 ! ! right hand side in en … … 384 371 ! ! TKE due to surface and internal wave breaking 385 372 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 386 IF( nn_etau == 1 ) THEN !* penetration throughout the water column373 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 387 374 DO jk = 2, jpkm1 388 375 DO jj = 2, jpjm1 389 376 DO ji = fs_2, fs_jpim1 ! vector opt. 390 377 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 391 & * ( 1. e0- fr_i(ji,jj) ) * tmask(ji,jj,jk)392 END DO 393 END DO 394 END DO 395 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) 378 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) 379 END DO 380 END DO 381 END DO 382 ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction) 396 383 DO jj = 2, jpjm1 397 384 DO ji = fs_2, fs_jpim1 ! vector opt. 398 385 jk = nmln(ji,jj) 399 386 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 400 & * ( 1. e0- fr_i(ji,jj) ) * tmask(ji,jj,jk)401 END DO 402 END DO 403 ELSEIF( nn_etau == 3 ) THEN !* penetration throughout the water column387 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) 388 END DO 389 END DO 390 ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability) 404 391 !CDIR NOVERRCHK 405 392 DO jk = 2, jpkm1 … … 410 397 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 411 398 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 412 ztau = 0.5 * SQRT( ztx2 * ztx2 + zty2 * zty2 )! module of the mean stress413 zdif = taum(ji,jj) - ztau ! mean of the module - moduleof the mean414 zdif = r n_sclhft * MAX( 0.e0, zdif + rn_addhft )! apply some modifications...399 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! module of the mean stress 400 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 401 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 415 402 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 416 & * ( 1. e0- fr_i(ji,jj) ) * tmask(ji,jj,jk)403 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) 417 404 END DO 418 405 END DO … … 439 426 !! mxl = sqrt(2*en) / N 440 427 !! where N is the brunt-vaisala frequency, with a minimum value set 441 !! to r n_lmin (rn_lmin0) in the interior (surface) ocean.428 !! to rmxl_min (rn_mxl0) in the interior (surface) ocean. 442 429 !! The mixing and dissipative length scale are bound as follow : 443 430 !! nn_mxl=0 : mxl bounded by the distance to surface and bottom. … … 479 466 ! !* Buoyancy length scale: l=sqrt(2*e/n**2) 480 467 ! 481 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 482 !!gm this should be useless 483 zmxlm(:,:,1) = 0.e0 484 !!gm end 485 zraug = vkarmn * 2.e5 / ( rau0 * grav ) 486 DO jj = 2, jpjm1 487 DO ji = fs_2, fs_jpim1 ! vector opt. 488 zmxlm(ji,jj,1) = MAX( rn_lmin0, zraug * taum(ji,jj) ) 489 END DO 490 END DO 491 ELSE ! surface set to the minimum value 492 zmxlm(:,:,1) = rn_lmin0 468 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 469 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 470 zmxlm(:,:,1) = MAX( rn_mxl0, zraug * taum(:,:) ) 471 ELSE ! surface set to the minimum value 472 zmxlm(:,:,1) = rn_mxl0 493 473 ENDIF 494 zmxlm(:,:,jpk) = rn_lmin ! bottomset to the interior minium value495 ! 496 !CDIR NOVERRCHK 497 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n**2)474 zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value 475 ! 476 !CDIR NOVERRCHK 477 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2) 498 478 !CDIR NOVERRCHK 499 479 DO jj = 2, jpjm1 … … 501 481 DO ji = fs_2, fs_jpim1 ! vector opt. 502 482 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 503 zmxlm(ji,jj,jk) = MAX( rn_lmin, SQRT( 2. * en(ji,jj,jk) / zrn2 ) ) 504 END DO 505 END DO 506 END DO 507 ! 508 !!gm CAUTION: to be added here the bottom boundary condition on zmxlm 483 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 484 END DO 485 END DO 486 END DO 487 ! 509 488 ! 510 489 ! !* Physical limits for the mixing length 511 490 ! 512 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the minimumvalue513 zmxld(:,:,jpk) = r n_lmin ! bottomset to the minimum value491 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the zmxlm value 492 zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value 514 493 ! 515 494 SELECT CASE ( nn_mxl ) … … 520 499 DO ji = fs_2, fs_jpim1 ! vector opt. 521 500 zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk), & 522 & fsdepw(ji,jj,mb athy(ji,jj)) - fsdepw(ji,jj,jk) )501 & fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 523 502 zmxlm(ji,jj,jk) = zemxl 524 503 zmxld(ji,jj,jk) = zemxl … … 625 604 DO jj = 2, jpjm1 626 605 DO ji = fs_2, fs_jpim1 ! vector opt. 627 zcoef = 0.5 / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk))606 zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 628 607 ! ! shear 629 608 zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) ) & … … 632 611 & + avmv(ji,jj ,jk) * ( vn(ji,jj ,jk-1) - vn(ji,jj ,jk) ) * ( vb(ji,jj ,jk-1) - vb(ji,jj ,jk) ) 633 612 ! ! local Richardson number 634 zri = MAX( rn2b(ji,jj,jk), 0. ) * avm(ji,jj,jk) / (zcoef * (zdku + zdkv + rn_bshear ))635 zpdlr = MAX( 0.1 , 0.2 / MAX( 0.2 , zri ) )613 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear ) 614 zpdlr = MAX( 0.1_wp, 0.2 / MAX( 0.2 , zri ) ) 636 615 !!gm and even better with the use of the "true" ri_crit=0.22222... (this change the results!) 637 !!gm zpdlr = MAX( 0.1 , ri_crit / MAX( ri_crit , zri ) )616 !!gm zpdlr = MAX( 0.1_wp, ri_crit / MAX( ri_crit , zri ) ) 638 617 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 639 618 # if defined key_c1d … … 656 635 657 636 658 SUBROUTINE tke_init637 SUBROUTINE zdf_tke_init 659 638 !!---------------------------------------------------------------------- 660 !! *** ROUTINE tke_init ***639 !! *** ROUTINE zdf_tke_init *** 661 640 !! 662 641 !! ** Purpose : Initialization of the vertical eddy diffivity and … … 672 651 INTEGER :: ji, jj, jk ! dummy loop indices 673 652 !! 674 NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , & 675 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 676 & rn_lmin , rn_lmin0 , nn_pdl , nn_etau , & 677 & nn_htau , rn_efr , rn_addhft, rn_sclhft, & 678 & ln_lc , rn_lc 653 NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , & 654 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 655 & rn_mxl0 , nn_pdl , ln_lc , rn_lc , & 656 & nn_etau , nn_htau , rn_efr 679 657 !!---------------------------------------------------------------------- 680 658 … … 682 660 READ ( numnam, namzdf_tke ) 683 661 684 ri_cri = 2. / ( 2. + rn_ediss / rn_ediff ) ! resulting critical Richardson number 662 ri_cri = 2._wp / ( 2._wp + rn_ediss / rn_ediff ) ! resulting critical Richardson number 663 rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) ) ! resulting minimum length to recover molecular viscosity 685 664 686 665 IF(lwp) THEN !* Control print 687 666 WRITE(numout,*) 688 WRITE(numout,*) 'zdf_tke : tke turbulent closure scheme - initialisation'689 WRITE(numout,*) '~~~~~~~~ '667 WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation' 668 WRITE(numout,*) '~~~~~~~~~~~~' 690 669 WRITE(numout,*) ' Namelist namzdf_tke : set tke mixing parameters' 691 670 WRITE(numout,*) ' coef. to compute avt rn_ediff = ', rn_ediff … … 698 677 WRITE(numout,*) ' prandl number flag nn_pdl = ', nn_pdl 699 678 WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0 700 WRITE(numout,*) ' surface mixing length minimum value rn_lmin0 = ', rn_lmin0 701 WRITE(numout,*) ' interior mixing length minimum value rn_lmin0 = ', rn_lmin0 679 WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0 680 WRITE(numout,*) ' flag to take into acc. Langmuir circ. ln_lc = ', ln_lc 681 WRITE(numout,*) ' coef to compute verticla velocity of LC rn_lc = ', rn_lc 702 682 WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau 703 683 WRITE(numout,*) ' flag for computation of exp. tke profile nn_htau = ', nn_htau 704 684 WRITE(numout,*) ' fraction of en which pene. the thermocline rn_efr = ', rn_efr 705 WRITE(numout,*) ' add offset applied to HF tau rn_addhft = ', rn_addhft706 WRITE(numout,*) ' scale factor applied to HF tau rn_sclhft = ', rn_sclhft707 WRITE(numout,*) ' flag to take into acc. Langmuir circ. ln_lc = ', ln_lc708 WRITE(numout,*) ' coef to compute verticla velocity of LC rn_lc = ', rn_lc709 685 WRITE(numout,*) 710 686 WRITE(numout,*) ' critical Richardson nb with your parameters ri_cri = ', ri_cri … … 712 688 713 689 ! !* Check of some namelist values 714 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 715 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 716 IF( nn_htau < 0 .OR. nn_htau > 2 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 717 IF( rn_lc < 0.15 .OR. rn_lc > 0.2 ) CALL ctl_stop( 'bad value: rn_lc must be between 0.15 and 0.2 ' ) 718 IF( rn_sclhft == 0. .AND. nn_etau == 3 ) CALL ctl_stop( 'force null HF tau to penetrate the thermocline...' ) 719 IF( .NOT. lhftau .AND. nn_etau == 3 ) CALL ctl_stop( 'bad flag: nn_etau == 3 must be used with HF tau' ) 720 690 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 691 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 692 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 693 #if ! key_coupled 694 IF( nn_etau == 3 ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 695 #endif 696 697 IF( ln_mxl0 ) THEN 698 IF(lwp) WRITE(numout,*) ' use a surface mixing length = F(stress) : set rn_mxl0 = rmxl_min' 699 rn_mxl0 = rmxl_min 700 ENDIF 701 721 702 IF( nn_etau == 2 ) CALL zdf_mxl( nit000 ) ! Initialization of nmln 722 703 … … 724 705 IF( nn_etau /= 0 ) THEN 725 706 SELECT CASE( nn_htau ) ! Choice of the depth of penetration 726 CASE( 0 ) ! constant depth penetration (here 10 meters) 727 htau(:,:) = 10.e0 728 CASE( 1 ) ! F(latitude) : 0.5m to 30m at high lat. 729 DO jj = 1, jpj 730 DO ji = 1, jpi 731 htau(ji,jj) = MAX( 0.5, 3./4. * MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) ) ) 732 END DO 733 END DO 734 CASE( 2 ) ! constant depth penetration (here 30 meters) 735 htau(:,:) = 30.e0 707 CASE( 0 ) ! constant depth penetration (here 10 meters) 708 htau(:,:) = 10._wp 709 CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees 710 htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) ) 736 711 END SELECT 737 712 ENDIF … … 744 719 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 745 720 END DO 746 dissl(:,:,:) = 1.e-12 721 dissl(:,:,:) = 1.e-12_wp 747 722 ! !* read or initialize all required files 748 723 CALL tke_rst( nit000, 'READ' ) 749 724 ! 750 END SUBROUTINE tke_init725 END SUBROUTINE zdf_tke_init 751 726 752 727 … … 825 800 LOGICAL, PUBLIC, PARAMETER :: lk_zdftke = .FALSE. !: TKE flag 826 801 CONTAINS 827 SUBROUTINE zdf_tke( kt ) ! Empty routine 802 SUBROUTINE zdf_tke_init ! Dummy routine 803 END SUBROUTINE zdf_tke_init 804 SUBROUTINE zdf_tke( kt ) ! Dummy routine 828 805 WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt 829 806 END SUBROUTINE zdf_tke
Note: See TracChangeset
for help on using the changeset viewer.
