- Timestamp:
- 2017-12-01T18:44:09+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r7813 r8882 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 30 !! - ! 2017-05 (G. Madec) add top/bottom friction as boundary condition (ln_drg) 29 31 !!---------------------------------------------------------------------- 30 #if defined key_zdftke 31 !!---------------------------------------------------------------------- 32 !! 'key_zdftke' TKE vertical physics 32 33 33 !!---------------------------------------------------------------------- 34 34 !! zdf_tke : update momentum and tracer Kz from a tke scheme … … 44 44 USE sbc_oce ! surface boundary condition: ocean 45 45 USE zdf_oce ! vertical physics: ocean variables 46 USE zdfdrg ! vertical physics: top/bottom drag coef. 46 47 USE zdfmxl ! vertical physics: mixed layer 47 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 48 USE prtctl ! Print control 48 ! 49 49 USE in_out_manager ! I/O manager 50 50 USE iom ! I/O manager library 51 51 USE lib_mpp ! MPP library 52 USE wrk_nemo ! work arrays 52 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 53 USE prtctl ! Print control 53 54 USE timing ! Timing 54 55 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 55 #if defined key_agrif56 USE agrif_opa_interp57 USE agrif_opa_update58 #endif59 56 60 57 IMPLICIT NONE … … 64 61 PUBLIC zdf_tke_init ! routine called in opa module 65 62 PUBLIC tke_rst ! routine called in step module 66 67 LOGICAL , PUBLIC, PARAMETER :: lk_zdftke = .TRUE. !: TKE vertical mixing flag68 63 69 64 ! !!** Namelist namzdf_tke ** … … 78 73 REAL(wp) :: rn_emin0 ! surface minimum value of tke [m2/s2] 79 74 REAL(wp) :: rn_bshear ! background shear (>0) currently a numerical threshold (do not change it) 75 LOGICAL :: ln_drg ! top/bottom friction forcing flag 80 76 INTEGER :: nn_etau ! type of depth penetration of surface tke (=0/1/2/3) 81 INTEGER :: nn_htau ! type of tke profile of penetration (=0/1)82 REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean77 INTEGER :: nn_htau ! type of tke profile of penetration (=0/1) 78 REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean 83 79 LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not 84 REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells80 REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells 85 81 86 82 REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values) … … 89 85 REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3) 90 86 91 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 92 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 93 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation 94 #if defined key_c1d 95 ! !!** 1D cfg only ** ('key_c1d') 96 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_dis, e_mix !: dissipation and mixing turbulent lengh scales 97 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers 98 #endif 87 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 88 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 89 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation 99 90 100 91 !! * Substitutions … … 111 102 !! *** FUNCTION zdf_tke_alloc *** 112 103 !!---------------------------------------------------------------------- 113 ALLOCATE( & 114 #if defined key_c1d 115 & e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) , & 116 & e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) , & 117 #endif 118 & htau (jpi,jpj) , dissl(jpi,jpj,jpk) , & 119 & apdlr(jpi,jpj,jpk) , STAT= zdf_tke_alloc ) 120 ! 104 ALLOCATE( htau(jpi,jpj) , dissl(jpi,jpj,jpk) , apdlr(jpi,jpj,jpk) , STAT= zdf_tke_alloc ) 105 ! 121 106 IF( lk_mpp ) CALL mpp_sum ( zdf_tke_alloc ) 122 107 IF( zdf_tke_alloc /= 0 ) CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays') … … 125 110 126 111 127 SUBROUTINE zdf_tke( kt )112 SUBROUTINE zdf_tke( kt, p_sh2, p_avm, p_avt ) 128 113 !!---------------------------------------------------------------------- 129 114 !! *** ROUTINE zdf_tke *** … … 162 147 !! 163 148 !! ** Action : compute en (now turbulent kinetic energy) 164 !! update avt, avm u, avmv(before vertical eddy coef.)149 !! update avt, avm (before vertical eddy coef.) 165 150 !! 166 151 !! References : Gaspar et al., JGR, 1990, … … 170 155 !! Bruchard OM 2002 171 156 !!---------------------------------------------------------------------- 172 INTEGER, INTENT(in) :: kt ! ocean time step 173 !!---------------------------------------------------------------------- 174 ! 175 #if defined key_agrif 176 ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2) 177 IF( .NOT.Agrif_Root() ) CALL Agrif_Tke 178 #endif 179 ! 180 IF( kt /= nit000 ) THEN ! restore before value to compute tke 181 avt (:,:,:) = avt_k (:,:,:) 182 avm (:,:,:) = avm_k (:,:,:) 183 avmu(:,:,:) = avmu_k(:,:,:) 184 avmv(:,:,:) = avmv_k(:,:,:) 185 ENDIF 186 ! 187 CALL tke_tke ! now tke (en) 188 ! 189 CALL tke_avn ! now avt, avm, avmu, avmv 190 ! 191 avt_k (:,:,:) = avt (:,:,:) 192 avm_k (:,:,:) = avm (:,:,:) 193 avmu_k(:,:,:) = avmu(:,:,:) 194 avmv_k(:,:,:) = avmv(:,:,:) 195 ! 196 #if defined key_agrif 197 ! Update child grid f => parent grid 198 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only 199 #endif 200 ! 157 INTEGER , INTENT(in ) :: kt ! ocean time step 158 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 159 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 160 !!---------------------------------------------------------------------- 161 ! 162 CALL tke_tke( gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt ) ! now tke (en) 163 ! 164 CALL tke_avn( gdepw_n, e3t_n, e3w_n, p_avm, p_avt ) ! now avt, avm, dissl 165 ! 201 166 END SUBROUTINE zdf_tke 202 167 203 168 204 SUBROUTINE tke_tke 169 SUBROUTINE tke_tke( pdepw, p_e3t, p_e3w, p_sh2, p_avm, p_avt ) 205 170 !!---------------------------------------------------------------------- 206 171 !! *** ROUTINE tke_tke *** … … 210 175 !! ** Method : - TKE surface boundary condition 211 176 !! - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T) 212 !! - source term due to shear ( saved in avmu, avmv arrays)177 !! - source term due to shear (= Kz dz[Ub] * dz[Un] ) 213 178 !! - Now TKE : resolution of the TKE equation by inverting 214 179 !! a tridiagonal linear system by a "methode de chasse" 215 180 !! - increase TKE due to surface and internal wave breaking 181 !! NB: when sea-ice is present, both LC parameterization 182 !! and TKE penetration are turned off when the ice fraction 183 !! is smaller than 0.25 216 184 !! 217 185 !! ** Action : - en : now turbulent kinetic energy) 218 !! - avmu, avmv : production of TKE by shear at u and v-points219 !! (= Kz dz[Ub] * dz[Un] )220 186 !! --------------------------------------------------------------------- 221 INTEGER :: ji, jj, jk ! dummy loop arguments 222 !!bfr INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! temporary scalar 223 !!bfr INTEGER :: ikbt, ikbumm1, ikbvmm1 ! temporary scalar 224 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 225 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient 226 REAL(wp) :: zbbrau, zesh2 ! temporary scalars 227 REAL(wp) :: zfact1, zfact2, zfact3 ! - - 228 REAL(wp) :: ztx2 , zty2 , zcof ! - - 229 REAL(wp) :: ztau , zdif ! - - 230 REAL(wp) :: zus , zwlc , zind ! - - 231 REAL(wp) :: zzd_up, zzd_lw ! - - 232 !!bfr REAL(wp) :: zebot ! - - 233 INTEGER , POINTER, DIMENSION(:,: ) :: imlc 234 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc 235 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw, z3du, z3dv 236 REAL(wp) :: zri ! local Richardson number 187 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdepw ! depth of w-points 188 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_e3t, p_e3w ! level thickness (t- & w-points) 189 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_sh2 ! shear production term 190 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 191 ! 192 INTEGER :: ji, jj, jk ! dummy loop arguments 193 REAL(wp) :: zetop, zebot, zmsku, zmskv ! local scalars 194 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 195 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient 196 REAL(wp) :: zbbrau, zri ! local scalars 197 REAL(wp) :: zfact1, zfact2, zfact3 ! - - 198 REAL(wp) :: ztx2 , zty2 , zcof ! - - 199 REAL(wp) :: ztau , zdif ! - - 200 REAL(wp) :: zus , zwlc , zind ! - - 201 REAL(wp) :: zzd_up, zzd_lw ! - - 202 INTEGER , DIMENSION(jpi,jpj) :: imlc 203 REAL(wp), DIMENSION(jpi,jpj) :: zhlc 204 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw 237 205 !!-------------------------------------------------------------------- 238 206 ! 239 IF( nn_timing == 1 ) CALL timing_start('tke_tke') 240 ! 241 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, z3du, z3dv ) 207 IF( ln_timing ) CALL timing_start('tke_tke') 244 208 ! 245 209 zbbrau = rn_ebb / rau0 ! Local constant initialisation … … 248 212 zfact3 = 0.5_wp * rn_ediss 249 213 ! 250 ! 251 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 252 ! ! Surface boundary condition on tke 253 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 214 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 215 ! ! Surface/top/bottom boundary condition on tke 216 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 217 218 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 219 DO ji = fs_2, fs_jpim1 ! vector opt. 220 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 221 END DO 222 END DO 254 223 IF ( ln_isfcav ) THEN 255 224 DO jj = 2, jpjm1 ! en(mikt(ji,jj)) = rn_emin … … 258 227 END DO 259 228 END DO 260 END IF 261 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 262 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 267 !!bfr - start commented area 229 ENDIF 230 ! 268 231 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 269 232 ! ! Bottom boundary condition on tke 270 233 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 271 234 ! 272 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 273 ! Tests to date have found the bottom boundary condition on tke to have very little effect. 274 ! The condition is coded here for completion but commented out until there is proof that the 275 ! computational cost is justified 276 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 277 ! en(bot) = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 278 !! DO jj = 2, jpjm1 279 !! DO ji = fs_2, fs_jpim1 ! vector opt. 280 !! ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + & 281 !! bfrua(ji ,jj) * ub(ji ,jj,mbku(ji ,jj) ) 282 !! zty2 = bfrva(ji,jj ) * vb(ji,jj ,mbkv(ji,jj )) + & 283 !! bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) ) 284 !! zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. 285 !! en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1) 286 !! END DO 287 !! END DO 288 !!bfr - end commented area 289 ! 290 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 291 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke (Axell JGR 2002) 235 ! en(bot) = (ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 236 ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75) 237 ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2 238 ! 239 IF( ln_drg ) THEN !== friction used as top/bottom boundary condition on TKE 240 ! 241 DO jj = 2, jpjm1 ! bottom friction 242 DO ji = fs_2, fs_jpim1 ! vector opt. 243 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 244 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 245 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 246 zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 & 247 & + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 ) 248 en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 249 END DO 250 END DO 251 IF( ln_isfcav ) THEN ! top friction 252 DO jj = 2, jpjm1 253 DO ji = fs_2, fs_jpim1 ! vector opt. 254 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 255 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 256 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 257 zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 & 258 & + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 ) 259 en(ji,jj,mikt(ji,jj)) = MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) ! masked at ocean surface 260 END DO 261 END DO 262 ENDIF 263 ! 264 ENDIF 265 ! 266 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 267 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke ! (Axell JGR 2002) 292 268 ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 293 269 ! 294 270 ! !* total energy produce by LC : cumulative sum over jk 295 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw_n(:,:,1) * e3w_n(:,:,1)271 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1) 296 272 DO jk = 2, jpk 297 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw_n(:,:,jk) * e3w_n(:,:,jk)273 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk) 298 274 END DO 299 275 ! !* finite Langmuir Circulation depth … … 311 287 DO jj = 1, jpj 312 288 DO ji = 1, jpi 313 zhlc(ji,jj) = gdepw_n(ji,jj,imlc(ji,jj))289 zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 314 290 END DO 315 291 END DO … … 320 296 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 321 297 ! ! vertical velocity due to LC 322 zind = 0.5 - SIGN( 0.5, gdepw_n(ji,jj,jk) - zhlc(ji,jj) )323 zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) )298 zind = 0.5 - SIGN( 0.5, pdepw(ji,jj,jk) - zhlc(ji,jj) ) 299 zwlc = zind * rn_lc * zus * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) 324 300 ! ! TKE Langmuir circulation source term 325 en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) &301 en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) & 326 302 & / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 327 303 END DO … … 338 314 ! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 339 315 ! 340 DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form) 341 DO jj = 1, jpjm1 342 DO ji = 1, fs_jpim1 ! vector opt. 343 z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji+1,jj,jk) ) & 344 & * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) & 345 & * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) & 346 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 347 z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji,jj+1,jk) ) & 348 & * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) & 349 & * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) & 350 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 351 END DO 352 END DO 353 END DO 354 ! 355 IF( nn_pdl == 1 ) THEN !* Prandtl number case: compute apdlr 356 ! Note that zesh2 is also computed in the next loop. 357 ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops 316 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 358 317 DO jk = 2, jpkm1 359 318 DO jj = 2, jpjm1 360 DO ji = fs_2, fs_jpim1 ! vector opt. 361 ! ! shear prod. at w-point weightened by mask 362 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 363 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 364 ! ! local Richardson number 365 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zesh2 + rn_bshear ) 319 DO ji = 2, jpim1 320 ! ! local Richardson number 321 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 322 ! ! inverse of Prandtl number 366 323 apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) ) 367 368 END DO 369 END DO 370 END DO 371 ! 324 END DO 325 END DO 326 END DO 372 327 ENDIF 373 328 ! … … 376 331 DO ji = fs_2, fs_jpim1 ! vector opt. 377 332 zcof = zfact1 * tmask(ji,jj,jk) 378 # if defined key_zdftmx_new 379 ! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability) 380 zzd_up = zcof * MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp ) & ! upper diagonal 381 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk ) ) 382 zzd_lw = zcof * MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp ) & ! lower diagonal 383 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 384 # else 385 zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal 386 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk ) ) 387 zzd_lw = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & ! lower diagonal 388 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 389 # endif 390 ! ! shear prod. at w-point weightened by mask 391 zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 392 & + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 333 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 334 ! ! eddy coefficient (ensure numerical stability) 335 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 336 & / ( p_e3t(ji,jj,jk ) * p_e3w(ji,jj,jk ) ) 337 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 338 & / ( p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk ) ) 393 339 ! 394 340 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 395 341 zd_lw(ji,jj,jk) = zzd_lw 396 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)342 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 397 343 ! 398 344 ! ! right hand side in en 399 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 400 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) & 401 & * wmask(ji,jj,jk) 345 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( p_sh2(ji,jj,jk) & ! shear 346 & - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 347 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation 348 & ) * wmask(ji,jj,jk) 402 349 END DO 403 350 END DO … … 442 389 END DO 443 390 END DO 444 391 ! 445 392 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 446 393 ! ! TKE due to surface and internal wave breaking 447 394 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 448 395 !!gm BUG : in the exp remove the depth of ssh !!! 396 !!gm i.e. use gde3w in argument (pdepw) 449 397 450 398 … … 453 401 DO jj = 2, jpjm1 454 402 DO ji = fs_2, fs_jpim1 ! vector opt. 455 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( - gdepw_n(ji,jj,jk) / htau(ji,jj) ) &456 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)403 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 404 & * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 457 405 END DO 458 406 END DO … … 462 410 DO ji = fs_2, fs_jpim1 ! vector opt. 463 411 jk = nmln(ji,jj) 464 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( - gdepw_n(ji,jj,jk) / htau(ji,jj) ) &465 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)412 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 413 & * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 466 414 END DO 467 415 END DO … … 475 423 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 476 424 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 477 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) ) & 478 & * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 479 END DO 480 END DO 481 END DO 482 ENDIF 483 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 484 ! 485 CALL wrk_dealloc( jpi,jpj, imlc ) ! integer 486 CALL wrk_dealloc( jpi,jpj, zhlc ) 487 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 488 ! 489 IF( nn_timing == 1 ) CALL timing_stop('tke_tke') 425 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 426 & * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 427 END DO 428 END DO 429 END DO 430 ENDIF 431 ! 432 IF( ln_timing ) CALL timing_stop('tke_tke') 490 433 ! 491 434 END SUBROUTINE tke_tke 492 435 493 436 494 SUBROUTINE tke_avn 437 SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt ) 495 438 !!---------------------------------------------------------------------- 496 439 !! *** ROUTINE tke_avn *** … … 524 467 !! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise. 525 468 !! 526 !! ** Action : - avt : now vertical eddy diffusivity (w-point) 527 !! - avmu, avmv : now vertical eddy viscosity at uw- and vw-points 528 !!---------------------------------------------------------------------- 469 !! ** Action : - avt, avm : now vertical eddy diffusivity and viscosity (w-point) 470 !!---------------------------------------------------------------------- 471 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdepw ! depth (w-points) 472 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_e3t, p_e3w ! level thickness (t- & w-points) 473 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 474 ! 529 475 INTEGER :: ji, jj, jk ! dummy loop indices 530 REAL(wp) :: zrn2, zraug, zcoef, zav 531 REAL(wp) :: zdku, zri, zsqen! - -532 REAL(wp) :: z dkv, zemxl, zemlm, zemlp! - -533 REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld476 REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars 477 REAL(wp) :: zdku, zdkv, zsqen ! - - 478 REAL(wp) :: zemxl, zemlm, zemlp ! - - 479 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmxlm, zmxld ! 3D workspace 534 480 !!-------------------------------------------------------------------- 535 481 ! 536 IF( nn_timing == 1 ) CALL timing_start('tke_avn') 537 538 CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 539 482 IF( ln_timing ) CALL timing_start('tke_avn') 483 ! 540 484 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 541 485 ! ! Mixing length … … 549 493 ! 550 494 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 495 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 551 496 DO jj = 2, jpjm1 552 497 DO ji = fs_2, fs_jpim1 553 zraug = vkarmn * 2.e5_wp / ( rau0 * grav )554 498 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 555 499 END DO … … 576 520 ! 577 521 !!gm Not sure of that coding for ISF.... 578 ! where wmask = 0 set zmxlm == e3w_n522 ! where wmask = 0 set zmxlm == p_e3w 579 523 CASE ( 0 ) ! bounded by the distance to surface and bottom 580 524 DO jk = 2, jpkm1 581 525 DO jj = 2, jpjm1 582 526 DO ji = fs_2, fs_jpim1 ! vector opt. 583 zemxl = MIN( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), &584 & gdepw_n(ji,jj,mbkt(ji,jj)+1) - gdepw_n(ji,jj,jk) )527 zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 528 & pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) 585 529 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 586 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))587 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))530 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 531 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 588 532 END DO 589 533 END DO … … 594 538 DO jj = 2, jpjm1 595 539 DO ji = fs_2, fs_jpim1 ! vector opt. 596 zemxl = MIN( e3w_n(ji,jj,jk), zmxlm(ji,jj,jk) )540 zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 597 541 zmxlm(ji,jj,jk) = zemxl 598 542 zmxld(ji,jj,jk) = zemxl … … 605 549 DO jj = 2, jpjm1 606 550 DO ji = fs_2, fs_jpim1 ! vector opt. 607 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )551 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 608 552 END DO 609 553 END DO … … 612 556 DO jj = 2, jpjm1 613 557 DO ji = fs_2, fs_jpim1 ! vector opt. 614 zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )558 zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 615 559 zmxlm(ji,jj,jk) = zemxl 616 560 zmxld(ji,jj,jk) = zemxl … … 623 567 DO jj = 2, jpjm1 624 568 DO ji = fs_2, fs_jpim1 ! vector opt. 625 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )569 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 626 570 END DO 627 571 END DO … … 630 574 DO jj = 2, jpjm1 631 575 DO ji = fs_2, fs_jpim1 ! vector opt. 632 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )576 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 633 577 END DO 634 578 END DO … … 647 591 END SELECT 648 592 ! 649 # if defined key_c1d 650 e_dis(:,:,:) = zmxld(:,:,:) ! c1d configuration : save mixing and dissipation turbulent length scales 651 e_mix(:,:,:) = zmxlm(:,:,:) 652 # endif 653 654 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 655 ! ! Vertical eddy viscosity and diffusivity (avmu, avmv, avt) 593 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 594 ! ! Vertical eddy viscosity and diffusivity (avm and avt) 656 595 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 657 596 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points … … 660 599 zsqen = SQRT( en(ji,jj,jk) ) 661 600 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 662 avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk)663 avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)601 p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 602 p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 664 603 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 665 604 END DO 666 605 END DO 667 606 END DO 668 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 669 ! 670 DO jk = 2, jpkm1 !* vertical eddy viscosity at wu- and wv-points 671 DO jj = 2, jpjm1 672 DO ji = fs_2, fs_jpim1 ! vector opt. 673 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 674 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 675 END DO 676 END DO 677 END DO 678 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 607 ! 679 608 ! 680 609 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt … … 682 611 DO jj = 2, jpjm1 683 612 DO ji = fs_2, fs_jpim1 ! vector opt. 684 avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 685 # if defined key_c1d 686 e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 687 !!gm bug NO zri here.... 688 !!gm remove the specific diag for c1d ! 689 e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk) ! c1d config. : save Ri 690 # endif 613 p_avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 691 614 END DO 692 615 END DO 693 616 END DO 694 617 ENDIF 695 CALL lbc_lnk( avt, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged) 696 618 ! 697 619 IF(ln_ctl) THEN 698 CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 699 CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke - u: ', mask1=umask, & 700 & tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk ) 701 ENDIF 702 ! 703 CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 704 ! 705 IF( nn_timing == 1 ) CALL timing_stop('tke_avn') 620 CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 621 CALL prt_ctl( tab3d_1=avm, clinfo1=' tke - m: ', ovlap=1, kdim=jpk ) 622 ENDIF 623 ! 624 IF( ln_timing ) CALL timing_stop('tke_avn') 706 625 ! 707 626 END SUBROUTINE tke_avn … … 727 646 NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , & 728 647 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 729 & rn_mxl0 , nn_pdl , ln_ lc, rn_lc , &648 & rn_mxl0 , nn_pdl , ln_drg , ln_lc , rn_lc , & 730 649 & nn_etau , nn_htau , rn_efr 731 650 !!---------------------------------------------------------------------- … … 741 660 ! 742 661 ri_cri = 2._wp / ( 2._wp + rn_ediss / rn_ediff ) ! resulting critical Richardson number 743 # if defined key_zdftmx_new744 ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used745 rn_emin = 1.e-10_wp746 rmxl_min = 1.e-03_wp747 IF(lwp) THEN ! Control print748 WRITE(numout,*)749 WRITE(numout,*) 'zdf_tke_init : New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 '750 WRITE(numout,*) '~~~~~~~~~~~~'751 ENDIF752 # else753 rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) ) ! resulting minimum length to recover molecular viscosity754 # endif755 662 ! 756 663 IF(lwp) THEN !* Control print … … 764 671 WRITE(numout,*) ' minimum value of tke rn_emin = ', rn_emin 765 672 WRITE(numout,*) ' surface minimum value of tke rn_emin0 = ', rn_emin0 673 WRITE(numout,*) ' prandl number flag nn_pdl = ', nn_pdl 766 674 WRITE(numout,*) ' background shear (>0) rn_bshear = ', rn_bshear 767 675 WRITE(numout,*) ' mixing length type nn_mxl = ', nn_mxl 768 WRITE(numout,*) ' prandl number flag nn_pdl = ', nn_pdl769 WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0770 WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0771 WRITE(numout,*) ' flag to take into acc. Langmuir circ.ln_lc = ', ln_lc772 WRITE(numout,*) ' coef to compute verticla velocity of LC rn_lc= ', rn_lc676 WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0 677 WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0 678 WRITE(numout,*) ' top/bottom friction forcing flag ln_drg = ', ln_drg 679 WRITE(numout,*) ' Langmuir cells parametrization ln_lc = ', ln_lc 680 WRITE(numout,*) ' coef to compute vertical velocity of LC rn_lc = ', rn_lc 773 681 WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau 774 WRITE(numout,*) ' flag for computation of exp. tke profilenn_htau = ', nn_htau775 WRITE(numout,*) ' fraction of en which pene. the thermoclinern_efr = ', rn_efr682 WRITE(numout,*) ' type of tke penetration profile nn_htau = ', nn_htau 683 WRITE(numout,*) ' fraction of TKE that penetrates rn_efr = ', rn_efr 776 684 WRITE(numout,*) 777 WRITE(numout,*) ' critical Richardson nb with your parameters ri_cri = ', ri_cri 685 IF( ln_drg ) THEN 686 WRITE(numout,*) ' Namelist namdrg_top/_bot: used values:' 687 WRITE(numout,*) ' top ocean cavity roughness (m) rn_z0(_top)= ', r_z0_top 688 WRITE(numout,*) ' Bottom seafloor roughness (m) rn_z0(_bot)= ', r_z0_bot 689 ENDIF 690 WRITE(numout,*) 691 WRITE(numout,*) 692 WRITE(numout,*) ' ==>> critical Richardson nb with your parameters ri_cri = ', ri_cri 693 WRITE(numout,*) 694 ENDIF 695 ! 696 IF( ln_zdfiwm ) THEN ! Internal wave-driven mixing 697 rn_emin = 1.e-10_wp ! specific values of rn_emin & rmxl_min are used 698 rmxl_min = 1.e-03_wp ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s) 699 IF(lwp) WRITE(numout,*) ' Internal wave-driven mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 ' 700 ELSE ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s) 701 rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) ) ! resulting minimum length to recover molecular viscosity 702 IF(lwp) WRITE(numout,*) ' minimum mixing length with your parameters rmxl_min = ', rmxl_min 778 703 ENDIF 779 704 ! … … 805 730 ! !* set vertical eddy coef. to the background value 806 731 DO jk = 1, jpk 807 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 808 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 809 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 810 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 732 avt(:,:,jk) = avtb(jk) * wmask(:,:,jk) 733 avm(:,:,jk) = avmb(jk) * wmask(:,:,jk) 811 734 END DO 812 735 dissl(:,:,:) = 1.e-12_wp … … 818 741 819 742 SUBROUTINE tke_rst( kt, cdrw ) 820 !!--------------------------------------------------------------------- 821 !! *** ROUTINE tke_rst *** 822 !! 823 !! ** Purpose : Read or write TKE file (en) in restart file 824 !! 825 !! ** Method : use of IOM library 826 !! if the restart does not contain TKE, en is either 827 !! set to rn_emin or recomputed 828 !!---------------------------------------------------------------------- 829 INTEGER , INTENT(in) :: kt ! ocean time-step 830 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 831 ! 832 INTEGER :: jit, jk ! dummy loop indices 833 INTEGER :: id1, id2, id3, id4, id5, id6 ! local integers 834 !!---------------------------------------------------------------------- 835 ! 836 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 837 ! ! --------------- 838 IF( ln_rstart ) THEN !* Read the restart file 839 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 840 id2 = iom_varid( numror, 'avt' , ldstop = .FALSE. ) 841 id3 = iom_varid( numror, 'avm' , ldstop = .FALSE. ) 842 id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. ) 843 id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. ) 844 id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 845 ! 846 IF( id1 > 0 ) THEN ! 'en' exists 847 CALL iom_get( numror, jpdom_autoglo, 'en', en ) 848 IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN ! all required arrays exist 849 CALL iom_get( numror, jpdom_autoglo, 'avt' , avt ) 850 CALL iom_get( numror, jpdom_autoglo, 'avm' , avm ) 851 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu ) 852 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv ) 853 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 854 ELSE ! one at least array is missing 855 CALL tke_avn ! compute avt, avm, avmu, avmv and dissl (approximation) 856 ENDIF 857 ELSE ! No TKE array found: initialisation 858 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop' 859 en (:,:,:) = rn_emin * tmask(:,:,:) 860 CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation) 861 ! 862 avt_k (:,:,:) = avt (:,:,:) 863 avm_k (:,:,:) = avm (:,:,:) 864 avmu_k(:,:,:) = avmu(:,:,:) 865 avmv_k(:,:,:) = avmv(:,:,:) 866 ! 867 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO 868 ENDIF 869 ELSE !* Start from rest 870 en(:,:,:) = rn_emin * tmask(:,:,:) 871 DO jk = 1, jpk ! set the Kz to the background value 872 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 873 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 874 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 875 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 876 END DO 877 ENDIF 878 ! 879 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 880 ! ! ------------------- 881 IF(lwp) WRITE(numout,*) '---- tke-rst ----' 882 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 883 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 884 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 885 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 886 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 887 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 888 ! 889 ENDIF 890 ! 743 !!--------------------------------------------------------------------- 744 !! *** ROUTINE tke_rst *** 745 !! 746 !! ** Purpose : Read or write TKE file (en) in restart file 747 !! 748 !! ** Method : use of IOM library 749 !! if the restart does not contain TKE, en is either 750 !! set to rn_emin or recomputed 751 !!---------------------------------------------------------------------- 752 INTEGER , INTENT(in) :: kt ! ocean time-step 753 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 754 ! 755 INTEGER :: jit, jk ! dummy loop indices 756 INTEGER :: id1, id2, id3, id4 ! local integers 757 !!---------------------------------------------------------------------- 758 ! 759 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 760 ! ! --------------- 761 IF( ln_rstart ) THEN !* Read the restart file 762 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 763 id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 764 id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 765 id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 766 ! 767 IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! fields exist 768 CALL iom_get( numror, jpdom_autoglo, 'en', en ) 769 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 770 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 771 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 772 ELSE ! start TKE from rest 773 IF(lwp) WRITE(numout,*) ' ==>> previous run without TKE scheme, set en to background values' 774 en(:,:,:) = rn_emin * wmask(:,:,:) 775 ! avt_k, avm_k already set to the background value in zdf_phy_init 776 ENDIF 777 ELSE !* Start from rest 778 IF(lwp) WRITE(numout,*) ' ==>> start from rest: set en to the background value' 779 en(:,:,:) = rn_emin * wmask(:,:,:) 780 ! avt_k, avm_k already set to the background value in zdf_phy_init 781 ENDIF 782 ! 783 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 784 ! ! ------------------- 785 IF(lwp) WRITE(numout,*) '---- tke-rst ----' 786 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 787 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 788 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 789 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 790 ! 791 ENDIF 792 ! 891 793 END SUBROUTINE tke_rst 892 893 #else894 !!----------------------------------------------------------------------895 !! Dummy module : NO TKE scheme896 !!----------------------------------------------------------------------897 LOGICAL, PUBLIC, PARAMETER :: lk_zdftke = .FALSE. !: TKE flag898 CONTAINS899 SUBROUTINE zdf_tke_init ! Dummy routine900 END SUBROUTINE zdf_tke_init901 SUBROUTINE zdf_tke( kt ) ! Dummy routine902 WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt903 END SUBROUTINE zdf_tke904 SUBROUTINE tke_rst( kt, cdrw )905 CHARACTER(len=*) :: cdrw906 WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr907 END SUBROUTINE tke_rst908 #endif909 794 910 795 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
