- Timestamp:
- 2019-08-06T15:36:27+02:00 (5 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex
r11388 r11407 1329 1329 \label{eq:Eqn_zad_Aimp_partition} 1330 1330 Cu_{min} &= 0.15 \nonumber \\ 1331 Cu_{max} &= 0. 27\nonumber \\1331 Cu_{max} &= 0.3 \nonumber \\ 1332 1332 Cu_{cut} &= 2Cu_{max} - Cu_{min} \nonumber \\ 1333 1333 Fcu &= 4Cu_{max}*(Cu_{max}-Cu_{min}) \nonumber \\ 1334 C f &=1334 C\kern-0.14em f &= 1335 1335 \begin{cases} 1336 1336 0.0 &\text{if $Cu \leq Cu_{min}$} \\ … … 1340 1340 \end{align} 1341 1341 1342 \noindent With these settings the coefficient ($Cf$) is shown in \autoref{fig:zad_Aimp_coeff}1343 1344 1342 \begin{figure}[!t] 1345 1343 \begin{center} … … 1347 1345 \caption{ 1348 1346 \protect\label{fig:zad_Aimp_coeff} 1349 The value of the partitioning coefficient ($C f$) used to partition vertical velocities into parts to1347 The value of the partitioning coefficient ($C\kern-0.14em f$) used to partition vertical velocities into parts to 1350 1348 be treated implicitly and explicitly for a range of typical Courant numbers (\forcode{ln_zad_Aimp=.true.}) 1351 1349 } … … 1359 1357 \begin{align} 1360 1358 \label{eq:Eqn_zad_Aimp_partition2} 1361 w_{i_{ijk}} &= C f_{ijk} w_{n_{ijk}} \nonumber \\1362 w_{n_{ijk}} &= (1-C f_{ijk}) w_{n_{ijk}}1359 w_{i_{ijk}} &= C\kern-0.14em f_{ijk} w_{n_{ijk}} \nonumber \\ 1360 w_{n_{ijk}} &= (1-C\kern-0.14em f_{ijk}) w_{n_{ijk}} 1363 1361 \end{align} 1364 1362 1365 1363 \noindent Note that the coefficient is such that the treatment is never fully implicit; 1366 1364 the three cases from \autoref{eq:Eqn_zad_Aimp_partition} can be considered as: 1367 fully-explicit; mixed explicit/implicit and mostly-implicit. The $w_i$ component is added 1368 to the implicit solvers for the vertical mixing in \mdl{dynzdf} and \mdl{trazdf} in a 1369 similar way to \citep{shchepetkin_OM15}. This is sufficient for the flux-limited 1370 advection scheme (\forcode{ln_traadv_mus}) but further intervention is required when using 1371 the flux-corrected scheme (\forcode{ln_traadv_fct}). For these schemes the implicit 1372 upstream fluxes must be added to both the monotonic guess and to the higher order solution 1373 when calculating the antidiffusive fluxes. The implicit vertical fluxes are then removed 1374 since they are added by the implicit solver later on. 1365 fully-explicit; mixed explicit/implicit and mostly-implicit. With the settings shown the 1366 coefficient ($C\kern-0.14em f$) varies as shown in \autoref{fig:zad_Aimp_coeff}. Note with these values 1367 the $Cu_{cut}$ boundary between the mixed implicit-explicit treatment and 'mostly 1368 implicit' is 0.45 which is just below the stability limited given in 1369 \autoref{tab:zad_Aimp_CFLcrit} for a 3rd order scheme. 1370 1371 The $w_i$ component is added to the implicit solvers for the vertical mixing in 1372 \mdl{dynzdf} and \mdl{trazdf} in a similar way to \citep{shchepetkin_OM15}. This is 1373 sufficient for the flux-limited advection scheme (\forcode{ln_traadv_mus}) but further 1374 intervention is required when using the flux-corrected scheme (\forcode{ln_traadv_fct}). 1375 For these schemes the implicit upstream fluxes must be added to both the monotonic guess 1376 and to the higher order solution when calculating the antidiffusive fluxes. The implicit 1377 vertical fluxes are then removed since they are added by the implicit solver later on. 1375 1378 1376 1379 The adaptive-implicit vertical advection option is new to NEMO at v4.0 and has yet to be 1377 1380 used in a wide range of simulations. The following test simulation, however, does illustrate 1378 1381 the potential benefits and will hopefully encourage further testing and feedback from users: 1382 1383 \begin{figure}[!t] 1384 \begin{center} 1385 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_overflow_frames} 1386 \caption{ 1387 \protect\label{fig:zad_Aimp_overflow_frames} 1388 A time-series of temperature vertical cross-sections for the OVERFLOW test case. These results are for the default 1389 settings with \forcode{nn_rdt=10.0} and without adaptive implicit vertical advection (\forcode{ln_zad_Aimp=.false.}). 1390 } 1391 \end{center} 1392 \end{figure} 1379 1393 1380 1394 \subsection{Adaptive-implicit vertical advection in the OVERFLOW test-case} … … 1389 1403 ln_zdfnpc = .true. 1390 1404 ln_traadv_fct = .true. 1391 nn_fct_h = 41392 nn_fct_v = 41405 nn_fct_h = 2 1406 nn_fct_v = 2 1393 1407 \end{verbatim} 1394 1408 … … 1404 1418 candidate, therefore, for use of the adaptive-implicit vertical advection scheme. 1405 1419 1406 The results with \forcode{ln_zad_Aimp=.true.} and a variety of model timesteps are shown 1407 in \autoref{fig:zad_Aimp_overflow_all_rdt} (together with the equivalent frames from the 1408 base run). In this simple example the use of the adaptive-implicit vertcal advection 1409 scheme has enabled a 12x increase in the model timestep without significantly altering the 1410 solution (although at this extreme the plume is more diffuse and has not travelled so far). 1411 Notably, the solution with and without the scheme is slightly different even with 1412 \forcode{nn_rdt = 10.}; suggesting that the base run was close enough to instability to 1413 trigger the scheme despite completing successfully. To assist in diagnosing how active 1414 the scheme is, in both location and time, the 3D implicit and explicit components of the 1415 vertical velocity are available via XIOS as \texttt{wimp} and \texttt{wexp} respectively. 1416 Likewise, the partitioning coefficient ($Cf$) is also available as \texttt{wi\_cff}. 1417 1418 \begin{figure}[!t] 1419 \begin{center} 1420 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_overflow_frames} 1421 \caption{ 1422 \protect\label{fig:zad_Aimp_overflow_frames} 1423 A time-series of temperature vertical cross-sections for the OVERFLOW test case. These results are for the default 1424 settings with \forcode{nn_rdt=10.0} and without adaptive implicit vertical advection (\forcode{ln_zad_Aimp=.false.}). 1425 } 1426 \end{center} 1427 \end{figure} 1420 The results with \forcode{ln_zad_Aimp=.true.} and a variety of model timesteps 1421 are shown in \autoref{fig:zad_Aimp_overflow_all_rdt} (together with the equivalent 1422 frames from the base run). In this simple example the use of the adaptive-implicit 1423 vertcal advection scheme has enabled a 12x increase in the model timestep without 1424 significantly altering the solution (although at this extreme the plume is more diffuse 1425 and has not travelled so far). Notably, the solution with and without the scheme is 1426 slightly different even with \forcode{nn_rdt = 10.}; suggesting that the base run was 1427 close enough to instability to trigger the scheme despite completing successfully. 1428 To assist in diagnosing how active the scheme is, in both location and time, the 3D 1429 implicit and explicit components of the vertical velocity are available via XIOS as 1430 \texttt{wimp} and \texttt{wexp} respectively. Likewise, the partitioning coefficient 1431 ($C\kern-0.14em f$) is also available as \texttt{wi\_cff}. For a quick oversight of 1432 the schemes activity the global maximum values of the absolute implicit component 1433 of the vertical velocity and the partitioning coefficient are written to the netCDF 1434 version of the run statistics file (\texttt{run.stat.nc}) if this is active (see 1435 \autoref{sec:MISC_opt} for activation details). 1436 1437 \autoref{fig:zad_Aimp_maxCf} shows examples of the maximum partitioning coefficient for 1438 the various overflow tests. Note that the adaptive-implicit vertical advection scheme is 1439 active even in the base run with \forcode{nn_rdt=10.0s} adding to the evidence that the 1440 test case is close to stability limits even with this value. At the larger timesteps, the 1441 vertical velocity is treated mostly implicitly at some location throughout the run. The 1442 oscillatory nature of this measure appears to be linked to the progress of the plume front 1443 as each cusp is associated with the location of the maximum shifting to the adjacent cell. 1444 This is illustrated in \autoref{fig:zad_Aimp_maxCf_loc} where the i- and k- locations of the 1445 maximum have been overlaid for the base run case. 1446 1447 \medskip 1448 \noindent Only limited tests have been performed in more realistic configurations. In the 1449 ORCA2\_ICE\_PISCES reference configuration the scheme does activate and passes 1450 restartability and reproducibility tests but it is unable to improve the model's stability 1451 enough to allow an increase in the model time-step. A view of the time-series of maximum 1452 partitioning coefficient (not shown here) suggests that the default time-step of 5400s is 1453 already pushing at stability limits, especially in the initial start-up phase. The 1454 time-series does not, however, exhibit any of the 'cuspiness' found with the overflow 1455 tests. 1456 1457 \medskip 1458 \noindent A short test with an eORCA1 configuration promises more since a test using a 1459 time-step of 3600s remains stable with \forcode{ln_zad_Aimp=.true.} whereas the 1460 time-step is limited to 2700s without. 1428 1461 1429 1462 \begin{figure}[!t] … … 1440 1473 \end{figure} 1441 1474 1475 \begin{figure}[!t] 1476 \begin{center} 1477 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_maxCf} 1478 \caption{ 1479 \protect\label{fig:zad_Aimp_maxCf} 1480 The maximum partitioning coefficient during a series of test runs with increasing model timestep length. 1481 At the larger timesteps, the vertical velocity is treated mostly implicitly at some location throughout 1482 the run. 1483 } 1484 \end{center} 1485 \end{figure} 1486 1487 \begin{figure}[!t] 1488 \begin{center} 1489 \includegraphics[width=\textwidth]{Fig_ZDF_zad_Aimp_maxCf_loc} 1490 \caption{ 1491 \protect\label{fig:zad_Aimp_maxCf_loc} 1492 The maximum partitioning coefficient for the \forcode{nn_rdt=10.0s} case overlaid with information on the gridcell i- and k- 1493 locations of the maximum value. 1494 } 1495 \end{center} 1496 \end{figure} 1442 1497 1443 1498 % ================================================================ -
NEMO/trunk/src/OCE/DYN/sshwzv.F90
r11293 r11407 286 286 REAL(wp) :: zCu, zcff, z1_e3t ! local scalars 287 287 REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters 288 REAL(wp) , PARAMETER :: Cu_max = 0. 27! local parameters288 REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters 289 289 REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters 290 290 REAL(wp) , PARAMETER :: Fcu = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters … … 346 346 wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 347 347 ! 348 Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient 348 Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl 349 349 END DO 350 350 END DO … … 353 353 ELSE 354 354 ! Fully explicit everywhere 355 Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient 355 Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl 356 356 wi (:,:,:) = 0._wp 357 357 ENDIF -
NEMO/trunk/src/OCE/TRA/traadv_fct.F90
r10425 r11407 21 21 USE diaar5 ! AR5 diagnostics 22 22 USE phycst , ONLY : rau0_rcp 23 USE zdf_oce , ONLY : ln_zad_Aimp 23 24 ! 24 25 USE in_out_manager ! I/O manager … … 86 87 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 87 88 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry 89 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup 90 LOGICAL :: ll_zAimp ! flag to apply adaptive implicit vertical advection 88 91 !!---------------------------------------------------------------------- 89 92 ! … … 97 100 l_hst = .FALSE. 98 101 l_ptr = .FALSE. 102 ll_zAimp = .FALSE. 99 103 IF( ( cdtype =='TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 100 104 IF( cdtype =='TRA' .AND. ln_diaptr ) l_ptr = .TRUE. … … 116 120 ! 117 121 zwi(:,:,:) = 0._wp 122 ! 123 ! If adaptive vertical advection, check if it is needed on this PE at this time 124 IF( ln_zad_Aimp ) THEN 125 IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE. 126 END IF 127 ! If active adaptive vertical advection, build tridiagonal matrix 128 IF( ll_zAimp ) THEN 129 ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 130 DO jk = 1, jpkm1 131 DO jj = 2, jpjm1 132 DO ji = fs_2, fs_jpim1 ! vector opt. (ensure same order of calculation as below if wi=0.) 133 zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) / e3t_a(ji,jj,jk) 134 zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t_a(ji,jj,jk) 135 zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t_a(ji,jj,jk) 136 END DO 137 END DO 138 END DO 139 END IF 118 140 ! 119 141 DO jn = 1, kjpt !== loop over the tracers ==! … … 169 191 END DO 170 192 END DO 193 194 IF ( ll_zAimp ) THEN 195 CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 ) 196 ! 197 ztu(:,:,1) = 0._wp; ztu(:,:,jpk) = 0._wp 198 DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask) 199 DO jj = 1, jpj 200 DO ji = 1, jpi 201 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 202 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 203 ztu(ji,jj,jk) = 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 204 END DO 205 END DO 206 END DO 207 DO jk = 1, jpkm1 208 DO jj = 2, jpjm1 209 DO ji = fs_2, fs_jpim1 ! vector opt. 210 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztu(ji,jj,jk) - ztu(ji ,jj ,jk+1) ) & 211 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 212 END DO 213 END DO 214 END DO 215 zwz(:,:,:) = zwz(:,:,:) + ztu(:,:,:) 216 ! 217 END IF 171 218 ! 172 219 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) … … 280 327 CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1., zwz, 'W', 1. ) 281 328 ! 329 ! 330 IF ( ll_zAimp ) THEN 331 DO jk = 1, jpkm1 !* trend and after field with monotonic scheme 332 DO jj = 2, jpjm1 333 DO ji = fs_2, fs_jpim1 ! vector opt. 334 ! ! total intermediate advective trends 335 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 336 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 337 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 338 ztu(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 339 END DO 340 END DO 341 END DO 342 ! 343 CALL tridia_solver( zwdia, zwsup, zwinf, ztu, ztu , 0 ) 344 ! 345 ztu(:,:,1) = 0._wp 346 DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask) 347 DO jj = 2, jpjm1 348 DO ji = fs_2, fs_jpim1 ! vector opt. 349 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 350 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 351 zwz(ji,jj,jk) = zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztu(ji,jj,jk) + zfm_wk * ztu(ji,jj,jk-1) ) * wmask(ji,jj,jk) 352 END DO 353 END DO 354 END DO 355 END IF 282 356 ! !== monotonicity algorithm ==! 283 357 ! … … 289 363 DO jj = 2, jpjm1 290 364 DO ji = fs_2, fs_jpim1 ! vector opt. 291 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 292 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 293 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) & 294 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 295 END DO 296 END DO 297 END DO 365 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 366 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 367 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 368 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk) 369 zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 370 END DO 371 END DO 372 END DO 373 ! 374 IF ( ll_zAimp ) THEN 375 ! 376 DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask) 377 DO jj = 2, jpjm1 378 DO ji = fs_2, fs_jpim1 ! vector opt. 379 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 380 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 381 zwz(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 382 END DO 383 END DO 384 END DO 385 DO jk = 1, jpkm1 386 DO jj = 2, jpjm1 387 DO ji = fs_2, fs_jpim1 ! vector opt. 388 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) & 389 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 390 END DO 391 END DO 392 END DO 393 END IF 298 394 ! 299 395 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport … … 318 414 END DO ! end of tracer loop 319 415 ! 416 IF ( ll_zAimp ) THEN 417 DEALLOCATE( zwdia, zwinf, zwsup ) 418 ENDIF 320 419 IF( l_trd .OR. l_hst ) THEN 321 420 DEALLOCATE( ztrdx, ztrdy, ztrdz ) -
NEMO/trunk/src/OCE/stpctl.F90
r10570 r11407 96 96 IF( ln_zad_Aimp ) THEN 97 97 istatus = NF90_DEF_VAR( idrun, 'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1 ) 98 istatus = NF90_DEF_VAR( idrun, 'C u_max', NF90_DOUBLE, (/ idtime /), idc1 )98 istatus = NF90_DEF_VAR( idrun, 'Cf_max', NF90_DOUBLE, (/ idtime /), idc1 ) 99 99 ENDIF 100 100 istatus = NF90_ENDDEF(idrun) … … 123 123 IF( ln_zad_Aimp ) THEN 124 124 zmax(8) = MAXVAL( ABS( wi(:,:,:) ) , mask = wmask(:,:,:) == 1._wp ) ! implicit vertical vel. max 125 zmax(9) = MAXVAL( Cu_adv(:,:,:) , mask = tmask(:,:,:) == 1._wp ) ! cell Courant no. max125 zmax(9) = MAXVAL( Cu_adv(:,:,:) , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max 126 126 ENDIF 127 127 !
Note: See TracChangeset
for help on using the changeset viewer.