Changeset 7412 for branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/LDF
- Timestamp:
- 2016-12-01T11:30:29+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
r6140 r7412 42 42 REAL(wp), PUBLIC :: rn_ahm_b !: lateral laplacian background eddy viscosity [m2/s] 43 43 REAL(wp), PUBLIC :: rn_bhm_0 !: lateral bilaplacian eddy viscosity [m4/s] 44 !! If nn_ahm_ijk_t = 32 a time and space varying Smagorinsky viscosity 45 !! will be computed. 46 REAL(wp), PUBLIC :: rn_csmc !: Smagorinsky constant of proportionality 47 REAL(wp), PUBLIC :: rn_minfac !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 48 REAL(wp), PUBLIC :: rn_maxfac !: Multiplicative factor of theorectical maximum Smagorinsky viscosity 44 49 45 50 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 46 51 47 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy diffusivity coef. at U- and V-points [m2/s or m4/s] 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dtensq !: horizontal tension squared (Smagorinsky only) 54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dshesq !: horizontal shearing strain squared (Smagorinsky only) 55 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: esqt, esqf !: Square of the local gridscale (e1e2/(e1+e2))**2 48 56 49 57 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 50 58 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 59 REAL(wp) :: r1_8 = 0.125_wp ! =1/8 51 60 REAL(wp) :: r1_288 = 1._wp / 288._wp ! =1/( 12^2 * 2 ) 52 61 … … 79 88 !! = 31 = F(i,j,k,t) = F(local velocity) ( |u|e /12 laplacian operator 80 89 !! or |u|e^3/12 bilaplacian operator ) 81 !!---------------------------------------------------------------------- 82 INTEGER :: jk ! dummy loop indices 90 !! = 32 = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 91 !! ( L^2|D| laplacian operator 92 !! or L^4|D|/8 bilaplacian operator ) 93 !!---------------------------------------------------------------------- 94 INTEGER :: ji, jj, jk ! dummy loop indices 83 95 INTEGER :: ierr, inum, ios ! local integer 84 96 REAL(wp) :: zah0 ! local scalar … … 86 98 NAMELIST/namdyn_ldf/ ln_dynldf_lap, ln_dynldf_blp, & 87 99 & ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso, & 88 & nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0 100 & nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0, & 101 & rn_csmc , rn_minfac, rn_maxfac 89 102 !!---------------------------------------------------------------------- 90 103 ! … … 115 128 WRITE(numout,*) ' coefficients :' 116 129 WRITE(numout,*) ' type of time-space variation nn_ahm_ijk_t = ', nn_ahm_ijk_t 117 WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0 _lap= ', rn_ahm_0, ' m2/s'130 WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0 = ', rn_ahm_0, ' m2/s' 118 131 WRITE(numout,*) ' background viscosity (iso case) rn_ahm_b = ', rn_ahm_b, ' m2/s' 119 WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_ahm_0_blp = ', rn_bhm_0, ' m4/s' 132 WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_bhm_0 = ', rn_bhm_0, ' m4/s' 133 WRITE(numout,*) ' smagorinsky settings (nn_ahm_ijk_t = 32) :' 134 WRITE(numout,*) ' Smagorinsky coefficient rn_csmc = ', rn_csmc 135 WRITE(numout,*) ' factor multiplier for theorectical lower limit for ' 136 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_minfac = ', rn_minfac 137 WRITE(numout,*) ' factor multiplier for theorectical lower upper for ' 138 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_maxfac = ', rn_maxfac 120 139 ENDIF 121 140 … … 208 227 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 209 228 ! 229 CASE( 32 ) !== time varying 3D field ==! 230 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth , time )' 231 IF(lwp) WRITE(numout,*) ' proportional to the local deformation rate and gridscale (Smagorinsky)' 232 IF(lwp) WRITE(numout,*) ' : L^2|D| or L^4|D|/8' 233 ! 234 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 235 ! 236 ! allocate arrays used in ldf_dyn. 237 ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr ) 238 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays') 239 ! 240 ! Set local gridscale values 241 DO jj = 2, jpjm1 242 DO ji = fs_2, fs_jpim1 243 esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2 244 esqf(ji,jj) = ( e1e2f(ji,jj) /( e1f(ji,jj) + e2f(ji,jj) ) )**2 245 END DO 246 END DO 247 ! 210 248 CASE DEFAULT 211 249 CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm') … … 232 270 !! nn_ahm_ijk_t = 31 ahmt, ahmf = F(i,j,k,t) = F(local velocity) 233 271 !! ( |u|e /12 or |u|e^3/12 for laplacian or bilaplacian operator ) 234 !! BLP case : sqrt of the eddy coef, since bilaplacian is en re-entrant laplacian235 272 !! 236 !! ** action : ahmt, ahmf update at each time step 273 !! nn_ahm_ijk_t = 32 ahmt, ahmf = F(i,j,k,t) = F(local deformation rate and gridscale) (D and L) (Smagorinsky) 274 !! ( L^2|D| or L^4|D|/8 for laplacian or bilaplacian operator ) 275 !! 276 !! ** note : in BLP cases the sqrt of the eddy coef is returned, since bilaplacian is en re-entrant laplacian 277 !! ** action : ahmt, ahmf updated at each time step 237 278 !!---------------------------------------------------------------------- 238 279 INTEGER, INTENT(in) :: kt ! time step index … … 240 281 INTEGER :: ji, jj, jk ! dummy loop indices 241 282 REAL(wp) :: zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax ! local scalar 283 REAL(wp) :: zcmsmag, zstabf_lo, zstabf_up, zdelta, zdb ! local scalar 242 284 !!---------------------------------------------------------------------- 243 285 ! … … 248 290 CASE( 31 ) !== time varying 3D field ==! = F( local velocity ) 249 291 ! 250 IF( ln_dynldf_lap ) THEN !laplacian operator : |u| e /12 = |u/144| e292 IF( ln_dynldf_lap ) THEN ! laplacian operator : |u| e /12 = |u/144| e 251 293 DO jk = 1, jpkm1 252 294 DO jj = 2, jpjm1 … … 280 322 CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. ) 281 323 ! 324 ! 325 CASE( 32 ) !== time varying 3D field ==! = F( local deformation rate and gridscale ) (Smagorinsky) 326 ! 327 IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! laplacian operator : (C_smag/pi)^2 L^2 |D| 328 ! 329 zcmsmag = (rn_csmc/rpi)**2 ! (C_smag/pi)^2 330 zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag ) ! lower limit stability factor scaling 331 zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rdt ) ! upper limit stability factor scaling 332 IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo ! provide |U|L^3/12 lower limit instead 333 ! ! of |U|L^3/16 in blp case 334 DO jk = 1, jpkm1 335 ! 336 DO jj = 2, jpj 337 DO ji = 2, jpi 338 zdb = ( ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) & 339 & * r1_e1t(ji,jj) * e2t(ji,jj) & 340 & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) & 341 & * r1_e2t(ji,jj) * e1t(ji,jj) ) * tmask(ji,jj,jk) 342 dtensq(ji,jj) = zdb*zdb 343 END DO 344 END DO 345 ! 346 DO jj = 1, jpjm1 347 DO ji = 1, jpim1 348 zdb = ( ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) & 349 & * r1_e2f(ji,jj) * e1f(ji,jj) & 350 & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) & 351 & * r1_e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,jk) 352 dshesq(ji,jj) = zdb*zdb 353 END DO 354 END DO 355 ! 356 DO jj = 2, jpjm1 357 DO ji = fs_2, fs_jpim1 358 ! 359 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 360 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 361 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 362 ! T-point value 363 zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2 364 ahmt(ji,jj,jk) = zdelta * sqrt( dtensq(ji,jj) + & 365 & r1_4 * ( dshesq(ji,jj) + dshesq(ji,jj-1) + & 366 & dshesq(ji-1,jj) + dshesq(ji-1,jj-1) ) ) 367 ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), & 368 & SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 369 ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 370 ! F-point value 371 zdelta = zcmsmag * esqf(ji,jj) ! L^2 * (C_smag/pi)^2 372 ahmf(ji,jj,jk) = zdelta * sqrt( dshesq(ji,jj) + & 373 & r1_4 * ( dtensq(ji,jj) + dtensq(ji,jj+1) + & 374 & dtensq(ji+1,jj) + dtensq(ji+1,jj+1) ) ) 375 ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), & 376 & SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 377 ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 378 ! 379 END DO 380 END DO 381 END DO 382 ! 383 ENDIF 384 ! 385 IF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8) 386 ! = sqrt( A_lap_smag L^2/8 ) 387 ! stability limits already applied to laplacian values 388 ! effective default limits are |U|L^3/12 < B_hm < L^4/(32*2dt) 389 ! 390 DO jk = 1, jpkm1 391 ! 392 DO jj = 2, jpjm1 393 DO ji = fs_2, fs_jpim1 394 ! 395 ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 396 ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 397 ! 398 END DO 399 END DO 400 END DO 401 ! 402 ENDIF 403 ! 404 CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. ) 405 ! 282 406 END SELECT 283 407 !
Note: See TracChangeset
for help on using the changeset viewer.