Changeset 11407


Ignore:
Timestamp:
2019-08-06T15:36:27+02:00 (15 months ago)
Author:
acc
Message:

Final documentation tweaks to the adaptive-implicit vertical advection section of chap_ZDF.tex and associated code changes. The code changes are mostly cosmetic or to align with the documentation but this does include a modified version of traadv_fct.F90 (thanks to Jerome) which correctly accounts for any implicit component of the vertical velocity when computing anti-diffusive fluxes for flux correction. The changes have null effect when ln_zad_Aimp is .false. and have been SETTE tested with ORCA2_ICE_PISCES with ln_zad_Aimp = .true.

Location:
NEMO/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex

    r11388 r11407  
    13291329  \label{eq:Eqn_zad_Aimp_partition} 
    13301330Cu_{min} &= 0.15 \nonumber \\ 
    1331 Cu_{max} &= 0.27 \nonumber \\ 
     1331Cu_{max} &= 0.3 \nonumber \\ 
    13321332Cu_{cut} &= 2Cu_{max} - Cu_{min} \nonumber \\ 
    13331333Fcu    &= 4Cu_{max}*(Cu_{max}-Cu_{min}) \nonumber \\ 
    1334 Cf &= 
     1334C\kern-0.14em f &= 
    13351335     \begin{cases} 
    13361336        0.0                                                        &\text{if $Cu \leq Cu_{min}$} \\ 
     
    13401340\end{align} 
    13411341 
    1342 \noindent With these settings the coefficient ($Cf$) is shown in \autoref{fig:zad_Aimp_coeff} 
    1343  
    13441342\begin{figure}[!t] 
    13451343  \begin{center} 
     
    13471345    \caption{ 
    13481346      \protect\label{fig:zad_Aimp_coeff} 
    1349       The value of the partitioning coefficient ($Cf$) used to partition vertical velocities into parts to 
     1347      The value of the partitioning coefficient ($C\kern-0.14em f$) used to partition vertical velocities into parts to 
    13501348      be treated implicitly and explicitly for a range of typical Courant numbers (\forcode{ln_zad_Aimp=.true.}) 
    13511349    } 
     
    13591357\begin{align} 
    13601358  \label{eq:Eqn_zad_Aimp_partition2} 
    1361     w_{i_{ijk}} &= Cf_{ijk} w_{n_{ijk}}     \nonumber \\ 
    1362     w_{n_{ijk}} &= (1-Cf_{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}}            
    13631361\end{align} 
    13641362 
    13651363\noindent Note that the coefficient is such that the treatment is never fully implicit; 
    13661364the 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. 
     1365fully-explicit; mixed explicit/implicit and mostly-implicit.  With the settings shown the 
     1366coefficient ($C\kern-0.14em f$) varies as shown in \autoref{fig:zad_Aimp_coeff}. Note with these values 
     1367the $Cu_{cut}$ boundary between the mixed implicit-explicit treatment and 'mostly 
     1368implicit' is 0.45 which is just below the stability limited given in 
     1369\autoref{tab:zad_Aimp_CFLcrit}  for a 3rd order scheme. 
     1370 
     1371The $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 
     1373sufficient for the flux-limited advection scheme (\forcode{ln_traadv_mus}) but further 
     1374intervention is required when using the flux-corrected scheme (\forcode{ln_traadv_fct}). 
     1375For these schemes the implicit upstream fluxes must be added to both the monotonic guess 
     1376and to the higher order solution when calculating the antidiffusive fluxes. The implicit 
     1377vertical fluxes are then removed since they are added by the implicit solver later on. 
    13751378 
    13761379The adaptive-implicit vertical advection option is new to NEMO at v4.0 and has yet to be  
    13771380used in a wide range of simulations. The following test simulation, however, does illustrate 
    13781381the 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} 
    13791393 
    13801394\subsection{Adaptive-implicit vertical advection in the OVERFLOW test-case} 
     
    13891403     ln_zdfnpc     = .true. 
    13901404     ln_traadv_fct = .true. 
    1391         nn_fct_h   =  4 
    1392         nn_fct_v   =  4 
     1405        nn_fct_h   =  2 
     1406        nn_fct_v   =  2 
    13931407\end{verbatim} 
    13941408 
     
    14041418candidate, therefore, for use of the adaptive-implicit vertical advection scheme. 
    14051419 
    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} 
     1420The results with \forcode{ln_zad_Aimp=.true.} and a variety of model timesteps 
     1421are shown in \autoref{fig:zad_Aimp_overflow_all_rdt} (together with the equivalent 
     1422frames from the base run).  In this simple example the use of the adaptive-implicit 
     1423vertcal advection scheme has enabled a 12x increase in the model timestep without 
     1424significantly altering the solution (although at this extreme the plume is more diffuse 
     1425and has not travelled so far).  Notably, the solution with and without the scheme is 
     1426slightly different even with \forcode{nn_rdt = 10.}; suggesting that the base run was 
     1427close enough to instability to trigger the scheme despite completing successfully. 
     1428To assist in diagnosing how active the scheme is, in both location and time, the 3D 
     1429implicit 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 
     1432the schemes activity the global maximum values of the absolute implicit component 
     1433of the vertical velocity and the partitioning coefficient are written to the netCDF 
     1434version 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 
     1438the various overflow tests.  Note that the adaptive-implicit vertical advection scheme is 
     1439active even in the base run with \forcode{nn_rdt=10.0s} adding to the evidence that the 
     1440test case is close to stability limits even with this value. At the larger timesteps, the 
     1441vertical velocity is treated mostly implicitly at some location throughout the run. The 
     1442oscillatory nature of this measure appears to be linked to the progress of the plume front 
     1443as each cusp is associated with the location of the maximum shifting to the adjacent cell. 
     1444This is illustrated in \autoref{fig:zad_Aimp_maxCf_loc} where the i- and k- locations of the 
     1445maximum 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 
     1449ORCA2\_ICE\_PISCES reference configuration the scheme does activate and passes 
     1450restartability and reproducibility tests but it is unable to improve the model's stability 
     1451enough to allow an increase in the model time-step. A view of the time-series of maximum 
     1452partitioning coefficient (not shown here)  suggests that the default time-step of 5400s is 
     1453already pushing at stability limits, especially in the initial start-up phase. The 
     1454time-series does not, however, exhibit any of the 'cuspiness' found with the overflow 
     1455tests. 
     1456 
     1457\medskip 
     1458\noindent A short test with an eORCA1 configuration promises more since a test using a 
     1459time-step of 3600s remains stable with \forcode{ln_zad_Aimp=.true.} whereas the 
     1460time-step is limited to 2700s without. 
    14281461 
    14291462\begin{figure}[!t] 
     
    14401473\end{figure} 
    14411474 
     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} 
    14421497 
    14431498% ================================================================ 
  • NEMO/trunk/src/OCE/DYN/sshwzv.F90

    r11293 r11407  
    286286      REAL(wp)             ::   zCu, zcff, z1_e3t                     ! local scalars 
    287287      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters 
    288       REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters 
     288      REAL(wp) , PARAMETER ::   Cu_max = 0.30_wp                      ! local parameters 
    289289      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters 
    290290      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters 
     
    346346                  wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 
    347347                  ! 
    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 
    349349               END DO 
    350350            END DO 
     
    353353      ELSE 
    354354         ! 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 
    356356         wi    (:,:,:) = 0._wp 
    357357      ENDIF 
  • NEMO/trunk/src/OCE/TRA/traadv_fct.F90

    r10425 r11407  
    2121   USE diaar5         ! AR5 diagnostics 
    2222   USE phycst  , ONLY : rau0_rcp 
     23   USE zdf_oce , ONLY : ln_zad_Aimp 
    2324   ! 
    2425   USE in_out_manager ! I/O manager 
     
    8687      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 
    8788      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 
    8891      !!---------------------------------------------------------------------- 
    8992      ! 
     
    97100      l_hst = .FALSE. 
    98101      l_ptr = .FALSE. 
     102      ll_zAimp = .FALSE. 
    99103      IF( ( cdtype =='TRA' .AND. l_trdtra  ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) )       l_trd = .TRUE. 
    100104      IF(   cdtype =='TRA' .AND. ln_diaptr )                                                l_ptr = .TRUE.  
     
    116120      ! 
    117121      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 
    118140      ! 
    119141      DO jn = 1, kjpt            !==  loop over the tracers  ==! 
     
    169191            END DO 
    170192         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 
    171218         !                 
    172219         IF( l_trd .OR. l_hst )  THEN             ! trend diagnostics (contribution of upstream fluxes) 
     
    280327         CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1.,  zwz, 'W',  1. ) 
    281328         ! 
     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 
    282356         !        !==  monotonicity algorithm  ==! 
    283357         ! 
     
    289363            DO jj = 2, jpjm1 
    290364               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          
    298394         ! 
    299395         IF( l_trd .OR. l_hst ) THEN   ! trend diagnostics // heat/salt transport 
     
    318414      END DO                     ! end of tracer loop 
    319415      ! 
     416      IF ( ll_zAimp ) THEN 
     417         DEALLOCATE( zwdia, zwinf, zwsup ) 
     418      ENDIF 
    320419      IF( l_trd .OR. l_hst ) THEN  
    321420         DEALLOCATE( ztrdx, ztrdy, ztrdz ) 
  • NEMO/trunk/src/OCE/stpctl.F90

    r10570 r11407  
    9696            IF( ln_zad_Aimp ) THEN 
    9797               istatus = NF90_DEF_VAR( idrun,   'abs_wi_max', NF90_DOUBLE, (/ idtime /), idw1  ) 
    98                istatus = NF90_DEF_VAR( idrun,       'Cu_max', NF90_DOUBLE, (/ idtime /), idc1  ) 
     98               istatus = NF90_DEF_VAR( idrun,       'Cf_max', NF90_DOUBLE, (/ idtime /), idc1  ) 
    9999            ENDIF 
    100100            istatus = NF90_ENDDEF(idrun) 
     
    123123      IF( ln_zad_Aimp ) THEN 
    124124         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. max 
     125         zmax(9) = MAXVAL(   Cu_adv(:,:,:)   , mask = tmask(:,:,:) == 1._wp ) ! partitioning coeff. max 
    126126      ENDIF 
    127127      ! 
Note: See TracChangeset for help on using the changeset viewer.