- Timestamp:
- 2017-12-13T15:58:53+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r7813 r9019 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 … … 43 43 USE domvvl ! domain: variable volume layer 44 44 USE sbc_oce ! surface boundary condition: ocean 45 USE zdf _oce ! vertical physics: ocean variables45 USE zdfdrg ! vertical physics: top/bottom drag coef. 46 46 USE zdfmxl ! vertical physics: mixed layer 47 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 48 USE prtctl ! Print control 47 ! 49 48 USE in_out_manager ! I/O manager 50 49 USE iom ! I/O manager library 51 50 USE lib_mpp ! MPP library 52 USE wrk_nemo ! work arrays 51 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 52 USE prtctl ! Print control 53 53 USE timing ! Timing 54 54 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 55 60 56 IMPLICIT NONE … … 64 60 PUBLIC zdf_tke_init ! routine called in opa module 65 61 PUBLIC tke_rst ! routine called in step module 66 67 LOGICAL , PUBLIC, PARAMETER :: lk_zdftke = .TRUE. !: TKE vertical mixing flag68 62 69 63 ! !!** Namelist namzdf_tke ** … … 78 72 REAL(wp) :: rn_emin0 ! surface minimum value of tke [m2/s2] 79 73 REAL(wp) :: rn_bshear ! background shear (>0) currently a numerical threshold (do not change it) 74 LOGICAL :: ln_drg ! top/bottom friction forcing flag 80 75 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 ocean76 INTEGER :: nn_htau ! type of tke profile of penetration (=0/1) 77 REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean 83 78 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 cells79 REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells 85 80 86 81 REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values) … … 89 84 REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3) 90 85 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 86 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 87 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 88 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation 99 89 100 90 !! * Substitutions … … 111 101 !! *** FUNCTION zdf_tke_alloc *** 112 102 !!---------------------------------------------------------------------- 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 ! 103 ALLOCATE( htau(jpi,jpj) , dissl(jpi,jpj,jpk) , apdlr(jpi,jpj,jpk) , STAT= zdf_tke_alloc ) 104 ! 121 105 IF( lk_mpp ) CALL mpp_sum ( zdf_tke_alloc ) 122 106 IF( zdf_tke_alloc /= 0 ) CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays') … … 125 109 126 110 127 SUBROUTINE zdf_tke( kt )111 SUBROUTINE zdf_tke( kt, p_sh2, p_avm, p_avt ) 128 112 !!---------------------------------------------------------------------- 129 113 !! *** ROUTINE zdf_tke *** … … 162 146 !! 163 147 !! ** Action : compute en (now turbulent kinetic energy) 164 !! update avt, avm u, avmv(before vertical eddy coef.)148 !! update avt, avm (before vertical eddy coef.) 165 149 !! 166 150 !! References : Gaspar et al., JGR, 1990, … … 170 154 !! Bruchard OM 2002 171 155 !!---------------------------------------------------------------------- 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 ! 156 INTEGER , INTENT(in ) :: kt ! ocean time step 157 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 158 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 159 !!---------------------------------------------------------------------- 160 ! 161 CALL tke_tke( gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt ) ! now tke (en) 162 ! 163 CALL tke_avn( gdepw_n, e3t_n, e3w_n, p_avm, p_avt ) ! now avt, avm, dissl 164 ! 201 165 END SUBROUTINE zdf_tke 202 166 203 167 204 SUBROUTINE tke_tke 168 SUBROUTINE tke_tke( pdepw, p_e3t, p_e3w, p_sh2, p_avm, p_avt ) 205 169 !!---------------------------------------------------------------------- 206 170 !! *** ROUTINE tke_tke *** … … 210 174 !! ** Method : - TKE surface boundary condition 211 175 !! - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T) 212 !! - source term due to shear ( saved in avmu, avmv arrays)176 !! - source term due to shear (= Kz dz[Ub] * dz[Un] ) 213 177 !! - Now TKE : resolution of the TKE equation by inverting 214 178 !! a tridiagonal linear system by a "methode de chasse" 215 179 !! - increase TKE due to surface and internal wave breaking 180 !! NB: when sea-ice is present, both LC parameterization 181 !! and TKE penetration are turned off when the ice fraction 182 !! is smaller than 0.25 216 183 !! 217 184 !! ** 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 185 !! --------------------------------------------------------------------- 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 186 USE zdf_oce , ONLY : en ! ocean vertical physics 187 !! 188 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdepw ! depth of w-points 189 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_e3t, p_e3w ! level thickness (t- & w-points) 190 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_sh2 ! shear production term 191 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 192 ! 193 INTEGER :: ji, jj, jk ! dummy loop arguments 194 REAL(wp) :: zetop, zebot, zmsku, zmskv ! local scalars 195 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 196 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient 197 REAL(wp) :: zbbrau, zri ! local scalars 198 REAL(wp) :: zfact1, zfact2, zfact3 ! - - 199 REAL(wp) :: ztx2 , zty2 , zcof ! - - 200 REAL(wp) :: ztau , zdif ! - - 201 REAL(wp) :: zus , zwlc , zind ! - - 202 REAL(wp) :: zzd_up, zzd_lw ! - - 203 INTEGER , DIMENSION(jpi,jpj) :: imlc 204 REAL(wp), DIMENSION(jpi,jpj) :: zhlc 205 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpelc, zdiag, zd_up, zd_lw 237 206 !!-------------------------------------------------------------------- 238 207 ! 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 ) 208 IF( ln_timing ) CALL timing_start('tke_tke') 244 209 ! 245 210 zbbrau = rn_ebb / rau0 ! Local constant initialisation … … 248 213 zfact3 = 0.5_wp * rn_ediss 249 214 ! 250 ! 251 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 252 ! ! Surface boundary condition on tke 253 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 215 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 216 ! ! Surface/top/bottom boundary condition on tke 217 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 218 219 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 220 DO ji = fs_2, fs_jpim1 ! vector opt. 221 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 222 END DO 223 END DO 254 224 IF ( ln_isfcav ) THEN 255 225 DO jj = 2, jpjm1 ! en(mikt(ji,jj)) = rn_emin … … 258 228 END DO 259 229 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 230 ENDIF 231 ! 268 232 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 269 233 ! ! Bottom boundary condition on tke 270 234 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 271 235 ! 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) 236 ! en(bot) = (ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 237 ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75) 238 ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2 239 ! 240 IF( ln_drg ) THEN !== friction used as top/bottom boundary condition on TKE 241 ! 242 DO jj = 2, jpjm1 ! bottom friction 243 DO ji = fs_2, fs_jpim1 ! vector opt. 244 zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 245 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 246 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 247 zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 & 248 & + ( zmskv*( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 ) 249 en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 250 END DO 251 END DO 252 IF( ln_isfcav ) THEN ! top friction 253 DO jj = 2, jpjm1 254 DO ji = fs_2, fs_jpim1 ! vector opt. 255 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 256 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 257 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 258 zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku*( ub(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 & 259 & + ( zmskv*( vb(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 ) 260 en(ji,jj,mikt(ji,jj)) = MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) ! masked at ocean surface 261 END DO 262 END DO 263 ENDIF 264 ! 265 ENDIF 266 ! 267 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 268 IF( ln_lc ) THEN ! Langmuir circulation source term added to tke ! (Axell JGR 2002) 292 269 ! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 293 270 ! 294 271 ! !* total energy produce by LC : cumulative sum over jk 295 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw_n(:,:,1) * e3w_n(:,:,1)272 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1) 296 273 DO jk = 2, jpk 297 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw_n(:,:,jk) * e3w_n(:,:,jk)274 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk) 298 275 END DO 299 276 ! !* finite Langmuir Circulation depth … … 311 288 DO jj = 1, jpj 312 289 DO ji = 1, jpi 313 zhlc(ji,jj) = gdepw_n(ji,jj,imlc(ji,jj))290 zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj)) 314 291 END DO 315 292 END DO … … 320 297 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift 321 298 ! ! 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) )299 zind = 0.5 - SIGN( 0.5, pdepw(ji,jj,jk) - zhlc(ji,jj) ) 300 zwlc = zind * rn_lc * zus * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) 324 301 ! ! 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 ) &302 en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) & 326 303 & / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 327 304 END DO … … 338 315 ! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal 339 316 ! 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 317 IF( nn_pdl == 1 ) THEN !* Prandtl number = F( Ri ) 358 318 DO jk = 2, jpkm1 359 319 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 ) 320 DO ji = 2, jpim1 321 ! ! local Richardson number 322 zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear ) 323 ! ! inverse of Prandtl number 366 324 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 ! 325 END DO 326 END DO 327 END DO 372 328 ENDIF 373 329 ! … … 376 332 DO ji = fs_2, fs_jpim1 ! vector opt. 377 333 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) ) 334 ! ! A minimum of 2.e-5 m2/s is imposed on TKE vertical 335 ! ! eddy coefficient (ensure numerical stability) 336 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 337 & / ( p_e3t(ji,jj,jk ) * p_e3w(ji,jj,jk ) ) 338 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 339 & / ( p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk ) ) 393 340 ! 394 341 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 395 342 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)343 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk) 397 344 ! 398 345 ! ! 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) 346 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( p_sh2(ji,jj,jk) & ! shear 347 & - p_avt(ji,jj,jk) * rn2(ji,jj,jk) & ! stratification 348 & + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk) & ! dissipation 349 & ) * wmask(ji,jj,jk) 402 350 END DO 403 351 END DO … … 442 390 END DO 443 391 END DO 444 392 ! 445 393 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 446 394 ! ! TKE due to surface and internal wave breaking 447 395 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 448 396 !!gm BUG : in the exp remove the depth of ssh !!! 397 !!gm i.e. use gde3w in argument (pdepw) 449 398 450 399 … … 453 402 DO jj = 2, jpjm1 454 403 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)404 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 405 & * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 457 406 END DO 458 407 END DO … … 462 411 DO ji = fs_2, fs_jpim1 ! vector opt. 463 412 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)413 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 414 & * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 466 415 END DO 467 416 END DO … … 475 424 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 476 425 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') 426 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) ) & 427 & * MAX(0.,1._wp - 4.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 428 END DO 429 END DO 430 END DO 431 ENDIF 432 ! 433 IF( ln_timing ) CALL timing_stop('tke_tke') 490 434 ! 491 435 END SUBROUTINE tke_tke 492 436 493 437 494 SUBROUTINE tke_avn 438 SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt ) 495 439 !!---------------------------------------------------------------------- 496 440 !! *** ROUTINE tke_avn *** … … 524 468 !! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise. 525 469 !! 526 !! ** Action : - avt : now vertical eddy diffusivity (w-point) 527 !! - avmu, avmv : now vertical eddy viscosity at uw- and vw-points 528 !!---------------------------------------------------------------------- 470 !! ** Action : - avt, avm : now vertical eddy diffusivity and viscosity (w-point) 471 !!---------------------------------------------------------------------- 472 USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d ! ocean vertical physics 473 !! 474 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdepw ! depth (w-points) 475 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_e3t, p_e3w ! level thickness (t- & w-points) 476 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 477 ! 529 478 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, zmxld479 REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars 480 REAL(wp) :: zdku, zdkv, zsqen ! - - 481 REAL(wp) :: zemxl, zemlm, zemlp ! - - 482 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmxlm, zmxld ! 3D workspace 534 483 !!-------------------------------------------------------------------- 535 484 ! 536 IF( nn_timing == 1 ) CALL timing_start('tke_avn') 537 538 CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 539 485 IF( ln_timing ) CALL timing_start('tke_avn') 486 ! 540 487 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 541 488 ! ! Mixing length … … 549 496 ! 550 497 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 498 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 551 499 DO jj = 2, jpjm1 552 500 DO ji = fs_2, fs_jpim1 553 zraug = vkarmn * 2.e5_wp / ( rau0 * grav )554 501 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 555 502 END DO … … 576 523 ! 577 524 !!gm Not sure of that coding for ISF.... 578 ! where wmask = 0 set zmxlm == e3w_n525 ! where wmask = 0 set zmxlm == p_e3w 579 526 CASE ( 0 ) ! bounded by the distance to surface and bottom 580 527 DO jk = 2, jpkm1 581 528 DO jj = 2, jpjm1 582 529 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) )530 zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 531 & pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) ) 585 532 ! 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))533 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk)) 534 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 535 END DO 589 536 END DO … … 594 541 DO jj = 2, jpjm1 595 542 DO ji = fs_2, fs_jpim1 ! vector opt. 596 zemxl = MIN( e3w_n(ji,jj,jk), zmxlm(ji,jj,jk) )543 zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 597 544 zmxlm(ji,jj,jk) = zemxl 598 545 zmxld(ji,jj,jk) = zemxl … … 605 552 DO jj = 2, jpjm1 606 553 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) )554 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 608 555 END DO 609 556 END DO … … 612 559 DO jj = 2, jpjm1 613 560 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) )561 zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 615 562 zmxlm(ji,jj,jk) = zemxl 616 563 zmxld(ji,jj,jk) = zemxl … … 623 570 DO jj = 2, jpjm1 624 571 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) )572 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 626 573 END DO 627 574 END DO … … 630 577 DO jj = 2, jpjm1 631 578 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) )579 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 633 580 END DO 634 581 END DO … … 647 594 END SELECT 648 595 ! 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) 596 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 597 ! ! Vertical eddy viscosity and diffusivity (avm and avt) 656 598 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 657 599 DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points … … 660 602 zsqen = SQRT( en(ji,jj,jk) ) 661 603 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)604 p_avm(ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 605 p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 664 606 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 665 607 END DO 666 608 END DO 667 609 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 610 ! 679 611 ! 680 612 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt … … 682 614 DO jj = 2, jpjm1 683 615 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 616 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 617 END DO 692 618 END DO 693 619 END DO 694 620 ENDIF 695 CALL lbc_lnk( avt, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged) 696 621 ! 697 622 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') 623 CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=p_avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 624 CALL prt_ctl( tab3d_1=p_avm, clinfo1=' tke - m: ', ovlap=1, kdim=jpk ) 625 ENDIF 626 ! 627 IF( ln_timing ) CALL timing_stop('tke_avn') 706 628 ! 707 629 END SUBROUTINE tke_avn … … 722 644 !! ** Action : Increase by 1 the nstop flag is setting problem encounter 723 645 !!---------------------------------------------------------------------- 646 USE zdf_oce , ONLY : ln_zdfiwm ! Internal Wave Mixing flag 647 !! 724 648 INTEGER :: ji, jj, jk ! dummy loop indices 725 649 INTEGER :: ios 726 650 !! 727 NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , &728 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , &729 & rn_mxl0 , nn_pdl , ln_ lc , rn_lc, &651 NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , & 652 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 653 & rn_mxl0 , nn_pdl , ln_drg , ln_lc , rn_lc, & 730 654 & nn_etau , nn_htau , rn_efr 731 655 !!---------------------------------------------------------------------- … … 741 665 ! 742 666 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 667 ! 756 668 IF(lwp) THEN !* Control print … … 764 676 WRITE(numout,*) ' minimum value of tke rn_emin = ', rn_emin 765 677 WRITE(numout,*) ' surface minimum value of tke rn_emin0 = ', rn_emin0 678 WRITE(numout,*) ' prandl number flag nn_pdl = ', nn_pdl 766 679 WRITE(numout,*) ' background shear (>0) rn_bshear = ', rn_bshear 767 680 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_lc681 WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0 682 WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0 683 WRITE(numout,*) ' top/bottom friction forcing flag ln_drg = ', ln_drg 684 WRITE(numout,*) ' Langmuir cells parametrization ln_lc = ', ln_lc 685 WRITE(numout,*) ' coef to compute vertical velocity of LC rn_lc = ', rn_lc 773 686 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_efr687 WRITE(numout,*) ' type of tke penetration profile nn_htau = ', nn_htau 688 WRITE(numout,*) ' fraction of TKE that penetrates rn_efr = ', rn_efr 776 689 WRITE(numout,*) 777 WRITE(numout,*) ' critical Richardson nb with your parameters ri_cri = ', ri_cri 690 IF( ln_drg ) THEN 691 WRITE(numout,*) ' Namelist namdrg_top/_bot: used values:' 692 WRITE(numout,*) ' top ocean cavity roughness (m) rn_z0(_top)= ', r_z0_top 693 WRITE(numout,*) ' Bottom seafloor roughness (m) rn_z0(_bot)= ', r_z0_bot 694 ENDIF 695 WRITE(numout,*) 696 WRITE(numout,*) 697 WRITE(numout,*) ' ==>> critical Richardson nb with your parameters ri_cri = ', ri_cri 698 WRITE(numout,*) 699 ENDIF 700 ! 701 IF( ln_zdfiwm ) THEN ! Internal wave-driven mixing 702 rn_emin = 1.e-10_wp ! specific values of rn_emin & rmxl_min are used 703 rmxl_min = 1.e-03_wp ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s) 704 IF(lwp) WRITE(numout,*) ' Internal wave-driven mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 ' 705 ELSE ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s) 706 rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) ) ! resulting minimum length to recover molecular viscosity 707 IF(lwp) WRITE(numout,*) ' minimum mixing length with your parameters rmxl_min = ', rmxl_min 778 708 ENDIF 779 709 ! … … 786 716 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 787 717 IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 788 718 ! 789 719 IF( ln_mxl0 ) THEN 790 720 IF(lwp) WRITE(numout,*) ' use a surface mixing length = F(stress) : set rn_mxl0 = rmxl_min' … … 803 733 END SELECT 804 734 ENDIF 805 ! !* set vertical eddy coef. to the background value 806 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) 811 END DO 812 dissl(:,:,:) = 1.e-12_wp 813 ! 814 CALL tke_rst( nit000, 'READ' ) !* read or initialize all required files 735 ! !* read or initialize all required files 736 CALL tke_rst( nit000, 'READ' ) ! (en, avt_k, avm_k, dissl) 815 737 ! 816 738 END SUBROUTINE zdf_tke_init … … 818 740 819 741 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 ! 742 !!--------------------------------------------------------------------- 743 !! *** ROUTINE tke_rst *** 744 !! 745 !! ** Purpose : Read or write TKE file (en) in restart file 746 !! 747 !! ** Method : use of IOM library 748 !! if the restart does not contain TKE, en is either 749 !! set to rn_emin or recomputed 750 !!---------------------------------------------------------------------- 751 USE zdf_oce , ONLY : en, avt_k, avm_k ! ocean vertical physics 752 !! 753 INTEGER , INTENT(in) :: kt ! ocean time-step 754 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 755 ! 756 INTEGER :: jit, jk ! dummy loop indices 757 INTEGER :: id1, id2, id3, id4 ! local integers 758 !!---------------------------------------------------------------------- 759 ! 760 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 761 ! ! --------------- 762 IF( ln_rstart ) THEN !* Read the restart file 763 id1 = iom_varid( numror, 'en' , ldstop = .FALSE. ) 764 id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. ) 765 id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. ) 766 id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. ) 767 ! 768 IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! fields exist 769 CALL iom_get( numror, jpdom_autoglo, 'en', en ) 770 CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k ) 771 CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k ) 772 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl ) 773 ELSE ! start TKE from rest 774 IF(lwp) WRITE(numout,*) ' ==>> previous run without TKE scheme, set en to background values' 775 en (:,:,:) = rn_emin * wmask(:,:,:) 776 dissl(:,:,:) = 1.e-12_wp 777 ! avt_k, avm_k already set to the background value in zdf_phy_init 778 ENDIF 779 ELSE !* Start from rest 780 IF(lwp) WRITE(numout,*) ' ==>> start from rest: set en to the background value' 781 en (:,:,:) = rn_emin * wmask(:,:,:) 782 dissl(:,:,:) = 1.e-12_wp 783 ! avt_k, avm_k already set to the background value in zdf_phy_init 784 ENDIF 785 ! 786 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 787 ! ! ------------------- 788 IF(lwp) WRITE(numout,*) '---- tke-rst ----' 789 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 790 CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k ) 791 CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k ) 792 CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl ) 793 ! 794 ENDIF 795 ! 891 796 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 797 910 798 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.