New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 10377 – NEMO

Changeset 10377


Ignore:
Timestamp:
2018-12-10T08:45:39+01:00 (6 years ago)
Author:
smasson
Message:

dev_r10164_HPC09_ESIWACE_PREP_MERGE: merge with trunk@10376, see #2133

Location:
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE
Files:
15 edited
1 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_top_cfg

    r9516 r10377  
    9090/ 
    9191!----------------------------------------------------------------------- 
     92&namtrc_snk      !  sedimentation of particles 
     93!----------------------------------------------------------------------- 
     94/ 
     95!----------------------------------------------------------------------- 
    9296&namtrc_dmp      !   passive tracer newtonian damping    
    9397!----------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_top_cfg

    r9516 r10377  
    8989/ 
    9090!----------------------------------------------------------------------- 
     91&namtrc_snk      !  sedimentation of particles 
     92!----------------------------------------------------------------------- 
     93/ 
     94!----------------------------------------------------------------------- 
    9195&namtrc_dmp      !   passive tracer newtonian damping    
    9296!----------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/SHARED/namelist_pisces_ref

    r10368 r10377  
    5151   wsbio2max  =  50.      ! Big particles maximum sinking speed 
    5252   wsbio2scale =  5000.    ! Big particles length scale of sinking 
    53    niter1max  =  1        ! Maximum number of iterations for POC 
    54    niter2max  =  2        ! Maximum number of iterations for GOC 
    5553!                         !  ln_ligand enabled 
    5654   wfep       =  0.2      ! FeP sinking speed  
    5755   ldocp      =  1.E-4    ! Phyto ligand production per unit doc  
    5856   ldocz      =  1.E-4    ! Zoo ligand production per unit doc  
    59    lthet      =  0.5      ! Proportional loss of ligands due to Fe uptake  
     57   lthet      =  1.0      ! Proportional loss of ligands due to Fe uptake  
    6058!                         !  ln_p5z enabled 
    6159   no3rat3    =  0.182    ! N/C ratio in zooplankton 
    … …  
    310308   xlam1        =  0.005  ! scavenging rate of Iron 
    311309   xlamdust     =  150.0  ! Scavenging rate of dust 
    312    ligand       =  0.6E-9 ! Ligands concentration  
     310   ligand       =  0.7E-9 ! Ligands concentration  
    313311   kfep         =  0.01   ! Nanoparticle formation rate constant 
    314312/ 
    … …  
    322320   xsilab    =  0.5       ! Fraction of labile biogenic silica 
    323321   feratb    =  10.E-6    ! Fe/C quota in bacteria 
    324    xkferb    =  2.5E-10   ! Half-saturation constant for bacteria Fe/C 
     322   xkferb    =  3E-10     ! Half-saturation constant for bacteria Fe/C 
    325323!                         ! ln_p5z 
    326324   xremikc   =  0.25      ! remineralization rate of DOC 
    … …  
    372370   ln_ironsed  =  .true.   ! boolean for Fe input from sediments 
    373371   ln_ironice  =  .true.   ! boolean for Fe input from sea ice 
    374    ln_hydrofe  =  .false.  ! boolean for from hydrothermal vents 
     372   ln_hydrofe  =  .true.   ! boolean for from hydrothermal vents 
    375373   sedfeinput  =  2.e-9    ! Coastal release of Iron 
    376374   distcoast   =  5.e3     ! Distance off the coast for Iron from sediments 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/SHARED/namelist_top_ref

    r9526 r10377  
    9595/ 
    9696!----------------------------------------------------------------------- 
     97&namtrc_snk      !  Sedimentation of particles 
     98!----------------------------------------------------------------------- 
     99   nitermax      =  2   !  number of iterations for sedimentation 
     100/ 
     101!----------------------------------------------------------------------- 
    97102&namtrc_dmp      !   passive tracer newtonian damping                   (ln_trcdmp=T) 
    98103!----------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/main/NEMO_manual.bib

    r10124 r10377  
    547547PAGES = {1285--1297}, 
    548548DOI = {10.5194/gmd-8-1285-2015} 
     549} 
     550 
     551@ARTICLE{Breivik_al_JPO2014, 
     552   AUTHOR = {{\O}yvind Breivik and Peter A.E.M. Janssen and Jean-Raymond Bidlot}, 
     553   YEAR = {2014}, 
     554   TITLE = "{Approximate Stokes Drift Profiles in Deep Water}", 
     555   JOURNAL = {JPO}, 
     556   VOLUME = {44}, 
     557   NUMBER = {9}, 
     558   DOI = {10.1175/JPO-D-14-0020.1.}, 
     559   PAGES = {2433--2445, arXiv:1406.5039} 
    549560} 
    550561 
    … …  
    16551666} 
    16561667 
     1668@TECHREPORT{Janssen_al_TM13, 
     1669 author = {P.A.E.M. Janssen and {\O}. Breivik and K. Mogensen  
     1670           and F. Vitart and  M. Balmaseda and J.B. Bidlot and  
     1671           S. Keeley and M. Leut-becher and L. Magnusson and F. Molteni}, 
     1672 title = {Air-Sea Interaction and Surface Waves}, 
     1673 year = {2013}, 
     1674 volume = {712}, 
     1675 institution = {ECMWF}, 
     1676} 
     1677 
    16571678@ARTICLE{Jayne_St_Laurent_GRL01, 
    16581679  author = {S.R. Jayne and L.C. {St. Laurent}}, 
    … …  
    29923013  volume = {359}, 
    29933014  pages = {123--129} 
     3015} 
     3016 
     3017@INCOLLECTION{Stokes_1847, 
     3018  author    = {G.G. Stokes}, 
     3019  title     = {On the theory of oscillatory waves}, 
     3020  booktitle = {Transactions of the Cambridge Philosophy Society}, 
     3021  year      = {1847}, 
     3022  volume    = {8}, 
     3023  Pages     = {441--455} 
    29943024} 
    29953025 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/chap_SBC.tex

    r10368 r10377  
    3232\end{itemize} 
    3333 
    34 Five different ways to provide the first six fields to the ocean are available which are controlled by 
     34Four different ways to provide the first six fields to the ocean are available which are controlled by 
    3535namelist \ngn{namsbc} variables: 
    3636an analytical formulation (\np{ln\_ana}\forcode{ = .true.}), 
    3737a flux formulation (\np{ln\_flx}\forcode{ = .true.}), 
    3838a bulk formulae formulation (CORE (\np{ln\_blk\_core}\forcode{ = .true.}), 
    39 CLIO (\np{ln\_blk\_clio}\forcode{ = .true.}) or 
    40 MFS \footnote { Note that MFS bulk formulae compute fluxes only for the ocean component} 
    41 (\np{ln\_blk\_mfs}\forcode{ = .true.}) bulk formulae) and 
     39CLIO (\np{ln\_blk\_clio}\forcode{ = .true.}) bulk formulae) and 
    4240a coupled or mixed forced/coupled formulation (exchanges with a atmospheric model via the OASIS coupler) 
    4341(\np{ln\_cpl} or \np{ln\_mixcpl}\forcode{ = .true.}).  
    … …  
    7674  the transformation of the solar radiation (if provided as daily mean) into a diurnal cycle 
    7775  (\np{ln\_dm2dc}\forcode{ = .true.}); 
    78   and a neutral drag coefficient can be read from an external wave model (\np{ln\_cdgw}\forcode{ = .true.}).  
     76\item 
     77  a neutral drag coefficient can be read from an external wave model (\np{ln\_cdgw}\forcode{ = .true.}); 
     78\item 
     79  the Stokes drift rom an external wave model can be accounted (\np{ln\_sdw}\forcode{ = .true.});  
     80\item 
     81  the Stokes-Coriolis term can be included (\np{ln\_stcor}\forcode{ = .true.}); 
     82\item 
     83  the surface stress felt by the ocean can be modified by surface waves (\np{ln\_tauwoc}\forcode{ = .true.}). 
    7984\end{itemize} 
    80 The latter option is possible only in case core or mfs bulk formulas are selected. 
    8185 
    8286In this chapter, we first discuss where the surface boundary condition appears in the model equations. 
    … …  
    593597% Bulk formulation 
    594598% ================================================================ 
    595 \section[Bulk formulation {(\textit{sbcblk\{\_core,\_clio,\_mfs\}.F90})}] 
    596          {Bulk formulation {(\protect\mdl{sbcblk\_core}, \protect\mdl{sbcblk\_clio}, \protect\mdl{sbcblk\_mfs})}} 
     599\section[Bulk formulation {(\textit{sbcblk\{\_core,\_clio\}.F90})}] 
     600                        {Bulk formulation {(\protect\mdl{sbcblk\_core}, \protect\mdl{sbcblk\_clio})}} 
    597601\label{sec:SBC_blk} 
    598602 
    … …  
    600604 
    601605The atmospheric fields used depend on the bulk formulae used. 
    602 Three bulk formulations are available: 
    603 the CORE, the CLIO and the MFS bulk formulea. 
     606Two bulk formulations are available: 
     607the CORE and the CLIO bulk formulea. 
    604608The choice is made by setting to true one of the following namelist variable: 
    605 \np{ln\_core} ; \np{ln\_clio} or  \np{ln\_mfs}. 
     609\np{ln\_core} or \np{ln\_clio}. 
    606610 
    607611Note: 
    … …  
    712716the namsbc\_blk\_core or namsbc\_blk\_clio namelist (see \autoref{subsec:SBC_fldread}).  
    713717 
    714 % ------------------------------------------------------------------------------------------------------------- 
    715 %        MFS Bulk formulae 
    716 % ------------------------------------------------------------------------------------------------------------- 
    717 \subsection{MFS formulea (\protect\mdl{sbcblk\_mfs}, \protect\np{ln\_mfs}\forcode{ = .true.})} 
    718 \label{subsec:SBC_blk_mfs} 
    719 %------------------------------------------namsbc_mfs---------------------------------------------------- 
    720 % 
    721 %\nlst{namsbc_mfs} 
    722 %---------------------------------------------------------------------------------------------------------- 
    723  
    724 The MFS (Mediterranean Forecasting System) bulk formulae have been developed by \citet{Castellari_al_JMS1998}.  
    725 They have been designed to handle the ECMWF operational data and are currently in use in 
    726 the MFS operational system \citep{Tonani_al_OS08}, \citep{Oddo_al_OS09}. 
    727 The wind stress computation uses a drag coefficient computed according to \citet{Hellerman_Rosenstein_JPO83}. 
    728 The surface boundary condition for temperature involves the balance between 
    729 surface solar radiation, net long-wave radiation, the latent and sensible heat fluxes. 
    730 Solar radiation is dependent on cloud cover and is computed by means of an astronomical formula \citep{Reed_JPO77}. 
    731 Albedo monthly values are from \citet{Payne_JAS72} as means of the values at $40^{o}N$ and $30^{o}N$ for 
    732 the Atlantic Ocean (hence the same latitudinal band of the Mediterranean Sea). 
    733 The net long-wave radiation flux \citep{Bignami_al_JGR95} is a function of 
    734 air temperature, sea-surface temperature, cloud cover and relative humidity. 
    735 Sensible heat and latent heat fluxes are computed by classical bulk formulae parameterised according to 
    736 \citet{Kondo1975}. 
    737 Details on the bulk formulae used can be found in \citet{Maggiore_al_PCE98} and \citet{Castellari_al_JMS1998}. 
    738  
    739 Options are defined through the \ngn{namsbc\_mfs} namelist variables. 
    740 The required 7 input fields must be provided on the model Grid-T and are: 
    741 \begin{itemize} 
    742 \item          Zonal Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windi}) 
    743 \item          Meridional Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windj}) 
    744 \item          Total Claud Cover (\%)  (\np{sn\_clc}) 
    745 \item          2m Air Temperature ($K$) (\np{sn\_tair}) 
    746 \item          2m Dew Point Temperature ($K$)  (\np{sn\_rhm}) 
    747 \item          Total Precipitation ${Kg} m^{-2} s^{-1}$ (\np{sn\_prec}) 
    748 \item          Mean Sea Level Pressure (${Pa}$) (\np{sn\_msl}) 
    749 \end{itemize} 
    750 % ------------------------------------------------------------------------------------------------------------- 
    751718% ================================================================ 
    752719% Coupled formulation 
    … …  
    12031170since its trajectory data may be spread across multiple files. 
    12041171 
     1172% ------------------------------------------------------------------------------------------------------------- 
     1173%        Interactions with waves (sbcwave.F90, ln_wave) 
     1174% ------------------------------------------------------------------------------------------------------------- 
     1175\section{Interactions with waves (\protect\mdl{sbcwave}, \protect\np{ln\_wave})} 
     1176\label{sec:SBC_wave} 
     1177%------------------------------------------namsbc_wave-------------------------------------------------------- 
     1178 
     1179\nlst{namsbc_wave} 
     1180%------------------------------------------------------------------------------------------------------------- 
     1181 
     1182Ocean waves represent the interface between the ocean and the atmosphere, so NEMO is extended to incorporate  
     1183physical processes related to ocean surface waves, namely the surface stress modified by growth and  
     1184dissipation of the oceanic wave field, the Stokes-Coriolis force and the Stokes drift impact on mass and  
     1185tracer advection; moreover the neutral surface drag coefficient from a wave model can be used to evaluate  
     1186the wind stress. 
     1187 
     1188Physical processes related to ocean surface waves can be accounted by setting the logical variable  
     1189\np{ln\_wave}\forcode{= .true.} in \ngn{namsbc} namelist. In addition, specific flags accounting for  
     1190different porcesses should be activated as explained in the following sections. 
     1191 
     1192Wave fields can be provided either in forced or coupled mode: 
     1193\begin{description} 
     1194\item[forced mode]: wave fields should be defined through the \ngn{namsbc\_wave} namelist  
     1195for external data names, locations, frequency, interpolation and all the miscellanous options allowed by  
     1196Input Data generic Interface (see \autoref{sec:SBC_input}).  
     1197\item[coupled mode]: NEMO and an external wave model can be coupled by setting \np{ln\_cpl} \forcode{= .true.}  
     1198in \ngn{namsbc} namelist and filling the \ngn{namsbc\_cpl} namelist. 
     1199\end{description} 
     1200 
     1201 
     1202% ================================================================ 
     1203% Neutral drag coefficient from wave model (ln_cdgw) 
     1204 
     1205% ================================================================ 
     1206\subsection{Neutral drag coefficient from wave model (\protect\np{ln\_cdgw})} 
     1207\label{subsec:SBC_wave_cdgw} 
     1208 
     1209The neutral surface drag coefficient provided from an external data source ($i.e.$ a wave 
     1210model),  
     1211can be used by setting the logical variable \np{ln\_cdgw} \forcode{= .true.} in \ngn{namsbc} namelist.  
     1212Then using the routine \rou{turb\_ncar} and starting from the neutral drag coefficent provided,  
     1213the drag coefficient is computed according to the stable/unstable conditions of the  
     1214air-sea interface following \citet{Large_Yeager_Rep04}.  
     1215 
     1216 
     1217% ================================================================ 
     1218% 3D Stokes Drift (ln_sdw, nn_sdrift) 
     1219% ================================================================ 
     1220\subsection{3D Stokes Drift (\protect\np{ln\_sdw, nn\_sdrift})} 
     1221\label{subsec:SBC_wave_sdw} 
     1222 
     1223The Stokes drift is a wave driven mechanism of mass and momentum transport \citep{Stokes_1847}.  
     1224It is defined as the difference between the average velocity of a fluid parcel (Lagrangian velocity)  
     1225and the current measured at a fixed point (Eulerian velocity).  
     1226As waves travel, the water particles that make up the waves travel in orbital motions but  
     1227without a closed path. Their movement is enhanced at the top of the orbit and slowed slightly  
     1228at the bottom so the result is a net forward motion of water particles, referred to as the Stokes drift.  
     1229An accurate evaluation of the Stokes drift and the inclusion of related processes may lead to improved  
     1230representation of surface physics in ocean general circulation models. 
     1231The Stokes drift velocity $\mathbf{U}_{st}$ in deep water can be computed from the wave spectrum and may be written as:  
     1232 
     1233\begin{equation} \label{eq:sbc_wave_sdw} 
     1234\mathbf{U}_{st} = \frac{16{\pi^3}} {g}  
     1235                \int_0^\infty \int_{-\pi}^{\pi} (cos{\theta},sin{\theta}) {f^3} 
     1236                \mathrm{S}(f,\theta) \mathrm{e}^{2kz}\,\mathrm{d}\theta {d}f 
     1237\end{equation} 
     1238 
     1239where: ${\theta}$ is the wave direction, $f$ is the wave intrinsic frequency,  
     1240$\mathrm{S}($f$,\theta)$ is the 2D frequency-direction spectrum,  
     1241$k$ is the mean wavenumber defined as:  
     1242$k=\frac{2\pi}{\lambda}$ (being $\lambda$ the wavelength). \\ 
     1243 
     1244In order to evaluate the Stokes drift in a realistic ocean wave field the wave spectral shape is required  
     1245and its computation quickly becomes expensive as the 2D spectrum must be integrated for each vertical level.  
     1246To simplify, it is customary to use approximations to the full Stokes profile. 
     1247Three possible parameterizations for the calculation for the approximate Stokes drift velocity profile  
     1248are included in the code through the \np{nn\_sdrift} parameter once provided the surface Stokes drift  
     1249$\mathbf{U}_{st |_{z=0}}$ which is evaluated by an external wave model that accurately reproduces the wave spectra  
     1250and makes possible the estimation of the surface Stokes drift for random directional waves in  
     1251realistic wave conditions: 
     1252 
     1253\begin{description} 
     1254\item[\np{nn\_sdrift} = 0]: exponential integral profile parameterization proposed by  
     1255\citet{Breivik_al_JPO2014}: 
     1256 
     1257\begin{equation} \label{eq:sbc_wave_sdw_0a} 
     1258\mathbf{U}_{st} \cong \mathbf{U}_{st |_{z=0}} \frac{\mathrm{e}^{-2k_ez}} {1-8k_ez}  
     1259\end{equation} 
     1260 
     1261where $k_e$ is the effective wave number which depends on the Stokes transport $T_{st}$ defined as follows: 
     1262 
     1263\begin{equation} \label{eq:sbc_wave_sdw_0b} 
     1264k_e = \frac{|\mathbf{U}_{\left.st\right|_{z=0}}|} {|T_{st}|}  
     1265\quad \text{and }\ 
     1266T_{st} = \frac{1}{16} \bar{\omega} H_s^2  
     1267\end{equation} 
     1268 
     1269where $H_s$ is the significant wave height and $\omega$ is the wave frequency. 
     1270 
     1271\item[\np{nn\_sdrift} = 1]: velocity profile based on the Phillips spectrum which is considered to be a  
     1272reasonable estimate of the part of the spectrum most contributing to the Stokes drift velocity near the surface 
     1273\citep{Breivik_al_OM2016}: 
     1274 
     1275\begin{equation} \label{eq:sbc_wave_sdw_1} 
     1276\mathbf{U}_{st} \cong \mathbf{U}_{st |_{z=0}} \Big[exp(2k_pz)-\beta \sqrt{-2 \pi k_pz}  
     1277\textit{ erf } \Big(\sqrt{-2 k_pz}\Big)\Big] 
     1278\end{equation} 
     1279 
     1280where $erf$ is the complementary error function and $k_p$ is the peak wavenumber. 
     1281 
     1282\item[\np{nn\_sdrift} = 2]: velocity profile based on the Phillips spectrum as for \np{nn\_sdrift} = 1  
     1283but using the wave frequency from a wave model. 
     1284 
     1285\end{description} 
     1286 
     1287The Stokes drift enters the wave-averaged momentum equation, as well as the tracer advection equations  
     1288and its effect on the evolution of the sea-surface height ${\eta}$ is considered as follows:  
     1289 
     1290\begin{equation} \label{eq:sbc_wave_eta_sdw} 
     1291\frac{\partial{\eta}}{\partial{t}} =  
     1292-\nabla_h \int_{-H}^{\eta} (\mathbf{U} + \mathbf{U}_{st}) dz  
     1293\end{equation} 
     1294 
     1295The tracer advection equation is also modified in order for Eulerian ocean models to properly account  
     1296for unresolved wave effect. The divergence of the wave tracer flux equals the mean tracer advection  
     1297that is induced by the three-dimensional Stokes velocity.  
     1298The advective equation for a tracer $c$ combining the effects of the mean current and sea surface waves  
     1299can be formulated as follows:  
     1300 
     1301\begin{equation} \label{eq:sbc_wave_tra_sdw} 
     1302\frac{\partial{c}}{\partial{t}} =  
     1303- (\mathbf{U} + \mathbf{U}_{st}) \cdot \nabla{c} 
     1304\end{equation} 
     1305 
     1306 
     1307% ================================================================ 
     1308% Stokes-Coriolis term (ln_stcor) 
     1309% ================================================================ 
     1310\subsection{Stokes-Coriolis term (\protect\np{ln\_stcor})} 
     1311\label{subsec:SBC_wave_stcor} 
     1312 
     1313In a rotating ocean, waves exert a wave-induced stress on the mean ocean circulation which results  
     1314in a force equal to $\mathbf{U}_{st}$×$f$, where $f$ is the Coriolis parameter.  
     1315This additional force may have impact on the Ekman turning of the surface current.  
     1316In order to include this term, once evaluated the Stokes drift (using one of the 3 possible  
     1317approximations described in \autoref{subsec:SBC_wave_sdw}),  
     1318\np{ln\_stcor}\forcode{ = .true.} has to be set. 
     1319 
     1320 
     1321% ================================================================ 
     1322% Waves modified stress (ln_tauwoc, ln_tauw) 
     1323% ================================================================ 
     1324\subsection{Wave modified sress (\protect\np{ln\_tauwoc, ln\_tauw})}  
     1325\label{subsec:SBC_wave_tauw} 
     1326 
     1327The surface stress felt by the ocean is the atmospheric stress minus the net stress going  
     1328into the waves \citep{Janssen_al_TM13}. Therefore, when waves are growing, momentum and energy is spent and is not  
     1329available for forcing the mean circulation, while in the opposite case of a decaying sea  
     1330state more momentum is available for forcing the ocean.  
     1331Only when the sea state is in equilibrium the ocean is forced by the atmospheric stress,  
     1332but in practice an equilibrium sea state is a fairly rare event.  
     1333So the atmospheric stress felt by the ocean circulation $\tau_{oc,a}$ can be expressed as:  
     1334 
     1335\begin{equation} \label{eq:sbc_wave_tauoc} 
     1336\tau_{oc,a} = \tau_a - \tau_w 
     1337\end{equation} 
     1338 
     1339where $\tau_a$ is the atmospheric surface stress; 
     1340$\tau_w$ is the atmospheric stress going into the waves defined as: 
     1341 
     1342\begin{equation} \label{eq:sbc_wave_tauw} 
     1343\tau_w = \rho g \int {\frac{dk}{c_p} (S_{in}+S_{nl}+S_{diss})} 
     1344\end{equation} 
     1345 
     1346where: $c_p$ is the phase speed of the gravity waves, 
     1347$S_{in}$, $S_{nl}$ and $S_{diss}$ are three source terms that represent  
     1348the physics of ocean waves. The first one, $S_{in}$, describes the generation  
     1349of ocean waves by wind and therefore represents the momentum and energy transfer  
     1350from air to ocean waves; the second term $S_{nl}$ denotes  
     1351the nonlinear transfer by resonant four-wave interactions; while the third term $S_{diss}$  
     1352describes the dissipation of waves by processes such as white-capping, large scale breaking  
     1353eddy-induced damping. 
     1354 
     1355The wave stress derived from an external wave model can be provided either through the normalized  
     1356wave stress into the ocean by setting \np{ln\_tauwoc}\forcode{ = .true.}, or through the zonal and  
     1357meridional stress components by setting \np{ln\_tauw}\forcode{ = .true.}. 
     1358 
    12051359 
    12061360% ================================================================ 
    … …  
    14211575\end{description} 
    14221576 
    1423 % ------------------------------------------------------------------------------------------------------------- 
    1424 %        Neutral Drag Coefficient from external wave model 
    1425 % ------------------------------------------------------------------------------------------------------------- 
    1426 \subsection[Neutral drag coeff. from external wave model (\protect\mdl{sbcwave})] 
    1427             {Neutral drag coefficient from external wave model (\protect\mdl{sbcwave})} 
    1428 \label{subsec:SBC_wave} 
    1429 %------------------------------------------namwave---------------------------------------------------- 
    1430  
    1431 \nlst{namsbc_wave} 
    1432 %------------------------------------------------------------------------------------------------------------- 
    1433  
    1434 In order to read a neutral drag coefficient, from an external data source ($i.e.$ a wave model), 
    1435 the logical variable \np{ln\_cdgw} in \ngn{namsbc} namelist must be set to \forcode{.true.}. 
    1436 The \mdl{sbcwave} module containing the routine \np{sbc\_wave} reads the namelist \ngn{namsbc\_wave} 
    1437 (for external data names, locations, frequency, interpolation and all the miscellanous options allowed by 
    1438 Input Data generic Interface see \autoref{sec:SBC_input}) and 
    1439 a 2D field of neutral drag coefficient. 
    1440 Then using the routine TURB\_CORE\_1Z or TURB\_CORE\_2Z, and starting from the neutral drag coefficent provided,  
    1441 the drag coefficient is computed according to stable/unstable conditions of the air-sea interface following 
    1442 \citet{Large_Yeager_Rep04}. 
    14431577 
    14441578 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/introduction.tex

    r10368 r10377  
    311311\item new definition of configurations ; 
    312312\item bulk formulation ; 
     313\item NEMO-wave large scale interactions ; 
    313314\item ... ;  
    314315\end{enumerate} 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/IOM/iom_nf90.F90

    r10358 r10377  
    691691      INTEGER, DIMENSION(4) :: idimid               ! dimensions id 
    692692      CHARACTER(LEN=256)    :: clinfo               ! info character 
    693       CHARACTER(LEN= 12), DIMENSION(4) :: cltmp     ! temporary character 
     693      CHARACTER(LEN= 12), DIMENSION(5) :: cltmp     ! temporary character 
    694694      INTEGER               :: if90id               ! nf90 file identifier 
    695695      INTEGER               :: idmy                 ! dummy variable 
    … …  
    716716         ENDIF 
    717717         ! define the dimension variables if it is not already done 
    718          IF(iom_file(kiomid)%nlev == jpk ) THEN 
    719           cltmp = (/ 'nav_lon     ', 'nav_lat     ', 'nav_lev     ', 'time_counter' /) 
    720          ELSE 
    721           cltmp = (/ 'nav_lon     ', 'nav_lat     ', 'numcat      ', 'time_counter' /) 
    722          ENDIF 
     718         cltmp = (/ 'nav_lon', 'nav_lat', 'nav_lev', 'time_counter', 'numcat' /) 
    723719         CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(1)), NF90_FLOAT , (/ 1, 2 /), iom_file(kiomid)%nvid(1) ), clinfo) 
    724720         CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(2)), NF90_FLOAT , (/ 1, 2 /), iom_file(kiomid)%nvid(2) ), clinfo) 
    … …  
    728724         iom_file(kiomid)%nvars       = 4 
    729725         iom_file(kiomid)%luld(1:4)   = (/ .FALSE., .FALSE., .FALSE., .TRUE. /) 
    730          iom_file(kiomid)%cn_var(1:4) = cltmp 
    731          iom_file(kiomid)%ndims(1:4)  = (/ 2, 2, 1, 1 /)   
     726         iom_file(kiomid)%cn_var(1:4) = cltmp(1:4) 
     727         iom_file(kiomid)%ndims(1:4)  = (/ 2, 2, 1, 1 /) 
     728         IF( NF90_INQ_DIMID( if90id, 'numcat', idmy ) == nf90_noerr ) THEN   ! add a 5th variable corresponding to the 5th dimension 
     729            CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(5)), NF90_FLOAT , (/ 5 /), iom_file(kiomid)%nvid(5) ), clinfo) 
     730            iom_file(kiomid)%nvars     = 5 
     731            iom_file(kiomid)%luld(5)   = .FALSE. 
     732            iom_file(kiomid)%cn_var(5) = cltmp(5) 
     733            iom_file(kiomid)%ndims(5)  = 1 
     734         ENDIF 
    732735         ! trick: defined to 0 to say that dimension variables are defined but not yet written 
    733736         iom_file(kiomid)%dimsz(1, 1)  = 0    
    … …  
    841844               CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lat'     , idmy )         , clinfo ) 
    842845               CALL iom_nf90_check( NF90_PUT_VAR  ( if90id, idmy, gphit(ix1:ix2, iy1:iy2) ), clinfo ) 
    843                IF(iom_file(kiomid)%nlev == jpk ) THEN  
    844                   !NEMO 
    845                   CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lev'     , idmy ), clinfo ) 
    846                   CALL iom_nf90_check( NF90_PUT_VAR  ( if90id, idmy, gdept_1d       ), clinfo ) 
    847                ELSE 
    848                   CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'numcat'     , idmy ), clinfo) 
     846               CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lev'     , idmy ), clinfo ) 
     847               CALL iom_nf90_check( NF90_PUT_VAR  ( if90id, idmy, gdept_1d       ), clinfo ) 
     848               IF( NF90_INQ_VARID( if90id, 'numcat', idmy ) == nf90_noerr ) THEN 
    849849                  CALL iom_nf90_check( NF90_PUT_VAR  ( if90id, idmy, (/ (idlv, idlv = 1,iom_file(kiomid)%nlev) /)), clinfo ) 
    850850               ENDIF 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/geo2ocean.F90

    r10297 r10377  
    1010   !!            3.7  !  11-2015  (G. Madec)  remove the unused repere and repcmo routines 
    1111   !!---------------------------------------------------------------------- 
    12 #if defined key_agrif 
    13 !clem: these lines do not seem necessary anymore 
    14 !!DIR$ OPTIMIZE (-O 1)  ! cray formulation 
    15 # if defined __INTEL_COMPILER 
    16 !acc: still breaks on at least one Ivybridge cluster with ifort 17.0.4 without this directive 
    17 !DIR$ OPTIMIZE:1        ! intel formulation 
    18 # endif 
    19 #endif 
    2012   !!---------------------------------------------------------------------- 
    2113   !!   rot_rep       : Rotate the Repere: geographic grid <==> stretched coordinates grid 
    … …  
    8173         IF(lwp) WRITE(numout,*) ' ~~~~~~~~    ' 
    8274         ! 
    83          CALL angle       ! initialization of the transformation 
     75         CALL angle( glamt, gphit, glamu, gphiu, glamv, gphiv, glamf, gphif )       ! initialization of the transformation 
    8476         lmust_init = .FALSE. 
    8577      ENDIF 
    … …  
    126118 
    127119 
    128    SUBROUTINE angle 
     120   SUBROUTINE angle( plamt, pphit, plamu, pphiu, plamv, pphiv, plamf, pphif ) 
    129121      !!---------------------------------------------------------------------- 
    130122      !!                  ***  ROUTINE angle  *** 
    … …  
    138130      !! ** Action  : - gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf 
    139131      !!---------------------------------------------------------------------- 
     132      ! WARNING: for an unexplained reason, we need to pass all glam, gphi arrays as input parameters in 
     133      !          order to get AGRIF working with -03 compilation option 
     134      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: plamt, pphit, plamu, pphiu, plamv, pphiv, plamf, pphif   
     135      ! 
    140136      INTEGER  ::   ji, jj   ! dummy loop indices 
    141137      INTEGER  ::   ierr     ! local integer 
    … …  
    167163         DO ji = fs_2, jpi   ! vector opt. 
    168164            !                   
    169             zlam = glamt(ji,jj)     ! north pole direction & modulous (at t-point) 
    170             zphi = gphit(ji,jj) 
     165            zlam = plamt(ji,jj)     ! north pole direction & modulous (at t-point) 
     166            zphi = pphit(ji,jj) 
    171167            zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    172168            zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    173169            znnpt = zxnpt*zxnpt + zynpt*zynpt 
    174170            ! 
    175             zlam = glamu(ji,jj)     ! north pole direction & modulous (at u-point) 
    176             zphi = gphiu(ji,jj) 
     171            zlam = plamu(ji,jj)     ! north pole direction & modulous (at u-point) 
     172            zphi = pphiu(ji,jj) 
    177173            zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    178174            zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    179175            znnpu = zxnpu*zxnpu + zynpu*zynpu 
    180176            ! 
    181             zlam = glamv(ji,jj)     ! north pole direction & modulous (at v-point) 
    182             zphi = gphiv(ji,jj) 
     177            zlam = plamv(ji,jj)     ! north pole direction & modulous (at v-point) 
     178            zphi = pphiv(ji,jj) 
    183179            zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    184180            zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    185181            znnpv = zxnpv*zxnpv + zynpv*zynpv 
    186182            ! 
    187             zlam = glamf(ji,jj)     ! north pole direction & modulous (at f-point) 
    188             zphi = gphif(ji,jj) 
     183            zlam = plamf(ji,jj)     ! north pole direction & modulous (at f-point) 
     184            zphi = pphif(ji,jj) 
    189185            zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    190186            zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. ) 
    191187            znnpf = zxnpf*zxnpf + zynpf*zynpf 
    192188            ! 
    193             zlam = glamv(ji,jj  )   ! j-direction: v-point segment direction (around t-point) 
    194             zphi = gphiv(ji,jj  ) 
    195             zlan = glamv(ji,jj-1) 
    196             zphh = gphiv(ji,jj-1) 
     189            zlam = plamv(ji,jj  )   ! j-direction: v-point segment direction (around t-point) 
     190            zphi = pphiv(ji,jj  ) 
     191            zlan = plamv(ji,jj-1) 
     192            zphh = pphiv(ji,jj-1) 
    197193            zxvvt =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   & 
    198194               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 
    … …  
    202198            znvvt = MAX( znvvt, 1.e-14 ) 
    203199            ! 
    204             zlam = glamf(ji,jj  )   ! j-direction: f-point segment direction (around u-point) 
    205             zphi = gphif(ji,jj  ) 
    206             zlan = glamf(ji,jj-1) 
    207             zphh = gphif(ji,jj-1) 
     200            zlam = plamf(ji,jj  )   ! j-direction: f-point segment direction (around u-point) 
     201            zphi = pphif(ji,jj  ) 
     202            zlan = plamf(ji,jj-1) 
     203            zphh = pphif(ji,jj-1) 
    208204            zxffu =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   & 
    209205               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 
    … …  
    213209            znffu = MAX( znffu, 1.e-14 ) 
    214210            ! 
    215             zlam = glamf(ji  ,jj)   ! i-direction: f-point segment direction (around v-point) 
    216             zphi = gphif(ji  ,jj) 
    217             zlan = glamf(ji-1,jj) 
    218             zphh = gphif(ji-1,jj) 
     211            zlam = plamf(ji  ,jj)   ! i-direction: f-point segment direction (around v-point) 
     212            zphi = pphif(ji  ,jj) 
     213            zlan = plamf(ji-1,jj) 
     214            zphh = pphif(ji-1,jj) 
    219215            zxffv =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   & 
    220216               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 
    … …  
    224220            znffv = MAX( znffv, 1.e-14 ) 
    225221            ! 
    226             zlam = glamu(ji,jj+1)   ! j-direction: u-point segment direction (around f-point) 
    227             zphi = gphiu(ji,jj+1) 
    228             zlan = glamu(ji,jj  ) 
    229             zphh = gphiu(ji,jj  ) 
     222            zlam = plamu(ji,jj+1)   ! j-direction: u-point segment direction (around f-point) 
     223            zphi = pphiu(ji,jj+1) 
     224            zlan = plamu(ji,jj  ) 
     225            zphh = pphiu(ji,jj  ) 
    230226            zxuuf =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   & 
    231227               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. ) 
    … …  
    257253      DO jj = 2, jpjm1 
    258254         DO ji = fs_2, jpi   ! vector opt. 
    259             IF( MOD( ABS( glamv(ji,jj) - glamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN 
     255            IF( MOD( ABS( plamv(ji,jj) - plamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN 
    260256               gsint(ji,jj) = 0. 
    261257               gcost(ji,jj) = 1. 
    262258            ENDIF 
    263             IF( MOD( ABS( glamf(ji,jj) - glamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN 
     259            IF( MOD( ABS( plamf(ji,jj) - plamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN 
    264260               gsinu(ji,jj) = 0. 
    265261               gcosu(ji,jj) = 1. 
    266262            ENDIF 
    267             IF(      ABS( gphif(ji,jj) - gphif(ji-1,jj) )         < 1.e-8 ) THEN 
     263            IF(      ABS( pphif(ji,jj) - pphif(ji-1,jj) )         < 1.e-8 ) THEN 
    268264               gsinv(ji,jj) = 0. 
    269265               gcosv(ji,jj) = 1. 
    270266            ENDIF 
    271             IF( MOD( ABS( glamu(ji,jj) - glamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN 
     267            IF( MOD( ABS( plamu(ji,jj) - plamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN 
    272268               gsinf(ji,jj) = 0. 
    273269               gcosf(ji,jj) = 1. 
    … …  
    457453         IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched' 
    458454         IF(lwp) WRITE(numout,*) ' ~~~~~~~   coordinate transformation' 
    459          CALL angle       ! initialization of the transformation 
     455         CALL angle( glamt, gphit, glamu, gphiu, glamv, gphiv, glamf, gphif )       ! initialization of the transformation 
    460456         lmust_init = .FALSE. 
    461457      ENDIF 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zmeso.F90

    r10368 r10377  
    252252         DEALLOCATE( zw3d ) 
    253253      ENDIF 
     254      ! 
     255      IF (ln_ligand)  DEALLOCATE( zz2ligprod ) 
    254256      ! 
    255257      IF(ln_ctl)   THEN  ! print mean trends (used for debugging) 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zmicro.F90

    r10368 r10377  
    205205      ENDIF 
    206206      ! 
     207      IF (ln_ligand)  DEALLOCATE( zzligprod ) 
     208      ! 
    207209      IF(ln_ctl) THEN      ! print mean trends (used for debugging) 
    208210         WRITE(charout, FMT="('micro')") 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsink.F90

    r10368 r10377  
    1616   USE trc             !  passive tracers common variables  
    1717   USE sms_pisces      !  PISCES Source Minus Sink variables 
     18   USE trcsink         !  General routine to compute sedimentation 
    1819   USE prtctl_trc      !  print control for debugging 
    1920   USE iom             !  I/O manager 
    … …  
    5960      !!--------------------------------------------------------------------- 
    6061      INTEGER, INTENT(in) :: kt, knt 
    61       INTEGER  ::   ji, jj, jk, jit 
    62       INTEGER  ::   iiter1, iiter2 
    63       REAL(wp) ::   zagg1, zagg2, zagg3, zagg4 
    64       REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3 
    65       REAL(wp) ::   zfact, zwsmax, zmax 
     62      INTEGER  ::   ji, jj, jk 
    6663      CHARACTER (len=25) :: charout 
     64      REAL(wp) :: zmax, zfact 
    6765      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d 
    6866      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zw2d 
    … …  
    7068      ! 
    7169      IF( ln_timing )   CALL timing_start('p4z_sink') 
    72  
    7370 
    7471      ! Initialization of some global variables 
    … …  
    9794 
    9895      ! 
    99       ! OA This is (I hope) a temporary solution for the problem that may  
    100       ! OA arise in specific situation where the CFL criterion is broken  
    101       ! OA for vertical sedimentation of particles. To avoid this, a time 
    102       ! OA splitting algorithm has been coded. A specific maximum 
    103       ! OA iteration number is provided and may be specified in the namelist  
    104       ! OA This is to avoid very large iteration number when explicit free 
    105       ! OA surface is used (for instance). When niter?max is set to 1,  
    106       ! OA this computation is skipped. The crude old threshold method is  
    107       ! OA then applied. This also happens when niter exceeds nitermax. 
    108       IF( MAX( niter1max, niter2max ) == 1 ) THEN 
    109         iiter1 = 1 
    110         iiter2 = 1 
    111       ELSE 
    112         iiter1 = 1 
    113         iiter2 = 1 
    114         DO jk = 1, jpkm1 
    115           DO jj = 1, jpj 
    116              DO ji = 1, jpi 
    117                 IF( tmask(ji,jj,jk) == 1) THEN 
    118                    zwsmax =  0.5 * e3t_n(ji,jj,jk) / xstep 
    119                    iiter1 =  MAX( iiter1, INT( wsbio3(ji,jj,jk) / zwsmax ) ) 
    120                    iiter2 =  MAX( iiter2, INT( wsbio4(ji,jj,jk) / zwsmax ) ) 
    121                 ENDIF 
    122              END DO 
    123           END DO 
    124         END DO 
    125         IF( lk_mpp ) THEN 
    126            CALL mpp_max( 'p4zsink', iiter1 ) 
    127            CALL mpp_max( 'p4zsink', iiter2 ) 
    128         ENDIF 
    129         iiter1 = MIN( iiter1, niter1max ) 
    130         iiter2 = MIN( iiter2, niter2max ) 
    131       ENDIF 
    132  
    133       DO jk = 1,jpkm1 
    134          DO jj = 1, jpj 
    135             DO ji = 1, jpi 
    136                IF( tmask(ji,jj,jk) == 1 ) THEN 
    137                  zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep 
    138                  wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * REAL( iiter1, wp ) ) 
    139                  wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * REAL( iiter2, wp ) ) 
    140                ENDIF 
    141             END DO 
    142          END DO 
    143       END DO 
    144  
    14596      !  Initializa to zero all the sinking arrays  
    14697      !   ----------------------------------------- 
    … …  
    154105      !   Compute the sedimentation term using p4zsink2 for all the sinking particles 
    155106      !   ----------------------------------------------------- 
    156       DO jit = 1, iiter1 
    157         CALL p4z_sink2( wsbio3, sinking , jppoc, iiter1 ) 
    158         CALL p4z_sink2( wsbio3, sinkfer , jpsfe, iiter1 ) 
    159       END DO 
    160  
    161       DO jit = 1, iiter2 
    162         CALL p4z_sink2( wsbio4, sinking2, jpgoc, iiter2 ) 
    163         CALL p4z_sink2( wsbio4, sinkfer2, jpbfe, iiter2 ) 
    164         CALL p4z_sink2( wsbio4, sinksil , jpgsi, iiter2 ) 
    165         CALL p4z_sink2( wsbio4, sinkcal , jpcal, iiter2 ) 
    166       END DO 
     107      CALL trc_sink( kt, wsbio3, sinking , jppoc, rfact2 ) 
     108      CALL trc_sink( kt, wsbio3, sinkfer , jpsfe, rfact2 ) 
     109      CALL trc_sink( kt, wsbio4, sinking2, jpgoc, rfact2 ) 
     110      CALL trc_sink( kt, wsbio4, sinkfer2, jpbfe, rfact2 ) 
     111      CALL trc_sink( kt, wsbio4, sinksil , jpgsi, rfact2 ) 
     112      CALL trc_sink( kt, wsbio4, sinkcal , jpcal, rfact2 ) 
    167113 
    168114      IF( ln_p5z ) THEN 
    … …  
    174120         !   Compute the sedimentation term using p4zsink2 for all the sinking particles 
    175121         !   ----------------------------------------------------- 
    176          DO jit = 1, iiter1 
    177            CALL p4z_sink2( wsbio3, sinkingn , jppon, iiter1 ) 
    178            CALL p4z_sink2( wsbio3, sinkingp , jppop, iiter1 ) 
    179          END DO 
    180  
    181          DO jit = 1, iiter2 
    182            CALL p4z_sink2( wsbio4, sinking2n, jpgon, iiter2 ) 
    183            CALL p4z_sink2( wsbio4, sinking2p, jpgop, iiter2 ) 
    184          END DO 
     122         CALL trc_sink( kt, wsbio3, sinkingn , jppon, rfact2 ) 
     123         CALL trc_sink( kt, wsbio3, sinkingp , jppop, rfact2 ) 
     124         CALL trc_sink( kt, wsbio4, sinking2n, jpgon, rfact2 ) 
     125         CALL trc_sink( kt, wsbio4, sinking2p, jpgop, rfact2 ) 
    185126      ENDIF 
    186127 
    187128      IF( ln_ligand ) THEN 
    188129         wsfep (:,:,:) = wfep 
    189          DO jk = 1,jpkm1 
    190             DO jj = 1, jpj 
    191                DO ji = 1, jpi 
    192                   IF( tmask(ji,jj,jk) == 1 ) THEN 
    193                     zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep 
    194                     wsfep(ji,jj,jk) = MIN( wsfep(ji,jj,jk), zwsmax * REAL( iiter1, wp ) ) 
    195                   ENDIF 
    196                END DO 
    197             END DO 
    198          END DO 
    199130         ! 
    200131         sinkfep(:,:,:) = 0.e0 
    201          DO jit = 1, iiter1 
    202            CALL p4z_sink2( wsfep, sinkfep , jpfep, iiter1 ) 
    203          END DO 
     132         CALL trc_sink( kt, wsfep, sinkfep , jpfep, rfact2 ) 
    204133      ENDIF 
    205134 
    … …  
    281210   END SUBROUTINE p4z_sink_init 
    282211 
    283  
    284    SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter ) 
    285       !!--------------------------------------------------------------------- 
    286       !!                     ***  ROUTINE p4z_sink2  *** 
    287       !! 
    288       !! ** Purpose :   Compute the sedimentation terms for the various sinking 
    289       !!     particles. The scheme used to compute the trends is based 
    290       !!     on MUSCL. 
    291       !! 
    292       !! ** Method  : - this ROUTINE compute not exactly the advection but the 
    293       !!      transport term, i.e.  div(u*tra). 
    294       !!--------------------------------------------------------------------- 
    295       INTEGER , INTENT(in   )                         ::   jp_tra    ! tracer index index       
    296       INTEGER , INTENT(in   )                         ::   kiter     ! number of iterations for time-splitting  
    297       REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed 
    298       REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe 
    299       ! 
    300       INTEGER  ::   ji, jj, jk, jn 
    301       REAL(wp) ::   zigma,zew,zign, zflx, zstep 
    302       REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz, zwsink2, ztrb  
    303       !!--------------------------------------------------------------------- 
    304       ! 
    305       IF( ln_timing )   CALL timing_start('p4z_sink2') 
    306       ! 
    307       zstep = rfact2 / REAL( kiter, wp ) / 2. 
    308  
    309       ztraz(:,:,:) = 0.e0 
    310       zakz (:,:,:) = 0.e0 
    311       ztrb (:,:,:) = trb(:,:,:,jp_tra) 
    312  
    313       DO jk = 1, jpkm1 
    314          zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1)  
    315       END DO 
    316       zwsink2(:,:,1) = 0.e0 
    317  
    318  
    319       ! Vertical advective flux 
    320       DO jn = 1, 2 
    321          !  first guess of the slopes interior values 
    322          DO jk = 2, jpkm1 
    323             ztraz(:,:,jk) = ( trb(:,:,jk-1,jp_tra) - trb(:,:,jk,jp_tra) ) * tmask(:,:,jk) 
    324          END DO 
    325          ztraz(:,:,1  ) = 0.0 
    326          ztraz(:,:,jpk) = 0.0 
    327  
    328          ! slopes 
    329          DO jk = 2, jpkm1 
    330             DO jj = 1,jpj 
    331                DO ji = 1, jpi 
    332                   zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) ) 
    333                   zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign 
    334                END DO 
    335             END DO 
    336          END DO 
    337           
    338          ! Slopes limitation 
    339          DO jk = 2, jpkm1 
    340             DO jj = 1, jpj 
    341                DO ji = 1, jpi 
    342                   zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) *        & 
    343                      &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) ) 
    344                END DO 
    345             END DO 
    346          END DO 
    347           
    348          ! vertical advective flux 
    349          DO jk = 1, jpkm1 
    350             DO jj = 1, jpj       
    351                DO ji = 1, jpi     
    352                   zigma = zwsink2(ji,jj,jk+1) * zstep / e3w_n(ji,jj,jk+1) 
    353                   zew   = zwsink2(ji,jj,jk+1) 
    354                   psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep 
    355                END DO 
    356             END DO 
    357          END DO 
    358          ! 
    359          ! Boundary conditions 
    360          psinkflx(:,:,1  ) = 0.e0 
    361          psinkflx(:,:,jpk) = 0.e0 
    362           
    363          DO jk=1,jpkm1 
    364             DO jj = 1,jpj 
    365                DO ji = 1, jpi 
    366                   zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk) 
    367                   trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx 
    368                END DO 
    369             END DO 
    370          END DO 
    371  
    372       ENDDO 
    373  
    374       DO jk = 1,jpkm1 
    375          DO jj = 1,jpj 
    376             DO ji = 1, jpi 
    377                zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk) 
    378                ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx 
    379             END DO 
    380          END DO 
    381       END DO 
    382  
    383       trb(:,:,:,jp_tra) = ztrb(:,:,:) 
    384       psinkflx(:,:,:)   = 2. * psinkflx(:,:,:) 
    385       ! 
    386       IF( ln_timing )  CALL timing_stop('p4z_sink2') 
    387       ! 
    388    END SUBROUTINE p4z_sink2 
    389  
    390  
    391212   INTEGER FUNCTION p4z_sink_alloc() 
    392213      !!---------------------------------------------------------------------- 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsms.F90

    r10345 r10377  
    7474            CALL p4z_che                              ! initialize the chemical constants 
    7575            CALL ahini_for_at(hi)   !  set PH at kt=nit000 
    76             t_oce_co2_flx_cum = 0._wp 
    7776        ELSE 
    7877            CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields 
    … …  
    189188      !! 
    190189      NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale,    & 
    191          &                   niter1max, niter2max, wfep, ldocp, ldocz, lthet,  & 
     190         &                   wfep, ldocp, ldocz, lthet,  & 
    192191         &                   no3rat3, po4rat3 
    193192         ! 
    … …  
    223222         WRITE(numout,*) '      Big particles maximum sinking speed       wsbio2max   =', wsbio2max 
    224223         WRITE(numout,*) '      Big particles sinking speed length scale  wsbio2scale =', wsbio2scale 
    225          WRITE(numout,*) '      Maximum number of iterations for POC      niter1max   =', niter1max 
    226          WRITE(numout,*) '      Maximum number of iterations for GOC      niter2max   =', niter2max 
    227224         IF( ln_ligand ) THEN 
    228225            WRITE(numout,*) '      FeP sinking speed                              wfep   =', wfep 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/sms_pisces.F90

    r10368 r10377  
    3636 
    3737   !!*  Biological parameters  
    38    INTEGER  ::   niter1max, niter2max   !: Maximum number of iterations for sinking 
    3938   REAL(wp) ::   rno3              !: ??? 
    4039   REAL(wp) ::   o2ut              !: ??? 

  • NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/trcini.F90

    r10372 r10377  
    205205      USE trcldf , ONLY:  trc_ldf_ini 
    206206      USE trcrad , ONLY:  trc_rad_ini 
     207      USE trcsink, ONLY:  trc_sink_ini 
    207208      ! 
    208209      INTEGER :: ierr 
    … …  
    214215                       !                          ! vertical diffusion: always implicit time stepping scheme 
    215216                       CALL  trc_rad_ini          ! positivity of passive tracers  
     217                       CALL  trc_sink_ini         ! Vertical sedimentation of particles 
    216218      ! 
    217219   END SUBROUTINE trc_ini_trp 
Note: See TracChangeset for help on using the changeset viewer.