Changeset 6383
- Timestamp:
- 2016-03-10T17:06:52+01:00 (8 years ago)
- Location:
- branches/2016/dev_r6381_SIMPLIF_5/NEMOGCM
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_r6381_SIMPLIF_5/NEMOGCM/CONFIG/SHARED/namelist_ref
r6152 r6383 815 815 rn_aht_0 = 2000. ! lateral eddy diffusivity (lap. operator) [m2/s] 816 816 rn_bht_0 = 1.e+12 ! lateral eddy diffusivity (bilap. operator) [m4/s] 817 ! 818 ! Smagorinsky settings: 819 rn_csmc = 3.5 ! Smagorinsky constant of proportionality 820 rn_cfacmin = 1.0 ! multiplier of theorectical lower limit 821 rn_cfacmax = 1.0 ! multiplier of theorectical upper limit 817 822 / 818 823 !----------------------------------------------------------------------- -
branches/2016/dev_r6381_SIMPLIF_5/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
r6140 r6383 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_cfacmin !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 48 REAL(wp), PUBLIC :: rn_cfacmax !: 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_cfacmin, rn_cfacmax 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_cfacmin = ', rn_cfacmin 137 WRITE(numout,*) ' factor multiplier for theorectical lower upper for ' 138 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_cfacmax = ', rn_cfacmax 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 locate gradient arrays used in ldf_dyn. Allocated once here for efficiency. 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_cfacmin * rn_cfacmin / ( 2_.wp * 4._wp * zcmsmag ) ! lower limit stability factor scaling 331 zstabf_up = rn_cfacmax / ( 4._wp * zcmsmag * rdt ) ! upper limit stability factor scaling 332 ! 333 DO jk = 1, jpkm1 334 ! 335 DO jj = 2, jpj 336 DO ji = 2, jpi 337 zdb = ( ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) & 338 & * r1_e1t(ji,jj) * e2t(ji,jj) & 339 & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) & 340 & * r1_e2t(ji,jj) * e1t(ji,jj) ) * tmask(ji,jj,jk) 341 dtensq(ji,jj) = zdb*zdb 342 END DO 343 END DO 344 ! 345 DO jj = 1, jpjm1 346 DO ji = 1, jpim1 347 zdb = ( ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) & 348 & * r1_e2f(ji,jj) * e1f(ji,jj) & 349 & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) & 350 & * r1_e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,jk) 351 dshesq(ji,jj) = zdb*zdb 352 END DO 353 END DO 354 ! 355 DO jj = 2, jpjm1 356 DO ji = fs_2, fs_jpim1 357 ! 358 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 359 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 360 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 361 ! T-point value 362 zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2 363 ahmt(ji,jj,jk) = zdelta * sqrt( dtensq(ji,jj) + & 364 & r1_4 * ( dshesq(ji,jj) + dshesq(ji,jj-1) + & 365 & dshesq(ji-1,jj) + dshesq(ji-1,jj-1) ) ) 366 ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), & 367 & SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == cfacmin * |U|L/2 368 ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == cfacmax * L^2/(4dt) 369 ! F-point value 370 zdelta = zcmsmag * esqf(ji,jj) ! L^2 * (C_smag/pi)^2 371 ahmf(ji,jj,jk) = zdelta * sqrt( dshesq(ji,jj) + & 372 & r1_4 * ( dtensq(ji,jj) + dtensq(ji,jj+1) + & 373 & dtensq(ji+1,jj) + dtensq(ji+1,jj+1) ) ) 374 ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), & 375 & SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == cfacmin * |U|L/2 376 ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == cfacmax * L^2/(4dt) 377 ! 378 END DO 379 END DO 380 END DO 381 ! 382 ENDIF 383 ! 384 IF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8) 385 ! = sqrt( A_lap_smag L^2/8 ) 386 ! stability limits already applied to laplacian values 387 ! 388 DO jk = 1, jpkm1 389 ! 390 DO jj = 2, jpjm1 391 DO ji = fs_2, fs_jpim1 392 ! 393 ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 394 ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 395 ! 396 END DO 397 END DO 398 END DO 399 ! 400 ENDIF 401 ! 402 CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. ) 403 ! 282 404 END SELECT 283 405 !
Note: See TracChangeset
for help on using the changeset viewer.