Changeset 3101
- Timestamp:
- 2011-11-14T18:39:45+01:00 (13 years ago)
- Location:
- branches/2011/dev_NOC_UKMO_MERGE
- Files:
-
- 33 edited
- 2 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_NOC_UKMO_MERGE/DOC/TexFiles/Biblio/Biblio.bib
r3094 r3101 1275 1275 url = {http://dx.doi.org/10.1016/j.ocemod.2009.12.003}, 1276 1276 issn = {1463-5003}, 1277 } 1278 1279 @ARTICLE{HollowayOM86, 1280 author = {Greg Holloway}, 1281 title = {A Shelf Wave/Topographic Pump Drives Mean Coastal Circulation (part I)}, 1282 journal = OM, 1283 year = {1986}, 1284 volume = {68}, 1285 } 1286 1287 @ARTICLE{HollowayJPO92, 1288 author = {Greg Holloway}, 1289 title = {Representing Topographic Stress for Large-Scale Ocean Models}, 1290 journal = JPO, 1291 year = {1992}, 1292 volume = {22}, 1293 pages = {1033--1046}, 1294 } 1295 1296 @ARTICLE{HollowayJPO94, 1297 author = {Michael Eby and Greg Holloway}, 1298 title = {Sensitivity of a Large-Scale Ocean Model to a Parameterization of Topographic Stress}, 1299 journal = JPO, 1300 year = {1994}, 1301 volume = {24}, 1302 pages = {2577--2587}, 1303 } 1304 1305 @ARTICLE{HollowayJGR09, 1306 author = {Greg Holloway and Zeliang Wang}, 1307 title = {Representing eddy stress in an Arctic Ocean model}, 1308 journal = JGR, 1309 year = {2009}, 1310 doi = {10.1029/2008JC005169}, 1311 } 1312 1313 @ARTICLE{HollowayOM08, 1314 author = {Mathew Maltrud and Greg Holloway}, 1315 title = {Implementing biharmonic neptune in a global eddying ocean model}, 1316 journal = OM, 1317 year = {2008}, 1318 volume = {21}, 1319 pages = {22--34}, 1277 1320 } 1278 1321 -
branches/2011/dev_NOC_UKMO_MERGE/DOC/TexFiles/Chapters/Chap_DYN.tex
r2541 r3101 633 633 $\bullet$ Rotated axes scheme (rot) \citep{Thiem_Berntsen_OM06} (\np{ln\_dynhpg\_rot}=true) 634 634 635 Note that expression \eqref{Eq_dynhpg_sco} is used when the variable volume 635 $\bullet$ Pressure Jacobian scheme (prj) \citep{Thiem_Berntsen_OM06} (\np{ln\_dynhpg\_prj}=true) 636 637 Note that expression \eqref{Eq_dynhpg_sco} is commonly used when the variable volume 636 638 formulation is activated (\key{vvl}) because in that case, even with a flat bottom, 637 639 the coordinate surfaces are not horizontal but follow the free surface 638 \citep{Levier2007}. The other pressure gradient options are not yet available. 640 \citep{Levier2007}. Only the pressure jacobian scheme (\np{ln\_dynhpg\_prj}=true) is available as an 641 alternative to the default \np{ln\_dynhpg\_sco}=true when \key{vvl} is active. The pressure Jacobian scheme uses 642 a constrained cubic spline to reconstruct the density profile across the water column. This method 643 maintains the monotonicity between the density nodes and is of a higher order than the linear 644 interpolation method. The pressure can be calculated by analytical integration of the density profile and 645 a pressure Jacobian method is used to solve the horizontal pressure gradient. This method should 646 provide a more accurate calculation of the horizontal pressure gradient than the standard scheme. 639 647 640 648 %-------------------------------------------------------------------------------------------------------------- … … 1162 1170 1163 1171 % ================================================================ 1172 % Neptune effect 1173 % ================================================================ 1174 \section [Neptune effect (\textit{dynnept})] 1175 {Neptune effect (\mdl{dynnept})} 1176 \label{DYN_nept} 1177 1178 The "Neptune effect" (thus named in \citep{HollowayOM86}) is a 1179 parameterisation of the potentially large effect of topographic form stress 1180 (caused by eddies) in driving the ocean circulation. Originally developed for 1181 low-resolution models, in which it was applied via a Laplacian (second-order) 1182 diffusion-like term in the momentum equation, it can also be applied in eddy 1183 permitting or resolving models, in which a more scale-selective bilaplacian 1184 (fourth-order) implementation is preferred. This mechanism has a 1185 significant effect on boundary currents (including undercurrents), and the 1186 upwelling of deep water near continental shelves. 1187 1188 The theoretical basis for the method can be found in 1189 \citep{HollowayJPO92}, including the explanation of why form stress is not 1190 necessarily a drag force, but may actually drive the flow. 1191 \citep{HollowayJPO94} demonstrate the effects of the parameterisation in 1192 the GFDL-MOM model, at a horizontal resolution of about 1.8 degrees. 1193 \citep{HollowayOM08} demonstrate the biharmonic version of the 1194 parameterisation in a global run of the POP model, with an average horizontal 1195 grid spacing of about 32km. 1196 1197 The NEMO implementation is a simplified form of that supplied by 1198 Greg Holloway, the testing of which was described in \citep{HollowayJGR09}. 1199 The major simplification is that a time invariant Neptune velocity 1200 field is assumed. This is computed only once, during start-up, and 1201 made available to the rest of the code via a module. Vertical 1202 diffusive terms are also ignored, and the model topography itself 1203 is used, rather than a separate topographic dataset as in 1204 \citep{HollowayOM08}. This implementation is only in the iso-level 1205 formulation, as is the case anyway for the bilaplacian operator. 1206 1207 The velocity field is derived from a transport stream function given by: 1208 1209 \begin{equation} \label{Eq_dynnept_sf} 1210 \psi = -fL^2H 1211 \end{equation} 1212 1213 where $L$ is a latitude-dependant length scale given by: 1214 1215 \begin{equation} \label{Eq_dynnept_ls} 1216 L = l_1 + (l_2 -l_1)\left ( {1 + \cos 2\phi \over 2 } \right ) 1217 \end{equation} 1218 1219 where $\phi$ is latitude and $l_1$ and $l_2$ are polar and equatorial length scales respectively. 1220 Neptune velocity components, $u^*$, $v^*$ are derived from the stremfunction as: 1221 1222 \begin{equation} \label{Eq_dynnept_vel} 1223 u^* = -{1\over H} {\partial \psi \over \partial y}\ \ \ ,\ \ \ v^* = {1\over H} {\partial \psi \over \partial x} 1224 \end{equation} 1225 1226 \smallskip 1227 %----------------------------------------------namdom---------------------------------------------------- 1228 \namdisplay{namdyn_nept} 1229 %-------------------------------------------------------------------------------------------------------- 1230 \smallskip 1231 1232 The Neptune effect is enabled when \np{ln\_neptsimp}=true (default=false). 1233 \np{ln\_smooth\_neptvel} controls whether a scale-selective smoothing is applied 1234 to the Neptune effect flow field (default=false) (this smoothing method is as 1235 used by Holloway). \np{rn\_tslse} and \np{rn\_tslsp} are the equatorial and 1236 polar values respectively of the length-scale parameter $L$ used in determining 1237 the Neptune stream function \eqref{Eq_dynnept_sf} and \eqref{Eq_dynnept_ls}. 1238 Values at intermediate latitudes are given by a cosine fit, mimicking the 1239 variation of the deformation radius with latitude. The default values of 12km 1240 and 3km are those given in \citep{HollowayJPO94}, appropriate for a coarse 1241 resolution model. The finer resolution study of \citep{HollowayOM08} increased 1242 the values of L by a factor of $\sqrt 2$ to 17km and 4.2km, thus doubling the 1243 stream function for a given topography. 1244 1245 The simple formulation for ($u^*$, $v^*$) can give unacceptably large velocities 1246 in shallow water, and \citep{HollowayOM08} add an offset to the depth in the 1247 denominator to control this problem. In this implementation we offer instead (at 1248 the suggestion of G. Madec) the option of ramping down the Neptune flow field to 1249 zero over a finite depth range. The switch \np{ln\_neptramp} activates this 1250 option (default=false), in which case velocities at depths greater than 1251 \np{rn\_htrmax} are unaltered, but ramp down linearly with depth to zero at a 1252 depth of \np{rn\_htrmin} (and shallower). 1253 1254 % ================================================================ -
branches/2011/dev_NOC_UKMO_MERGE/DOC/TexFiles/Chapters/Chap_MISC.tex
r2541 r3101 253 253 Note this implementation may be sensitive to the optimization level. 254 254 255 \subsection{MPP scalability} 256 \label{MISC_mppsca} 257 258 The default method of communicating values across the north-fold in distributed memory applications 259 (\key{mpp\_mpi}) uses a \textsc{MPI\_ALLGATHER} function to exchange values from each processing 260 region in the northern row with every other processing region in the northern row. This enables a 261 global width array containing the top 4 rows to be collated on every northern row processor and then 262 folded with a simple algorithm. Although conceptually simple, this "All to All" communication will 263 hamper performance scalability for large numbers of northern row processors. From version 3.4 264 onwards an alternative method is available which only performs direct "Peer to Peer" communications 265 between each processor and its immediate "neighbours" across the fold line. This is achieved by 266 using the default \textsc{MPI\_ALLGATHER} method during initialisation to help identify the "active" 267 neighbours. Stored lists of these neighbours are then used in all subsequent north-fold exchanges to 268 restrict exchanges to those between associated regions. The collated global width array for each 269 region is thus only partially filled but is guaranteed to be set at all the locations actually 270 required by each individual for the fold operation. This alternative method should give identical 271 results to the default \textsc{ALLGATHER} method and is recommended for large values of \np{jpni}. 272 The new method is activated by setting \np{ln\_nnogather} to be true ({\bf nammpp}). The 273 reproducibility of results using the two methods should be confirmed for each new, non-reference 274 configuration. 255 275 256 276 % ================================================================ -
branches/2011/dev_NOC_UKMO_MERGE/DOC/TexFiles/Chapters/Chap_ZDF.tex
r2541 r3101 1 1 % ================================================================ 2 % Chapter ÑVertical Ocean Physics (ZDF)2 % Chapter Vertical Ocean Physics (ZDF) 3 3 % ================================================================ 4 4 \chapter{Vertical Ocean Physics (ZDF)} … … 539 539 the clipping factor is of crucial importance for the entrainment depth predicted in 540 540 stably stratified situations, and that its value has to be chosen in accordance 541 with the algebraic model for the turbulent ßuxes. The clipping is only activated541 with the algebraic model for the turbulent fluxes. The clipping is only activated 542 542 if \np{ln\_length\_lim}=true, and the $c_{lim}$ is set to the \np{rn\_clim\_galp} value. 543 543 … … 981 981 reduced as necessary to ensure stability; these changes are not reported. 982 982 983 Limits on the bottom friction coefficient are not imposed if the user has elected to 984 handle the bottom friction implicitly (see \S\ref{ZDF_bfr_imp}). The number of potential 985 breaches of the explicit stability criterion are still reported for information purposes. 986 987 % ------------------------------------------------------------------------------------------------------------- 988 % Implicit Bottom Friction 989 % ------------------------------------------------------------------------------------------------------------- 990 \subsection{Implicit Bottom Friction (\np{ln\_bfrimp}$=$\textit{T})} 991 \label{ZDF_bfr_imp} 992 993 An optional implicit form of bottom friction has been implemented to improve 994 model stability. We recommend this option for shelf sea and coastal ocean applications, especially 995 for split-explicit time splitting. This option can be invoked by setting \np{ln\_bfrimp} 996 to \textit{true} in the \textit{nambfr} namelist. This option requires \np{ln\_zdfexp} to be \textit{false} 997 in the \textit{namzdf} namelist. 998 999 This implementation is realised in \mdl{dynzdf\_imp} and \mdl{dynspg\_ts}. In \mdl{dynzdf\_imp}, the 1000 bottom boundary condition is implemented implicitly. 1001 1002 \begin{equation} \label{Eq_dynzdf_bfr} 1003 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{mbk} 1004 = \binom{c_{b}^{u}u^{n+1}_{mbk}}{c_{b}^{v}v^{n+1}_{mbk}} 1005 \end{equation} 1006 1007 where $mbk$ is the layer number of the bottom wet layer. superscript $n+1$ means the velocity used in the 1008 friction formula is to be calculated, so, it is implicit. 1009 1010 If split-explicit time splitting is used, care must be taken to avoid the double counting of 1011 the bottom friction in the 2-D barotropic momentum equations. As NEMO only updates the barotropic 1012 pressure gradient and Coriolis' forcing terms in the 2-D barotropic calculation, we need to remove 1013 the bottom friction induced by these two terms which has been included in the 3-D momentum trend 1014 and update it with the latest value. On the other hand, the bottom friction contributed by the 1015 other terms (e.g. the advection term, viscosity term) has been included in the 3-D momentum equations 1016 and should not be added in the 2-D barotropic mode. 1017 1018 The implementation of the implicit bottom friction in \mdl{dynspg\_ts} is done in two steps as the 1019 following: 1020 1021 \begin{equation} \label{Eq_dynspg_ts_bfr1} 1022 \frac{\textbf{U}_{med}-\textbf{U}^{m-1}}{2\Delta t}=-g\nabla\eta-f\textbf{k}\times\textbf{U}^{m}+c_{b} 1023 \left(\textbf{U}_{med}-\textbf{U}^{m-1}\right) 1024 \end{equation} 1025 \begin{equation} \label{Eq_dynspg_ts_bfr2} 1026 \frac{\textbf{U}^{m+1}-\textbf{U}_{med}}{2\Delta t}=\textbf{T}+ 1027 \left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{U}^{'}\right)- 1028 2\Delta t_{bc}c_{b}\left(g\nabla\eta^{'}+f\textbf{k}\times\textbf{u}_{b}\right) 1029 \end{equation} 1030 1031 where $\textbf{T}$ is the vertical integrated 3-D momentum trend. We assume the leap-frog time-stepping 1032 is used here. $\Delta t$ is the barotropic mode time step and $\Delta t_{bc}$ is the baroclinic mode time step. 1033 $c_{b}$ is the friction coefficient. $\eta$ is the sea surface level calculated in the barotropic loops 1034 while $\eta^{'}$ is the sea surface level used in the 3-D baroclinic mode. $\textbf{u}_{b}$ is the bottom 1035 layer horizontal velocity. 1036 1037 1038 1039 983 1040 % ------------------------------------------------------------------------------------------------------------- 984 1041 % Bottom Friction with split-explicit time splitting 985 1042 % ------------------------------------------------------------------------------------------------------------- 986 \subsection{Bottom Friction with split-explicit time splitting }1043 \subsection{Bottom Friction with split-explicit time splitting (\np{ln\_bfrimp}$=$\textit{F})} 987 1044 \label{ZDF_bfr_ts} 988 1045 … … 993 1050 {\key{dynspg\_flt}). Extra attention is required, however, when using 994 1051 split-explicit time stepping (\key{dynspg\_ts}). In this case the free surface 995 equation is solved with a small time step \np{ nn\_baro}*\np{rn\_rdt}, while the three996 dimensional prognostic variables are solved with a longer time step that is a997 multiple of \np{rn\_rdt}. The trend in the barotropic momentum due to bottom1052 equation is solved with a small time step \np{rn\_rdt}/\np{nn\_baro}, while the three 1053 dimensional prognostic variables are solved with the longer time step 1054 of \np{rn\_rdt} seconds. The trend in the barotropic momentum due to bottom 998 1055 friction appropriate to this method is that given by the selected parameterisation 999 1056 ($i.e.$ linear or non-linear bottom friction) computed with the evolving velocities … … 1018 1075 \end{enumerate} 1019 1076 1020 Note that the use of an implicit formulation 1077 Note that the use of an implicit formulation within the barotropic loop 1021 1078 for the bottom friction trend means that any limiting of the bottom friction coefficient 1022 1079 in \mdl{dynbfr} does not adversely affect the solution when using split-explicit time 1023 1080 splitting. This is because the major contribution to bottom friction is likely to come from 1024 the barotropic component which uses the unrestricted value of the coefficient. 1025 1026 The implicit formulation takes the form: 1081 the barotropic component which uses the unrestricted value of the coefficient. However, if the 1082 limiting is thought to be having a major effect (a more likely prospect in coastal and shelf seas 1083 applications) then the fully implicit form of the bottom friction should be used (see \S\ref{ZDF_bfr_imp} ) 1084 which can be selected by setting \np{ln\_bfrimp} $=$ \textit{true}. 1085 1086 Otherwise, the implicit formulation takes the form: 1027 1087 \begin{equation} \label{Eq_zdfbfr_implicitts} 1028 1088 \bar{U}^{t+ \rdt} = \; \left [ \bar{U}^{t-\rdt}\; + 2 \rdt\;RHS \right ] / \left [ 1 - 2 \rdt \;c_b^{u} / H_e \right ] … … 1091 1151 The essential goal of the parameterization is to represent the momentum 1092 1152 exchange between the barotropic tides and the unrepresented internal waves 1093 induced by the tidal ßow over rough topography in a stratified ocean.1153 induced by the tidal flow over rough topography in a stratified ocean. 1094 1154 In the current version of \NEMO, the map is built from the output of 1095 1155 the barotropic global ocean tide model MOG2D-G \citep{Carrere_Lyard_GRL03}. -
branches/2011/dev_NOC_UKMO_MERGE/DOC/TexFiles/Namelist/nambfr
r2540 r3101 9 9 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 10 10 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d=T) 11 ln_bfrimp = .false. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 11 12 / -
branches/2011/dev_NOC_UKMO_MERGE/DOC/TexFiles/Namelist/namdyn_hpg
r2540 r3101 9 9 ln_hpg_djc = .false. ! s-coordinate (Density Jacobian with Cubic polynomial) 10 10 ln_hpg_rot = .false. ! s-coordinate (ROTated axes scheme) 11 ln_hpg_prj = .false. ! s-coordinate (Pressure Jacobian scheme) 11 12 rn_gamma = 0.e0 ! weighting coefficient (wdj scheme) 12 13 ln_dynhpg_imp = .false. ! time stepping: semi-implicit time scheme (T) -
branches/2011/dev_NOC_UKMO_MERGE/DOC/TexFiles/Namelist/namtra_ldf
r2540 r3101 9 9 ln_traldf_hor = .false. ! horizontal (geopotential) (require "key_ldfslp" when ln_sco=T) 10 10 ln_traldf_iso = .true. ! iso-neutral (require "key_ldfslp") 11 ln_traldf_grif = .false. ! griffies skew flux formulation (require "key_ldfslp") ! UNDER TEST, DO NOT USE 12 ln_traldf_gdia = .false. ! griffies operator strfn diagnostics (require "key_ldfslp") ! UNDER TEST, DO NOT USE 11 ln_traldf_grif = .false. ! griffies skew flux formulation (require "key_ldfslp") 12 ln_traldf_gdia = .false. ! griffies operator strfn diagnostics (require "key_ldfslp") 13 ln_triad_iso = .false. ! griffies operator calculates triads twice => pure lateral mixing in ML (require "key_ldfslp") 14 ln_botmix_grif = .false. ! griffies operator with lateral mixing on bottom (require "key_ldfslp") 13 15 ! ! Coefficient 14 16 rn_aht_0 = 2000. ! horizontal eddy diffusivity for tracers [m2/s] -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/ARCH/arch-ALTIX_NAUTILUS4.fcm
r2364 r3101 22 22 # Note use of -Bstatic because the library root directories are not accessible to the back-end compute nodes 23 23 %NCDF_LIB -L%HDF5_HOME/lib -L%NCDF_HOME/lib -Bstatic -lnetcdf -lhdf5_fortran -lhdf5_hl -lhdf5 -Bdynamic -lz 24 %FC mpif9024 %FC ifort 25 25 %FCFLAGS -r8 -O3 -xT -ip -vec-report0 26 26 %FFLAGS -r8 -O3 -xT -ip -vec-report0 27 %LD mpif9027 %LD ifort 28 28 %FPPFLAGS -P -C -traditional 29 %LDFLAGS 29 %LDFLAGS -lmpi 30 30 %AR ar 31 31 %ARFLAGS -r -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/CONFIG/GYRE/EXP00/namelist
r3094 r3101 408 408 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 409 409 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d = .true.) 410 ln_bfrimp = .false. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 410 411 / 411 412 !----------------------------------------------------------------------- … … 454 455 ln_traadv_muscl2 = .false. ! MUSCL2 scheme + cen2 at boundaries 455 456 ln_traadv_ubs = .false. ! UBS scheme 456 ln_traadv_qck = .false. ! QU CIKEST scheme457 ln_traadv_qck = .false. ! QUICKEST scheme 457 458 / 458 459 !----------------------------------------------------------------------- … … 466 467 ln_traldf_hor = .false. ! horizontal (geopotential) (require "key_ldfslp" when ln_sco=T) 467 468 ln_traldf_iso = .true. ! iso-neutral (require "key_ldfslp") 468 ln_traldf_grif = .false. ! griffies skew flux formulation (require "key_ldfslp") ! UNDER TEST, DO NOT USE 469 ln_traldf_gdia = .false. ! griffies operator strfn diagnostics (require "key_ldfslp") ! UNDER TEST, DO NOT USE 469 ln_traldf_grif = .false. ! griffies skew flux formulation (require "key_ldfslp") 470 ln_traldf_gdia = .false. ! griffies operator strfn diagnostics (require "key_ldfslp") 471 ln_triad_iso = .false. ! griffies operator calculates triads twice => pure lateral mixing in ML (require "key_ldfslp") 472 ln_botmix_grif = .false. ! griffies operator with lateral mixing on bottom (require "key_ldfslp") 470 473 ! ! Coefficient 471 474 rn_aht_0 = 1000. ! horizontal eddy diffusivity for tracers [m2/s] … … 523 526 ln_hpg_djc = .false. ! s-coordinate (Density Jacobian with Cubic polynomial) 524 527 ln_hpg_rot = .false. ! s-coordinate (ROTated axes scheme) 528 ln_hpg_prj = .false. ! s-coordinate (Pressure Jacobian scheme) 525 529 rn_gamma = 0.e0 ! weighting coefficient (wdj scheme) 526 530 ln_dynhpg_imp = .false. ! time stepping: semi-implicit time scheme (T) … … 686 690 ! buffer blocking send or immediate non-blocking sends, resp. 687 691 nn_buffer = 0 ! size in bytes of exported buffer ('B' case), 0 no exportation 692 ln_nnogather= .false. ! activate code to avoid mpi_allgather use at the northfold 688 693 jpni = 0 ! jpni number of processors following i (set automatically if < 1) 689 694 jpnj = 0 ! jpnj number of processors following j (set automatically if < 1) … … 854 859 salfixmin = -9999 ! Minimum salinity after applying the increments 855 860 / 861 !----------------------------------------------------------------------- 862 &namdyn_nept ! Neptune effect (simplified: lateral and vertical diffusions removed) 863 !----------------------------------------------------------------------- 864 ! Suggested lengthscale values are those of Eby & Holloway (1994) for a coarse model 865 ln_neptsimp = .false. ! yes/no use simplified neptune 866 867 ln_smooth_neptvel = .false. ! yes/no smooth zunep, zvnep 868 rn_tslse = 1.2e4 ! value of lengthscale L at the equator 869 rn_tslsp = 3.0e3 ! value of lengthscale L at the pole 870 ! Specify whether to ramp down the Neptune velocity in shallow 871 ! water, and if so the depth range controlling such ramping down 872 ln_neptramp = .false. ! ramp down Neptune velocity in shallow water 873 rn_htrmin = 100.0 ! min. depth of transition range 874 rn_htrmax = 200.0 ! max. depth of transition range 875 / -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist
r3094 r3101 408 408 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 409 409 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d=T) 410 ln_bfrimp = .false. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 410 411 / 411 412 !----------------------------------------------------------------------- … … 523 524 ln_hpg_djc = .false. ! s-coordinate (Density Jacobian with Cubic polynomial) 524 525 ln_hpg_rot = .false. ! s-coordinate (ROTated axes scheme) 526 ln_hpg_prj = .false. ! s-coordinate (Pressure Jacobian scheme) 525 527 rn_gamma = 0.e0 ! weighting coefficient (wdj scheme) 526 528 ln_dynhpg_imp = .false. ! time stepping: semi-implicit time scheme (T) -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
r3094 r3101 408 408 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 409 409 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d=T) 410 ln_bfrimp = .false. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 410 411 / 411 412 !----------------------------------------------------------------------- … … 454 455 ln_traadv_muscl2 = .false. ! MUSCL2 scheme + cen2 at boundaries 455 456 ln_traadv_ubs = .false. ! UBS scheme 456 ln_traadv_qck = .false. ! QU CIKEST scheme457 ln_traadv_qck = .false. ! QUICKEST scheme 457 458 / 458 459 !----------------------------------------------------------------------- … … 466 467 ln_traldf_hor = .false. ! horizontal (geopotential) (require "key_ldfslp" when ln_sco=T) 467 468 ln_traldf_iso = .true. ! iso-neutral (require "key_ldfslp") 468 ln_traldf_grif = .false. ! griffies skew flux formulation (require "key_ldfslp") ! UNDER TEST, DO NOT USE 469 ln_traldf_gdia = .false. ! griffies operator strfn diagnostics (require "key_ldfslp") ! UNDER TEST, DO NOT USE 470 ! ! Coefficient 469 ln_traldf_grif = .false. ! griffies skew flux formulation (require "key_ldfslp") 470 ln_traldf_gdia = .false. ! griffies operator strfn diagnostics (require "key_ldfslp") 471 ln_triad_iso = .false. ! griffies operator calculates triads twice => pure lateral mixing in ML (require "key_ldfslp") 472 ln_botmix_grif = .false. ! griffies operator with lateral mixing on bottom (require "key_ldfslp") 473 ! Coefficient 471 474 rn_aht_0 = 2000. ! horizontal eddy diffusivity for tracers [m2/s] 472 475 rn_ahtb_0 = 0. ! background eddy diffusivity for ldf_iso [m2/s] … … 523 526 ln_hpg_djc = .false. ! s-coordinate (Density Jacobian with Cubic polynomial) 524 527 ln_hpg_rot = .false. ! s-coordinate (ROTated axes scheme) 528 ln_hpg_prj = .false. ! s-coordinate (Pressure Jacobian scheme) 525 529 rn_gamma = 0.e0 ! weighting coefficient (wdj scheme) 526 530 ln_dynhpg_imp = .false. ! time stepping: semi-implicit time scheme (T) … … 686 690 ! buffer blocking send or immediate non-blocking sends, resp. 687 691 nn_buffer = 0 ! size in bytes of exported buffer ('B' case), 0 no exportation 692 ln_nnogather= .false. ! activate code to avoid mpi_allgather use at the northfold 688 693 jpni = 0 ! jpni number of processors following i (set automatically if < 1) 689 694 jpnj = 0 ! jpnj number of processors following j (set automatically if < 1) … … 854 859 salfixmin = -9999 ! Minimum salinity after applying the increments 855 860 / 861 !----------------------------------------------------------------------- 862 &namdyn_nept ! Neptune effect (simplified: lateral and vertical diffusions removed) 863 !----------------------------------------------------------------------- 864 ! Suggested lengthscale values are those of Eby & Holloway (1994) for a coarse model 865 ln_neptsimp = .false. ! yes/no use simplified neptune 866 867 ln_smooth_neptvel = .false. ! yes/no smooth zunep, zvnep 868 rn_tslse = 1.2e4 ! value of lengthscale L at the equator 869 rn_tslsp = 3.0e3 ! value of lengthscale L at the pole 870 ! Specify whether to ramp down the Neptune velocity in shallow 871 ! water, and if so the depth range controlling such ramping down 872 ln_neptramp = .true. ! ramp down Neptune velocity in shallow water 873 rn_htrmin = 100.0 ! min. depth of transition range 874 rn_htrmax = 200.0 ! max. depth of transition range 875 / -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist
r3094 r3101 408 408 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 409 409 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d = .true.) 410 ln_bfrimp = .false. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 410 411 / 411 412 !----------------------------------------------------------------------- … … 454 455 ln_traadv_muscl2 = .false. ! MUSCL2 scheme + cen2 at boundaries 455 456 ln_traadv_ubs = .false. ! UBS scheme 456 ln_traadv_qck = .false. ! QU CIKEST scheme457 ln_traadv_qck = .false. ! QUICKEST scheme 457 458 / 458 459 !----------------------------------------------------------------------- … … 466 467 ln_traldf_hor = .false. ! horizontal (geopotential) (require "key_ldfslp" when ln_sco=T) 467 468 ln_traldf_iso = .true. ! iso-neutral (require "key_ldfslp") 468 ln_traldf_grif = .false. ! griffies skew flux formulation (require "key_ldfslp") ! UNDER TEST, DO NOT USE 469 ln_traldf_gdia = .false. ! griffies operator strfn diagnostics (require "key_ldfslp") ! UNDER TEST, DO NOT USE 470 ! ! Coefficient 469 ln_traldf_grif = .false. ! griffies skew flux formulation (require "key_ldfslp") 470 ln_traldf_gdia = .false. ! griffies operator strfn diagnostics (require "key_ldfslp") 471 ln_triad_iso = .false. ! griffies operator calculates triads twice => pure lateral mixing in ML (require "key_ldfslp") 472 ln_botmix_grif = .false. ! griffies operator with lateral mixing on bottom (require "key_ldfslp") 473 ! Coefficient 471 474 rn_aht_0 = 2000. ! horizontal eddy diffusivity for tracers [m2/s] 472 475 rn_ahtb_0 = 0. ! background eddy diffusivity for ldf_iso [m2/s] … … 524 527 ln_hpg_djc = .false. ! s-coordinate (Density Jacobian with Cubic polynomial) 525 528 ln_hpg_rot = .false. ! s-coordinate (ROTated axes scheme) 529 ln_hpg_prj = .false. ! s-coordinate (Pressure Jacobian scheme) 526 530 rn_gamma = 0.e0 ! weighting coefficient (wdj scheme) 527 531 ln_dynhpg_imp = .false. ! time stepping: semi-implicit time scheme (T) … … 700 704 ! buffer blocking send or immediate non-blocking sends, resp. 701 705 nn_buffer = 0 ! size in bytes of exported buffer ('B' case), 0 no exportation 706 ln_nnogather= .false. ! activate code to avoid mpi_allgather use at the northfold 702 707 jpni = 0 ! jpni number of processors following i (set automatically if < 1) 703 708 jpnj = 0 ! jpnj number of processors following j (set automatically if < 1) … … 869 874 salfixmin = -9999 ! Minimum salinity after applying the increments 870 875 / 876 !----------------------------------------------------------------------- 877 &namdyn_nept ! Neptune effect (simplified: lateral and vertical diffusions removed) 878 !----------------------------------------------------------------------- 879 ! Suggested lengthscale values are those of Eby & Holloway (1994) for a coarse model 880 ln_neptsimp = .false. ! yes/no use simplified neptune 881 882 ln_smooth_neptvel = .false. ! yes/no smooth zunep, zvnep 883 rn_tslse = 1.2e4 ! value of lengthscale L at the equator 884 rn_tslsp = 3.0e3 ! value of lengthscale L at the pole 885 ! Specify whether to ramp down the Neptune velocity in shallow 886 ! water, and if so the depth range controlling such ramping down 887 ln_neptramp = .false. ! ramp down Neptune velocity in shallow water 888 rn_htrmin = 100.0 ! min. depth of transition range 889 rn_htrmax = 200.0 ! max. depth of transition range 890 / -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/CONFIG/POMME/EXP00/namelist
r3094 r3101 408 408 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 409 409 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d=T) 410 ln_bfrimp = .false. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 410 411 / 411 412 !----------------------------------------------------------------------- … … 454 455 ln_traadv_muscl2 = .false. ! MUSCL2 scheme + cen2 at boundaries 455 456 ln_traadv_ubs = .false. ! UBS scheme 456 ln_traadv_qck = .false. ! QU CIKEST scheme457 ln_traadv_qck = .false. ! QUICKEST scheme 457 458 / 458 459 !----------------------------------------------------------------------- … … 466 467 ln_traldf_hor = .false. ! horizontal (geopotential) (require "key_ldfslp" when ln_sco=T) 467 468 ln_traldf_iso = .true. ! iso-neutral (require "key_ldfslp") 468 ln_traldf_grif = .false. ! griffies skew flux formulation (require "key_ldfslp") ! UNDER TEST, DO NOT USE 469 ln_traldf_gdia = .false. ! griffies operator strfn diagnostics (require "key_ldfslp") ! UNDER TEST, DO NOT USE 469 ln_traldf_grif = .false. ! griffies skew flux formulation (require "key_ldfslp") 470 ln_traldf_gdia = .false. ! griffies operator strfn diagnostics (require "key_ldfslp") 471 ln_triad_iso = .false. ! griffies operator calculates triads twice => pure lateral mixing in ML (require "key_ldfslp") 472 ln_botmix_grif = .false. ! griffies operator with lateral mixing on bottom (require "key_ldfslp") 470 473 ! ! Coefficient 471 474 rn_aht_0 = 300. ! horizontal eddy diffusivity for tracers [m2/s] … … 523 526 ln_hpg_djc = .false. ! s-coordinate (Density Jacobian with Cubic polynomial) 524 527 ln_hpg_rot = .false. ! s-coordinate (ROTated axes scheme) 528 ln_hpg_prj = .false. ! s-coordinate (Pressure Jacobian scheme) 525 529 rn_gamma = 0.e0 ! weighting coefficient (wdj scheme) 526 530 ln_dynhpg_imp = .true. ! time stepping: semi-implicit time scheme (T) … … 859 863 salfixmin = -9999 ! Minimum salinity after applying the increments 860 864 / 865 !----------------------------------------------------------------------- 866 &namdyn_nept ! Neptune effect (simplified: lateral and vertical diffusions removed) 867 !----------------------------------------------------------------------- 868 ! Suggested lengthscale values are those of Eby & Holloway (1994) for a coarse model 869 ln_neptsimp = .false. ! yes/no use simplified neptune 870 871 ln_smooth_neptvel = .false. ! yes/no smooth zunep, zvnep 872 rn_tslse = 1.2e4 ! value of lengthscale L at the equator 873 rn_tslsp = 3.0e3 ! value of lengthscale L at the pole 874 ! Specify whether to ramp down the Neptune velocity in shallow 875 ! water, and if so the depth range controlling such ramping down 876 ln_neptramp = .false. ! ramp down Neptune velocity in shallow water 877 rn_htrmin = 100.0 ! min. depth of transition range 878 rn_htrmax = 200.0 ! max. depth of transition range 879 / -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r2779 r3101 24 24 PRIVATE 25 25 26 PUBLIC dom_vvl ! called by domain.F9027 PUBLIC dom_vvl_ alloc ! called by nemogcm.F9028 29 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ee_t, ee_u, ee_v, ee_f !: ??? 30 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mut , muu , muv , muf !: ???26 PUBLIC dom_vvl ! called by domain.F90 27 PUBLIC dom_vvl_2 ! called by domain.F90 28 PUBLIC dom_vvl_alloc ! called by nemogcm.F90 29 30 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mut , muu , muv , muf !: 1/H_0 at t-,u-,v-,f-points 31 31 32 32 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: r2dt ! vertical profile time-step, = 2 rdttra … … 49 49 ! 50 50 ALLOCATE( mut (jpi,jpj,jpk) , muu (jpi,jpj,jpk) , muv (jpi,jpj,jpk) , muf (jpi,jpj,jpk) , & 51 & ee_t(jpi,jpj) , ee_u(jpi,jpj) , ee_v(jpi,jpj) , ee_f(jpi,jpj) , &52 51 & r2dt (jpk) , STAT=dom_vvl_alloc ) 53 52 ! … … 62 61 !! *** ROUTINE dom_vvl *** 63 62 !! 64 !! ** Purpose : compute coefficients muX at T-U-V-F points to spread 65 !! ssh over the whole water column (scale factors) 63 !! ** Purpose : compute mu coefficients at t-, u-, v- and f-points to 64 !! spread ssh over the whole water column (scale factors) 65 !! set the before and now ssh at u- and v-points 66 !! (also f-point in now case) 66 67 !!---------------------------------------------------------------------- 67 68 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 68 USE wrk_nemo, ONLY: z s_t => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_3! 2D workspace69 USE wrk_nemo, ONLY: zee_t => wrk_2d_1, zee_u => wrk_2d_2, zee_v => wrk_2d_3, zee_f => wrk_2d_4 ! 2D workspace 69 70 ! 70 71 INTEGER :: ji, jj, jk ! dummy loop indices 71 REAL(wp) :: zcoefu , zcoefv , zcoeff! local scalars72 REAL(wp) :: zv _t_ij, zv_t_ip1j, zv_t_ijp1, zv_t_ip1jp1 ! - -73 !!---------------------------------------------------------------------- 74 75 IF( wrk_in_use(2, 1,2,3 ) ) THEN72 REAL(wp) :: zcoefu, zcoefv , zcoeff ! local scalars 73 REAL(wp) :: zvt , zvt_ip1, zvt_jp1, zvt_ip1jp1 ! - - 74 !!---------------------------------------------------------------------- 75 76 IF( wrk_in_use(2, 1,2,3,4) ) THEN 76 77 CALL ctl_stop('dom_vvl: requested workspace arrays unavailable') ; RETURN 77 78 ENDIF … … 97 98 98 99 ! !== mu computation ==! 99 ee_t(:,:) = fse3t_0(:,:,1) ! Lower bound : thickness of the first model level100 ee_u(:,:) = fse3u_0(:,:,1)101 ee_v(:,:) = fse3v_0(:,:,1)102 ee_f(:,:) = fse3f_0(:,:,1)100 zee_t(:,:) = fse3t_0(:,:,1) ! Lower bound : thickness of the first model level 101 zee_u(:,:) = fse3u_0(:,:,1) 102 zee_v(:,:) = fse3v_0(:,:,1) 103 zee_f(:,:) = fse3f_0(:,:,1) 103 104 DO jk = 2, jpkm1 ! Sum of the masked vertical scale factors 104 ee_t(:,:) =ee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk)105 ee_u(:,:) =ee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk)106 ee_v(:,:) =ee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk)105 zee_t(:,:) = zee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 106 zee_u(:,:) = zee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 107 zee_v(:,:) = zee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 107 108 DO jj = 1, jpjm1 ! f-point : fmask=shlat at coasts, use the product of umask 108 ee_f(:,jj) =ee_f(:,jj) + fse3f_0(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)109 zee_f(:,jj) = zee_f(:,jj) + fse3f_0(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 109 110 END DO 110 111 END DO 111 112 ! ! Compute and mask the inverse of the local depth at T, U, V and F points 112 ee_t(:,:) = 1. /ee_t(:,:) * tmask(:,:,1)113 ee_u(:,:) = 1. /ee_u(:,:) * umask(:,:,1)114 ee_v(:,:) = 1. /ee_v(:,:) * vmask(:,:,1)113 zee_t(:,:) = 1._wp / zee_t(:,:) * tmask(:,:,1) 114 zee_u(:,:) = 1._wp / zee_u(:,:) * umask(:,:,1) 115 zee_v(:,:) = 1._wp / zee_v(:,:) * vmask(:,:,1) 115 116 DO jj = 1, jpjm1 ! f-point case fmask cannot be used 116 ee_f(:,jj) = 1. /ee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1)117 END DO 118 CALL lbc_lnk( ee_f, 'F', 1. )! lateral boundary condition on ee_f117 zee_f(:,jj) = 1._wp / zee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 118 END DO 119 CALL lbc_lnk( zee_f, 'F', 1. ) ! lateral boundary condition on ee_f 119 120 ! 120 121 DO jk = 1, jpk ! mu coefficients 121 mut(:,:,jk) = ee_t(:,:) * tmask(:,:,jk) ! T-point at T levels122 muu(:,:,jk) = ee_u(:,:) * umask(:,:,jk) ! U-point at T levels123 muv(:,:,jk) = ee_v(:,:) * vmask(:,:,jk) ! V-point at T levels122 mut(:,:,jk) = zee_t(:,:) * tmask(:,:,jk) ! T-point at T levels 123 muu(:,:,jk) = zee_u(:,:) * umask(:,:,jk) ! U-point at T levels 124 muv(:,:,jk) = zee_v(:,:) * vmask(:,:,jk) ! V-point at T levels 124 125 END DO 125 126 DO jk = 1, jpk ! F-point : fmask=shlat at coasts, use the product of umask 126 127 DO jj = 1, jpjm1 127 muf(:,jj,jk) = ee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk) ! at T levels128 END DO 129 muf(:,jpj,jk) = 0. e0128 muf(:,jj,jk) = zee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk) ! at T levels 129 END DO 130 muf(:,jpj,jk) = 0._wp 130 131 END DO 131 132 CALL lbc_lnk( muf, 'F', 1. ) ! lateral boundary condition … … 139 140 END DO 140 141 141 ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations142 ! for ssh and scale factors143 zs_t (:,:) = e1t(:,:) * e2t(:,:)144 zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )145 zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )146 147 142 DO jj = 1, jpjm1 ! initialise before and now Sea Surface Height at u-, v-, f-points 148 143 DO ji = 1, jpim1 ! NO vector opt. 149 zcoefu = umask(ji,jj,1) * zs_u_1(ji,jj)150 zcoefv = vmask(ji,jj,1) * zs_v_1(ji,jj)151 zcoeff = 0. 5 * umask(ji,jj,1) * umask(ji,jj+1,1) / ( e1f(ji,jj) * e2f(ji,jj))152 ! before fields153 zv _t_ij = zs_t(ji ,jj ) * sshb(ji ,jj )154 zv _t_ip1j = zs_t(ji+1,jj ) * sshb(ji+1,jj )155 zv _t_ijp1 = zs_t(ji ,jj+1) * sshb(ji ,jj+1)156 sshu_b(ji,jj) = zcoefu * ( zv _t_ij + zv_t_ip1j)157 sshv_b(ji,jj) = zcoefv * ( zv _t_ij + zv_t_ijp1 )158 ! now fields159 zv _t_ij = zs_t(ji ,jj ) * sshn(ji ,jj )160 zv _t_ip1j = zs_t(ji+1,jj ) * sshn(ji+1,jj )161 zv _t_ijp1 = zs_t(ji ,jj+1) * sshn(ji ,jj+1)162 zv _t_ip1jp1 = zs_t(ji ,jj+1) * sshn(ji,jj+1)163 sshu_n(ji,jj) = zcoefu * ( zv _t_ij + zv_t_ip1j)164 sshv_n(ji,jj) = zcoefv * ( zv _t_ij + zv_t_ijp1 )165 sshf_n(ji,jj) = zcoeff * ( zv _t_ij + zv_t_ip1j + zv_t_ijp1 + zv_t_ip1jp1 )144 zcoefu = 0.50_wp / ( e1u(ji,jj) * e2u(ji,jj) ) * umask(ji,jj,1) 145 zcoefv = 0.50_wp / ( e1v(ji,jj) * e2v(ji,jj) ) * vmask(ji,jj,1) 146 zcoeff = 0.25_wp / ( e1f(ji,jj) * e2f(ji,jj) ) * umask(ji,jj,1) * umask(ji,jj+1,1) 147 ! 148 zvt = e1e2t(ji ,jj ) * sshb(ji ,jj ) ! before fields 149 zvt_ip1 = e1e2t(ji+1,jj ) * sshb(ji+1,jj ) 150 zvt_jp1 = e1e2t(ji ,jj+1) * sshb(ji ,jj+1) 151 sshu_b(ji,jj) = zcoefu * ( zvt + zvt_ip1 ) 152 sshv_b(ji,jj) = zcoefv * ( zvt + zvt_jp1 ) 153 ! 154 zvt = e1e2t(ji ,jj ) * sshn(ji ,jj ) ! now fields 155 zvt_ip1 = e1e2t(ji+1,jj ) * sshn(ji+1,jj ) 156 zvt_jp1 = e1e2t(ji ,jj+1) * sshn(ji ,jj+1) 157 zvt_ip1jp1 = e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1) 158 sshu_n(ji,jj) = zcoefu * ( zvt + zvt_ip1 ) 159 sshv_n(ji,jj) = zcoefv * ( zvt + zvt_jp1 ) 160 sshf_n(ji,jj) = zcoeff * ( zvt + zvt_ip1 + zvt_jp1 + zvt_ip1jp1 ) 166 161 END DO 167 162 END DO … … 169 164 CALL lbc_lnk( sshv_n, 'V', 1. ) ; CALL lbc_lnk( sshv_b, 'V', 1. ) 170 165 CALL lbc_lnk( sshf_n, 'F', 1. ) 171 172 ! initialise before scale factors at (u/v)-points 173 ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 174 DO jk = 1, jpkm1 175 DO jj = 1, jpjm1 176 DO ji = 1, jpim1 177 zv_t_ij = zs_t(ji ,jj ) * fse3t_b(ji ,jj ,jk) 178 zv_t_ip1j = zs_t(ji+1,jj ) * fse3t_b(ji+1,jj ,jk) 179 zv_t_ijp1 = zs_t(ji ,jj+1) * fse3t_b(ji ,jj+1,jk) 180 fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 181 fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 182 END DO 183 END DO 184 END DO 185 CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. ) ! lateral boundary conditions 186 CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 187 ! Add initial scale factor to scale factor anomaly 188 fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 189 fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 190 ! 191 IF( wrk_not_released(2, 1,2,3) ) CALL ctl_stop('dom_vvl: failed to release workspace arrays') 166 ! 167 IF( wrk_not_released(2, 1,2,3,4) ) CALL ctl_stop('dom_vvl: failed to release workspace arrays') 192 168 ! 193 169 END SUBROUTINE dom_vvl 194 170 171 172 SUBROUTINE dom_vvl_2( kt, pe3u_b, pe3v_b ) 173 !!---------------------------------------------------------------------- 174 !! *** ROUTINE dom_vvl_2 *** 175 !! 176 !! ** Purpose : compute the vertical scale factors at u- and v-points 177 !! in variable volume case. 178 !! 179 !! ** Method : In variable volume case (non linear sea surface) the 180 !! the vertical scale factor at velocity points is computed 181 !! as the average of the cell surface weighted e3t. 182 !! It uses the sea surface heigth so it have to be initialized 183 !! after ssh is read/set 184 !!---------------------------------------------------------------------- 185 INTEGER , INTENT(in ) :: kt ! ocean time-step index 186 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe3u_b, pe3v_b ! before vertical scale factor at u- & v-pts 187 ! 188 INTEGER :: ji, jj, jk ! dummy loop indices 189 INTEGER :: iku, ikv ! local integers 190 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers 191 REAL(wp) :: zvt ! local scalars 192 !!---------------------------------------------------------------------- 193 194 IF( lwp .AND. kt == nit000 ) THEN 195 WRITE(numout,*) 196 WRITE(numout,*) 'dom_vvl_2 : Variable volume, fse3t_b initialization' 197 WRITE(numout,*) '~~~~~~~~~ ' 198 pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk) 199 pe3v_b(:,:,jpk) = fse3u_0(:,:,jpk) 200 ENDIF 201 202 DO jk = 1, jpkm1 ! set the before scale factors at u- & v-points 203 DO jj = 2, jpjm1 204 DO ji = fs_2, fs_jpim1 205 zvt = fse3t_b(ji,jj,jk) * e1e2t(ji,jj) 206 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1e2t(ji+1,jj) ) / ( e1u(ji,jj) * e2u(ji,jj) ) 207 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e1e2t(ji,jj+1) ) / ( e1v(ji,jj) * e2v(ji,jj) ) 208 END DO 209 END DO 210 END DO 211 212 ! Correct scale factors at locations that have been individually modified in domhgr 213 ! Such modifications break the relationship between e1e2t and e1u*e2u etc. Recompute 214 ! scale factors ignoring the modified metric. 215 ! ! ===================== 216 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 217 ! ! ===================== 218 IF( nn_cla == 0 ) THEN 219 ! 220 ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u was modified) 221 ij0 = 102 ; ij1 = 102 222 DO jk = 1, jpkm1 ! set the before scale factors at u-points 223 DO jj = mj0(ij0), mj1(ij1) 224 DO ji = mi0(ii0), mi1(ii1) 225 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 226 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 227 END DO 228 END DO 229 END DO 230 ! 231 ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u and e1v were modified) 232 ij0 = 88 ; ij1 = 88 233 DO jk = 1, jpkm1 ! set the before scale factors at u-points 234 DO jj = mj0(ij0), mj1(ij1) 235 DO ji = mi0(ii0), mi1(ii1) 236 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 237 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 238 END DO 239 END DO 240 END DO 241 DO jk = 1, jpkm1 ! set the before scale factors at v-points 242 DO jj = mj0(ij0), mj1(ij1) 243 DO ji = mi0(ii0), mi1(ii1) 244 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 245 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 246 END DO 247 END DO 248 END DO 249 ENDIF 250 251 ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u was modified) 252 ij0 = 116 ; ij1 = 116 253 DO jk = 1, jpkm1 ! set the before scale factors at u-points 254 DO jj = mj0(ij0), mj1(ij1) 255 DO ji = mi0(ii0), mi1(ii1) 256 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 257 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 258 END DO 259 END DO 260 END DO 261 ! 262 ENDIF 263 ! ! ===================== 264 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration 265 ! ! ===================== 266 267 ii0 = 281 ; ii1 = 282 ! Gibraltar Strait (e2u was modified) 268 ij0 = 200 ; ij1 = 200 269 DO jk = 1, jpkm1 ! set the before scale factors at u-points 270 DO jj = mj0(ij0), mj1(ij1) 271 DO ji = mi0(ii0), mi1(ii1) 272 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 273 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 274 END DO 275 END DO 276 END DO 277 278 ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u was modified) 279 ij0 = 208 ; ij1 = 208 280 DO jk = 1, jpkm1 ! set the before scale factors at u-points 281 DO jj = mj0(ij0), mj1(ij1) 282 DO ji = mi0(ii0), mi1(ii1) 283 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 284 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 285 END DO 286 END DO 287 END DO 288 289 ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v was modified) 290 ij0 = 124 ; ij1 = 125 291 DO jk = 1, jpkm1 ! set the before scale factors at v-points 292 DO jj = mj0(ij0), mj1(ij1) 293 DO ji = mi0(ii0), mi1(ii1) 294 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 295 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 296 END DO 297 END DO 298 END DO 299 300 ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v was modified) [closed from bathy_11 on] 301 ij0 = 124 ; ij1 = 125 302 DO jk = 1, jpkm1 ! set the before scale factors at v-points 303 DO jj = mj0(ij0), mj1(ij1) 304 DO ji = mi0(ii0), mi1(ii1) 305 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 306 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 307 END DO 308 END DO 309 END DO 310 311 ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v was modified) 312 ij0 = 124 ; ij1 = 125 313 DO jk = 1, jpkm1 ! set the before scale factors at v-points 314 DO jj = mj0(ij0), mj1(ij1) 315 DO ji = mi0(ii0), mi1(ii1) 316 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 317 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 318 END DO 319 END DO 320 END DO 321 322 ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v was modified) 323 ij0 = 124 ; ij1 = 125 324 DO jk = 1, jpkm1 ! set the before scale factors at v-points 325 DO jj = mj0(ij0), mj1(ij1) 326 DO ji = mi0(ii0), mi1(ii1) 327 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 328 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 329 END DO 330 END DO 331 END DO 332 333 ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v was modified) 334 ij0 = 141 ; ij1 = 142 335 DO jk = 1, jpkm1 ! set the before scale factors at v-points 336 DO jj = mj0(ij0), mj1(ij1) 337 DO ji = mi0(ii0), mi1(ii1) 338 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 339 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 340 END DO 341 END DO 342 END DO 343 344 ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v was modified) 345 ij0 = 141 ; ij1 = 142 346 DO jk = 1, jpkm1 ! set the before scale factors at v-points 347 DO jj = mj0(ij0), mj1(ij1) 348 DO ji = mi0(ii0), mi1(ii1) 349 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 350 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 351 END DO 352 END DO 353 END DO 354 355 ! 356 ENDIF 357 ! ! ====================== 358 IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration 359 ! ! ====================== 360 ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u was modified) 361 ij0 = 327 ; ij1 = 327 362 DO jk = 1, jpkm1 ! set the before scale factors at u-points 363 DO jj = mj0(ij0), mj1(ij1) 364 DO ji = mi0(ii0), mi1(ii1) 365 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 366 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 367 END DO 368 END DO 369 END DO 370 ! 371 ii0 = 627 ; ii1 = 628 ! Bosphore Strait (e2u was modified) 372 ij0 = 343 ; ij1 = 343 373 DO jk = 1, jpkm1 ! set the before scale factors at u-points 374 DO jj = mj0(ij0), mj1(ij1) 375 DO ji = mi0(ii0), mi1(ii1) 376 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 377 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 378 END DO 379 END DO 380 END DO 381 ! 382 ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u was modified) 383 ij0 = 232 ; ij1 = 232 384 DO jk = 1, jpkm1 ! set the before scale factors at u-points 385 DO jj = mj0(ij0), mj1(ij1) 386 DO ji = mi0(ii0), mi1(ii1) 387 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 388 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 389 END DO 390 END DO 391 END DO 392 ! 393 ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u was modified) 394 ij0 = 232 ; ij1 = 232 395 DO jk = 1, jpkm1 ! set the before scale factors at u-points 396 DO jj = mj0(ij0), mj1(ij1) 397 DO ji = mi0(ii0), mi1(ii1) 398 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 399 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 400 END DO 401 END DO 402 END DO 403 ! 404 ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u was modified) 405 ij0 = 270 ; ij1 = 270 406 DO jk = 1, jpkm1 ! set the before scale factors at u-points 407 DO jj = mj0(ij0), mj1(ij1) 408 DO ji = mi0(ii0), mi1(ii1) 409 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 410 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 411 END DO 412 END DO 413 END DO 414 ! 415 ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v was modified) 416 ij0 = 232 ; ij1 = 233 417 DO jk = 1, jpkm1 ! set the before scale factors at v-points 418 DO jj = mj0(ij0), mj1(ij1) 419 DO ji = mi0(ii0), mi1(ii1) 420 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 421 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 422 END DO 423 END DO 424 END DO 425 ! 426 ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v was modified) 427 ij0 = 276 ; ij1 = 276 428 DO jk = 1, jpkm1 ! set the before scale factors at v-points 429 DO jj = mj0(ij0), mj1(ij1) 430 DO ji = mi0(ii0), mi1(ii1) 431 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 432 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 433 END DO 434 END DO 435 END DO 436 ! 437 ENDIF 438 ! End of individual corrections to scale factors 439 440 IF( ln_zps ) THEN ! minimum of the e3t at partial cell level 441 DO jj = 2, jpjm1 442 DO ji = fs_2, fs_jpim1 443 iku = mbku(ji,jj) 444 ikv = mbkv(ji,jj) 445 pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj ,iku) ) 446 pe3v_b(ji,jj,ikv) = MIN( fse3t_b(ji,jj,ikv), fse3t_b(ji ,jj+1,ikv) ) 447 END DO 448 END DO 449 ENDIF 450 451 pe3u_b(:,:,:) = pe3u_b(:,:,:) - fse3u_0(:,:,:) ! anomaly to avoid zero along closed boundary/extra halos 452 pe3v_b(:,:,:) = pe3v_b(:,:,:) - fse3v_0(:,:,:) 453 CALL lbc_lnk( pe3u_b(:,:,:), 'U', 1. ) ! lateral boundary conditions 454 CALL lbc_lnk( pe3v_b(:,:,:), 'V', 1. ) 455 pe3u_b(:,:,:) = pe3u_b(:,:,:) + fse3u_0(:,:,:) ! recover the full scale factor 456 pe3v_b(:,:,:) = pe3v_b(:,:,:) + fse3v_0(:,:,:) 457 ! 458 END SUBROUTINE dom_vvl_2 459 195 460 #else 196 461 !!---------------------------------------------------------------------- … … 200 465 SUBROUTINE dom_vvl 201 466 END SUBROUTINE dom_vvl 467 SUBROUTINE dom_vvl_2(kdum, pudum, pvdum ) 468 USE par_kind 469 INTEGER , INTENT(in ) :: kdum 470 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pudum, pvdum 471 END SUBROUTINE dom_vvl_2 202 472 #endif 203 473 -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r2777 r3101 83 83 neuler = 1 ! Set time-step indicator at nit000 (leap-frog) 84 84 CALL rst_read ! Read the restart file 85 ! ! define e3u_b, e3v_b from e3t_b read in restart file 86 CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 85 87 CALL tra_swap ! swap 3D arrays (t,s) in a 4D array (ts) 86 88 CALL day_init ! model calendar (using both namelist and restart infos) … … 92 94 CALL day_init ! model calendar (using both namelist and restart infos) 93 95 ! ! Initialization of ocean to zero 94 ! before fields ! now fields 95 sshb (:,:) = 0.e0 ; sshn (:,:) = 0.e0 96 ub (:,:,:) = 0.e0 ; un (:,:,:) = 0.e0 97 vb (:,:,:) = 0.e0 ; vn (:,:,:) = 0.e0 98 rotb (:,:,:) = 0.e0 ; rotn (:,:,:) = 0.e0 99 hdivb(:,:,:) = 0.e0 ; hdivn(:,:,:) = 0.e0 96 ! before fields ! now fields 97 sshb (:,:) = 0._wp ; sshn (:,:) = 0._wp 98 ub (:,:,:) = 0._wp ; un (:,:,:) = 0._wp 99 vb (:,:,:) = 0._wp ; vn (:,:,:) = 0._wp 100 rotb (:,:,:) = 0._wp ; rotn (:,:,:) = 0._wp 101 hdivb(:,:,:) = 0._wp ; hdivn(:,:,:) = 0._wp 102 ! 103 ! ! define e3u_b, e3v_b from e3t_b initialized in domzgr 104 CALL dom_vvl_2( nit000, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 100 105 ! 101 106 IF( cp_cfg == 'eel' ) THEN -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
r2715 r3101 13 13 USE dom_oce ! ocean space and time domain variables 14 14 USE zdf_oce ! ocean vertical physics variables 15 USE zdfbfr ! ocean bottom friction variables 15 16 USE trdmod ! ocean active dynamics and tracers trends 16 17 USE trdmod_oce ! ocean variables trends … … 51 52 !!--------------------------------------------------------------------- 52 53 ! 53 zm1_2dt = - 1._wp / ( 2._wp * rdt ) 54 IF( .not. ln_bfrimp) THEN ! only for explicit bottom friction form 55 ! implicit bfr is implemented in dynzdf_imp 56 ! H. Liu, Sept. 2011 54 57 55 IF( l_trddyn ) THEN ! temporary save of ua and va trends 56 ztrduv(:,:,:,1) = ua(:,:,:) 57 ztrduv(:,:,:,2) = va(:,:,:) 58 ENDIF 58 zm1_2dt = - 1._wp / ( 2._wp * rdt ) 59 60 IF( l_trddyn ) THEN ! temporary save of ua and va trends 61 ztrduv(:,:,:,1) = ua(:,:,:) 62 ztrduv(:,:,:,2) = va(:,:,:) 63 ENDIF 64 59 65 60 66 # if defined key_vectopt_loop 61 DO jj = 1, 162 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)67 DO jj = 1, 1 68 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 63 69 # else 64 DO jj = 2, jpjm165 DO ji = 2, jpim170 DO jj = 2, jpjm1 71 DO ji = 2, jpim1 66 72 # endif 67 ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels68 ikbv = mbkv(ji,jj)69 !70 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)71 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu)72 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv)73 END DO74 END DO73 ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels 74 ikbv = mbkv(ji,jj) 75 ! 76 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 77 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) 78 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) 79 END DO 80 END DO 75 81 76 ! 77 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 78 ztrduv(:,:,:,1) = ua(:,:,:) - ztrduv(:,:,:,1) 79 ztrduv(:,:,:,2) = va(:,:,:) - ztrduv(:,:,:,2) 80 CALL trd_mod( ztrduv(:,:,:,1), ztrduv(:,:,:,2), jpdyn_trd_bfr, 'DYN', kt ) 81 ENDIF 82 ! ! print mean trends (used for debugging) 83 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr - Ua: ', mask1=umask, & 84 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 85 ! 82 ! 83 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 84 ztrduv(:,:,:,1) = ua(:,:,:) - ztrduv(:,:,:,1) 85 ztrduv(:,:,:,2) = va(:,:,:) - ztrduv(:,:,:,2) 86 CALL trd_mod( ztrduv(:,:,:,1), ztrduv(:,:,:,2), jpdyn_trd_bfr, 'DYN', kt ) 87 ENDIF 88 ! ! print mean trends (used for debugging) 89 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr - Ua: ', mask1=umask, & 90 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 91 ! 92 ENDIF ! end explicit bottom friction 86 93 END SUBROUTINE dyn_bfr 87 94 -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r2715 r3101 27 27 !! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial) 28 28 !! hpg_rot : s-coordinate (ROTated axes scheme) 29 !! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial) 29 30 !!---------------------------------------------------------------------- 30 31 USE oce ! ocean dynamics and tracers … … 52 53 LOGICAL , PUBLIC :: ln_hpg_djc = .FALSE. !: s-coordinate (Density Jacobian with Cubic polynomial) 53 54 LOGICAL , PUBLIC :: ln_hpg_rot = .FALSE. !: s-coordinate (ROTated axes scheme) 55 LOGICAL , PUBLIC :: ln_hpg_prj = .FALSE. !: s-coordinate (Pressure Jacobian scheme) 54 56 REAL(wp), PUBLIC :: rn_gamma = 0._wp !: weighting coefficient 55 57 LOGICAL , PUBLIC :: ln_dynhpg_imp = .FALSE. !: semi-implicite hpg flag 56 58 57 INTEGER :: nhpg = 0 ! = 0 to 6, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags)59 INTEGER :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) 58 60 59 61 !! * Substitutions … … 92 94 ENDIF 93 95 ! 94 SELECT CASE ( nhpg ) ! Hydr astatic pressure gradient computation96 SELECT CASE ( nhpg ) ! Hydrostatic pressure gradient computation 95 97 CASE ( 0 ) ; CALL hpg_zco ( kt ) ! z-coordinate 96 98 CASE ( 1 ) ; CALL hpg_zps ( kt ) ! z-coordinate plus partial steps (interpolation) … … 100 102 CASE ( 5 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial) 101 103 CASE ( 6 ) ; CALL hpg_rot ( kt ) ! s-coordinate (ROTated axes scheme) 104 CASE ( 7 ) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme) 102 105 END SELECT 103 106 ! … … 129 132 !! 130 133 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel, & 131 & ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, rn_gamma , ln_dynhpg_imp 134 & ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, ln_hpg_prj, & 135 & rn_gamma , ln_dynhpg_imp 132 136 !!---------------------------------------------------------------------- 133 137 ! … … 147 151 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 148 152 WRITE(numout,*) ' s-coord. (ROTated axes scheme) ln_hpg_rot = ', ln_hpg_rot 153 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj 149 154 WRITE(numout,*) ' weighting coeff. (wdj scheme) rn_gamma = ', rn_gamma 150 155 WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp 151 156 ENDIF 152 157 ! 153 IF( lk_vvl .AND. .NOT. ln_hpg_sco ) & 154 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl require the standard jacobian formulation hpg_sco') 158 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) ) & 159 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 160 & the standard jacobian formulation hpg_sco or & 161 & the pressure jacobian formulation hpg_prj') 155 162 ! 156 163 ! ! Set nhpg from ln_hpg_... flags … … 162 169 IF( ln_hpg_djc ) nhpg = 5 163 170 IF( ln_hpg_rot ) nhpg = 6 164 ! 165 ! ! Consitency check 171 IF( ln_hpg_prj ) nhpg = 7 172 ! 173 ! ! Consistency check 166 174 ioptio = 0 167 175 IF( ln_hpg_zco ) ioptio = ioptio + 1 … … 172 180 IF( ln_hpg_djc ) ioptio = ioptio + 1 173 181 IF( ln_hpg_rot ) ioptio = ioptio + 1 182 IF( ln_hpg_prj ) ioptio = ioptio + 1 174 183 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 175 184 ! … … 1005 1014 END SUBROUTINE hpg_rot 1006 1015 1016 1017 SUBROUTINE hpg_prj( kt ) 1018 !!--------------------------------------------------------------------- 1019 !! *** ROUTINE hpg_prj *** 1020 !! 1021 !! ** Method : s-coordinate case. 1022 !! A Pressure-Jacobian horizontal pressure gradient method 1023 !! based on the constrained cubic-spline interpolation for 1024 !! all vertical coordinate systems 1025 !! 1026 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 1027 !! - Save the trend (l_trddyn=T) 1028 !! 1029 !!---------------------------------------------------------------------- 1030 1031 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 1032 USE oce , ONLY: zdeptht => ta ! (ta,sa) used as 3D workspace 1033 USE oce , ONLY: zrhh => sa 1034 USE wrk_nemo, ONLY: zhpi => wrk_3d_3 1035 USE wrk_nemo, ONLY: zu => wrk_3d_4 1036 USE wrk_nemo, ONLY: zv => wrk_3d_5 1037 USE wrk_nemo, ONLY: fsp => wrk_3d_6 1038 USE wrk_nemo, ONLY: xsp => wrk_3d_7 1039 USE wrk_nemo, ONLY: asp => wrk_3d_8 1040 USE wrk_nemo, ONLY: bsp => wrk_3d_9 1041 USE wrk_nemo, ONLY: csp => wrk_3d_10 1042 USE wrk_nemo, ONLY: dsp => wrk_3d_11 1043 !! 1044 !!---------------------------------------------------------------------- 1045 !! 1046 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear 1047 INTEGER, INTENT(in) :: kt ! ocean time-step index 1048 !! 1049 INTEGER :: ji, jj, jk, jkk ! dummy loop indices 1050 REAL(wp) :: zcoef0, znad ! temporary scalars 1051 !! 1052 !! The local variables for the correction term 1053 INTEGER :: jk1, jis, jid, jjs, jjd 1054 REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps 1055 REAL(wp) :: zrhdt1 1056 REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 1057 INTEGER :: zbhitwe, zbhitns 1058 !!---------------------------------------------------------------------- 1059 1060 IF( wrk_in_use(3, 3,4,5,6,7,8,9,10,11) ) THEN 1061 CALL ctl_stop('dyn:hpg_prj: requested workspace arrays unavailable') ; RETURN 1062 ENDIF 1063 1064 IF( kt == nit000 ) THEN 1065 IF(lwp) WRITE(numout,*) 1066 IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 1067 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, cubic spline pressure Jacobian' 1068 ENDIF 1069 1070 !!---------------------------------------------------------------------- 1071 ! Local constant initialization 1072 zcoef0 = - grav 1073 znad = 0.0_wp 1074 IF( lk_vvl ) znad = 1._wp 1075 1076 ! Clean 3-D work arrays 1077 zhpi(:,:,:) = 0. 1078 1079 ! Preparing vertical density profile for hybrid-sco coordinate 1080 DO jj = 1, jpj 1081 DO ji = 1, jpi 1082 jk = mbathy(ji,jj) 1083 IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp 1084 ELSE IF(jk == 1) THEN; zrhh(ji,jj, jk+1:jpk) = rhd(ji,jj,jk) 1085 ELSE IF(jk < jpkm1) THEN 1086 DO jkk = jk+1, jpk 1087 zrhh(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1), & 1088 fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 1089 END DO 1090 ENDIF 1091 END DO 1092 END DO 1093 1094 DO jj = 1, jpj 1095 DO ji = 1, jpi 1096 zdeptht(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) 1097 zdeptht(ji,jj,1) = zdeptht(ji,jj,1) - sshn(ji,jj) * znad 1098 DO jk = 2, jpk 1099 zdeptht(ji,jj,jk) = zdeptht(ji,jj,jk-1) + fse3w(ji,jj,jk) 1100 END DO 1101 END DO 1102 END DO 1103 1104 DO jk = 1, jpkm1 1105 DO jj = 1, jpj 1106 DO ji = 1, jpi 1107 fsp(ji,jj,jk) = zrhh(ji,jj,jk) 1108 xsp(ji,jj,jk) = zdeptht(ji,jj,jk) 1109 END DO 1110 END DO 1111 END DO 1112 1113 ! Construct the vertical density profile with the 1114 ! constrained cubic spline interpolation 1115 CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type) 1116 1117 ! Calculate the hydrostatic pressure at T(ji,jj,1) 1118 DO jj = 2, jpj 1119 DO ji = 2, jpi 1120 zrhdt1 = zrhh(ji,jj,1) - interp3(zdeptht(ji,jj,1),asp(ji,jj,1), & 1121 bsp(ji,jj,1), csp(ji,jj,1), & 1122 dsp(ji,jj,1) ) * 0.5_wp * zdeptht(ji,jj,1) 1123 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0) ! no lighter than fresh water 1124 1125 ! assuming linear profile across the top half surface layer 1126 zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * zrhdt1 1127 END DO 1128 END DO 1129 1130 ! Calculate the pressure at T(ji,jj,2:jpkm1) 1131 DO jk = 2, jpkm1 1132 DO jj = 2, jpj 1133 DO ji = 2, jpi 1134 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 1135 integ2(zdeptht(ji,jj,jk-1), zdeptht(ji,jj,jk),& 1136 asp(ji,jj,jk-1), bsp(ji,jj,jk-1), & 1137 csp(ji,jj,jk-1), dsp(ji,jj,jk-1)) 1138 END DO 1139 END DO 1140 END DO 1141 1142 ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 1143 DO jj = 2, jpjm1 1144 DO ji = 2, jpim1 1145 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad) 1146 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad) 1147 END DO 1148 END DO 1149 1150 DO jk = 2, jpkm1 1151 DO jj = 2, jpjm1 1152 DO ji = 2, jpim1 1153 zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk) 1154 zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk) 1155 END DO 1156 END DO 1157 END DO 1158 1159 DO jk = 1, jpkm1 1160 DO jj = 2, jpjm1 1161 DO ji = 2, jpim1 1162 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk) 1163 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk) 1164 END DO 1165 END DO 1166 END DO 1167 1168 DO jk = 1, jpkm1 1169 DO jj = 2, jpjm1 1170 DO ji = 2, jpim1 1171 zpwes = 0._wp; zpwed = 0._wp 1172 zpnss = 0._wp; zpnsd = 0._wp 1173 zuijk = zu(ji,jj,jk) 1174 zvijk = zv(ji,jj,jk) 1175 1176 !!!!! for u equation 1177 IF( jk <= mbku(ji,jj) ) THEN 1178 IF( -zdeptht(ji+1,jj,mbku(ji,jj)) >= -zdeptht(ji,jj,mbku(ji,jj)) ) THEN 1179 jis = ji + 1; jid = ji 1180 ELSE 1181 jis = ji; jid = ji +1 1182 ENDIF 1183 1184 ! integrate the pressure on the shallow side 1185 jk1 = jk 1186 zbhitwe = 0 1187 DO WHILE ( -zdeptht(jis,jj,jk1) > zuijk ) 1188 IF( jk1 == mbku(ji,jj) ) THEN 1189 zbhitwe = 1 1190 EXIT 1191 ENDIF 1192 zdeps = MIN(zdeptht(jis,jj,jk1+1), -zuijk) 1193 zpwes = zpwes + & 1194 integ2(zdeptht(jis,jj,jk1), zdeps, & 1195 asp(jis,jj,jk1), bsp(jis,jj,jk1), & 1196 csp(jis,jj,jk1), dsp(jis,jj,jk1)) 1197 jk1 = jk1 + 1 1198 END DO 1199 1200 IF(zbhitwe == 1) THEN 1201 zuijk = -zdeptht(jis,jj,jk1) 1202 ENDIF 1203 1204 ! integrate the pressure on the deep side 1205 jk1 = jk 1206 zbhitwe = 0 1207 DO WHILE ( -zdeptht(jid,jj,jk1) < zuijk ) 1208 IF( jk1 == 1 ) THEN 1209 zbhitwe = 1 1210 EXIT 1211 ENDIF 1212 zdeps = MAX(zdeptht(jid,jj,jk1-1), -zuijk) 1213 zpwed = zpwed + & 1214 integ2(zdeps, zdeptht(jid,jj,jk1), & 1215 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 1216 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) 1217 jk1 = jk1 - 1 1218 END DO 1219 1220 IF( zbhitwe == 1 ) THEN 1221 zdeps = zdeptht(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad) 1222 zrhdt1 = zrhh(jid,jj,1) - interp3(zdeptht(jid,jj,1), asp(jid,jj,1), & 1223 bsp(jid,jj,1), csp(jid,jj,1), & 1224 dsp(jid,jj,1)) * zdeps 1225 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0) ! no lighter than fresh water 1226 zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 1227 ENDIF 1228 1229 ! update the momentum trends in u direction 1230 1231 zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 1232 IF( lk_vvl ) THEN 1233 zdpdx2 = zcoef0 / e1u(ji,jj) * & 1234 ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 1235 ELSE 1236 zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1237 ENDIF 1238 1239 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 1240 & umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1241 ENDIF 1242 1243 !!!!! for v equation 1244 IF( jk <= mbkv(ji,jj) ) THEN 1245 IF( -zdeptht(ji,jj+1,mbkv(ji,jj)) >= -zdeptht(ji,jj,mbkv(ji,jj)) ) THEN 1246 jjs = jj + 1; jjd = jj 1247 ELSE 1248 jjs = jj ; jjd = jj + 1 1249 ENDIF 1250 1251 ! integrate the pressure on the shallow side 1252 jk1 = jk 1253 zbhitns = 0 1254 DO WHILE ( -zdeptht(ji,jjs,jk1) > zvijk ) 1255 IF( jk1 == mbkv(ji,jj) ) THEN 1256 zbhitns = 1 1257 EXIT 1258 ENDIF 1259 zdeps = MIN(zdeptht(ji,jjs,jk1+1), -zvijk) 1260 zpnss = zpnss + & 1261 integ2(zdeptht(ji,jjs,jk1), zdeps, & 1262 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & 1263 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) 1264 jk1 = jk1 + 1 1265 END DO 1266 1267 IF(zbhitns == 1) THEN 1268 zvijk = -zdeptht(ji,jjs,jk1) 1269 ENDIF 1270 1271 ! integrate the pressure on the deep side 1272 jk1 = jk 1273 zbhitns = 0 1274 DO WHILE ( -zdeptht(ji,jjd,jk1) < zvijk ) 1275 IF( jk1 == 1 ) THEN 1276 zbhitns = 1 1277 EXIT 1278 ENDIF 1279 zdeps = MAX(zdeptht(ji,jjd,jk1-1), -zvijk) 1280 zpnsd = zpnsd + & 1281 integ2(zdeps, zdeptht(ji,jjd,jk1), & 1282 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 1283 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 1284 jk1 = jk1 - 1 1285 END DO 1286 1287 IF( zbhitns == 1 ) THEN 1288 zdeps = zdeptht(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad) 1289 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdeptht(ji,jjd,1), asp(ji,jjd,1), & 1290 bsp(ji,jjd,1), csp(ji,jjd,1), & 1291 dsp(ji,jjd,1) ) * zdeps 1292 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0) ! no lighter than fresh water 1293 zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps 1294 ENDIF 1295 1296 ! update the momentum trends in v direction 1297 1298 zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 1299 IF( lk_vvl ) THEN 1300 zdpdy2 = zcoef0 / e2v(ji,jj) * & 1301 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 1302 ELSE 1303 zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1304 ENDIF 1305 1306 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 1307 & vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1308 ENDIF 1309 1310 1311 END DO 1312 END DO 1313 END DO 1314 1315 ! 1316 IF( wrk_not_released(3, 3,4,5,6,7,8,9,10,11) ) & 1317 CALL ctl_stop('dyn:hpg_prj: failed to release workspace arrays') 1318 ! 1319 XXXXXXX 1320 END SUBROUTINE hpg_prj 1321 1322 SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 1323 !!---------------------------------------------------------------------- 1324 !! *** ROUTINE cspline *** 1325 !! 1326 !! ** Purpose : constrained cubic spline interpolation 1327 !! 1328 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 1329 !! Reference: K.W. Brodlie, A review of mehtods for curve and function 1330 !! drawing, 1980 1331 !! 1332 !!---------------------------------------------------------------------- 1333 IMPLICIT NONE 1334 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: fsp, xsp ! value and coordinate 1335 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of 1336 ! the interpoated function 1337 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 1338 ! 2: Linear 1339 1340 ! Local Variables 1341 INTEGER :: ji, jj, jk ! dummy loop indices 1342 INTEGER :: jpi, jpj, jpkm1 1343 REAL(wp) :: zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp 1344 REAL(wp) :: zdxtmp1, zdxtmp2, zalpha 1345 REAL(wp) :: zdf(size(fsp,3)) 1346 !!---------------------------------------------------------------------- 1347 1348 jpi = size(fsp,1) 1349 jpj = size(fsp,2) 1350 jpkm1 = size(fsp,3) - 1 1351 1352 ! Clean output arrays 1353 asp = 0.0_wp 1354 bsp = 0.0_wp 1355 csp = 0.0_wp 1356 dsp = 0.0_wp 1357 1358 DO ji = 1, jpi 1359 DO jj = 1, jpj 1360 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 1361 DO jk = 2, jpkm1-1 1362 zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 1363 zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1364 zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 1365 zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 1366 1367 zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 1368 1369 IF(zdf1 * zdf2 <= 0._wp) THEN 1370 zdf(jk) = 0._wp 1371 ELSE 1372 zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 1373 ENDIF 1374 END DO 1375 1376 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1377 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1378 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1379 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1380 & 0.5_wp * zdf(jpkm1 - 1) 1381 1382 DO jk = 1, jpkm1 - 1 1383 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1384 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1385 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1386 zddf1 = -2._wp * ztmp1 + ztmp2 1387 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1388 zddf2 = 2._wp * ztmp1 - ztmp2 1389 1390 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1391 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1392 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1393 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1394 & dsp(ji,jj,jk) * ( xsp(ji,jj,jk+1)**2 + & 1395 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk) + & 1396 & xsp(ji,jj,jk)**2 ) 1397 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) - & 1398 & csp(ji,jj,jk) * xsp(ji,jj,jk)**2 - & 1399 & dsp(ji,jj,jk) * xsp(ji,jj,jk)**3 1400 END DO 1401 1402 ELSE IF (polynomial_type == 2) THEN ! Linear 1403 1404 DO jk = 1, jpkm1-1 1405 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1406 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1407 1408 dsp(ji,jj,jk) = 0._wp 1409 csp(ji,jj,jk) = 0._wp 1410 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1411 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1412 END DO 1413 1414 ELSE 1415 CALL ctl_stop( 'invalid polynomial type in cspline' ) 1416 ENDIF 1417 1418 END DO 1419 END DO 1420 1421 END SUBROUTINE cspline 1422 1423 1424 FUNCTION interp1(x, xl, xr, fl, fr) RESULT(f) 1425 !!---------------------------------------------------------------------- 1426 !! *** ROUTINE interp1 *** 1427 !! 1428 !! ** Purpose : 1-d linear interpolation 1429 !! 1430 !! ** Method : 1431 !! interpolation is straight forward 1432 !! extrapolation is also permitted (no value limit) 1433 !! 1434 !! H.Liu, Jan 2009, POL 1435 !!---------------------------------------------------------------------- 1436 IMPLICIT NONE 1437 REAL(wp), INTENT(in) :: x, xl, xr, fl, fr 1438 REAL(wp) :: f ! result of the interpolation (extrapolation) 1439 REAL(wp) :: zdeltx 1440 !!---------------------------------------------------------------------- 1441 1442 zdeltx = xr - xl 1443 IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN 1444 f = 0.5_wp * (fl + fr) 1445 ELSE 1446 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 1447 ENDIF 1448 1449 END FUNCTION interp1 1450 1451 FUNCTION interp2(x, a, b, c, d) RESULT(f) 1452 !!---------------------------------------------------------------------- 1453 !! *** ROUTINE interp1 *** 1454 !! 1455 !! ** Purpose : 1-d constrained cubic spline interpolation 1456 !! 1457 !! ** Method : cubic spline interpolation 1458 !! 1459 !!---------------------------------------------------------------------- 1460 IMPLICIT NONE 1461 REAL(wp), INTENT(in) :: x, a, b, c, d 1462 REAL(wp) :: f ! value from the interpolation 1463 !!---------------------------------------------------------------------- 1464 1465 f = a + x* ( b + x * ( c + d * x ) ) 1466 1467 END FUNCTION interp2 1468 1469 1470 FUNCTION interp3(x, a, b, c, d) RESULT(f) 1471 !!---------------------------------------------------------------------- 1472 !! *** ROUTINE interp1 *** 1473 !! 1474 !! ** Purpose : Calculate the first order of deriavtive of 1475 !! a cubic spline function y=a+b*x+c*x^2+d*x^3 1476 !! 1477 !! ** Method : f=dy/dx=b+2*c*x+3*d*x^2 1478 !! 1479 !!---------------------------------------------------------------------- 1480 IMPLICIT NONE 1481 REAL(wp), INTENT(in) :: x, a, b, c, d 1482 REAL(wp) :: f ! value from the interpolation 1483 !!---------------------------------------------------------------------- 1484 1485 f = b + x * ( 2._wp * c + 3._wp * d * x) 1486 1487 END FUNCTION interp3 1488 1489 1490 FUNCTION integ2(xl, xr, a, b, c, d) RESULT(f) 1491 !!---------------------------------------------------------------------- 1492 !! *** ROUTINE interp1 *** 1493 !! 1494 !! ** Purpose : 1-d constrained cubic spline integration 1495 !! 1496 !! ** Method : integrate polynomial a+bx+cx^2+dx^3 from xl to xr 1497 !! 1498 !!---------------------------------------------------------------------- 1499 IMPLICIT NONE 1500 REAL(wp), INTENT(in) :: xl, xr, a, b, c, d 1501 REAL(wp) :: za1, za2, za3 1502 REAL(wp) :: f ! integration result 1503 !!---------------------------------------------------------------------- 1504 1505 za1 = 0.5_wp * b 1506 za2 = c / 3.0_wp 1507 za3 = 0.25_wp * d 1508 1509 f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & 1510 & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 1511 1512 END FUNCTION integ2 1513 1514 1007 1515 !!====================================================================== 1008 1516 END MODULE dynhpg -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r3094 r3101 92 92 !! un,vn now horizontal velocity of next time-step 93 93 !!---------------------------------------------------------------------- 94 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released95 94 USE oce , ONLY: ze3u_f => ta , ze3v_f => sa ! (ta,sa) used as 3D workspace 96 USE wrk_nemo, ONLY: zs_t => wrk_2d_1 , zs_u_1 => wrk_2d_2 , zs_v_1 => wrk_2d_397 95 ! 98 96 INTEGER, INTENT( in ) :: kt ! ocean time-step index 99 97 ! 100 98 INTEGER :: ji, jj, jk ! dummy loop indices 99 INTEGER :: iku, ikv ! local integers 101 100 #if ! defined key_dynspg_flt 102 101 REAL(wp) :: z2dt ! temporary scalar 103 102 #endif 104 REAL(wp) :: zue3a, zue3n, zue3b, zuf ! local scalars 105 REAL(wp) :: zve3a, zve3n, zve3b, zvf ! - - 106 REAL(wp) :: zec, zv_t_ij, zv_t_ip1j, zv_t_ijp1 103 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars 104 REAL(wp) :: zve3a, zve3n, zve3b, zvf ! - - 107 105 !!---------------------------------------------------------------------- 108 109 IF( wrk_in_use(2, 1,2,3) ) THEN110 CALL ctl_stop('dyn_nxt: requested workspace arrays unavailable') ; RETURN111 ENDIF112 106 113 107 IF( kt == nit000 ) THEN … … 215 209 ELSE ! Variable volume ! 216 210 ! ! ================! 217 ! Before scale factor at t-points 218 ! ------------------------------- 219 DO jk = 1, jpkm1 211 ! 212 DO jk = 1, jpkm1 ! Before scale factor at t-points 220 213 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) & 221 214 & + atfp * ( fse3t_b(:,:,jk) + fse3t_a(:,:,jk) & 222 & - 2.e0 * fse3t_n(:,:,jk) ) 223 ENDDO 224 ! Add volume filter correction only at the first level of t-point scale factors 225 zec = atfp * rdt / rau0 215 & - 2._wp * fse3t_n(:,:,jk) ) 216 END DO 217 zec = atfp * rdt / rau0 ! Add filter correction only at the 1st level of t-point scale factors 226 218 fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 227 ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations228 zs_t (:,:) = e1t(:,:) * e2t(:,:)229 zs_u_1(:,:) = 0.5 / ( e1u(:,:) * e2u(:,:) )230 zs_v_1(:,:) = 0.5 / ( e1v(:,:) * e2v(:,:) )231 219 ! 232 IF( ln_dynadv_vec ) THEN 233 ! Before scale factor at (u/v)-points 234 ! ----------------------------------- 235 ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 236 DO jk = 1, jpkm1 237 DO jj = 1, jpjm1 238 DO ji = 1, jpim1 239 zv_t_ij = zs_t(ji ,jj ) * fse3t_b(ji ,jj ,jk) 240 zv_t_ip1j = zs_t(ji+1,jj ) * fse3t_b(ji+1,jj ,jk) 241 zv_t_ijp1 = zs_t(ji ,jj+1) * fse3t_b(ji ,jj+1,jk) 242 fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 243 fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 244 END DO 245 END DO 246 END DO 247 ! lateral boundary conditions 248 CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. ) 249 CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 250 ! Add initial scale factor to scale factor anomaly 251 fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 252 fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 253 ! Leap-Frog - Asselin filter and swap: applied on velocity 254 ! ----------------------------------- 255 DO jk = 1, jpkm1 256 DO jj = 1, jpj 220 IF( ln_dynadv_vec ) THEN ! vector invariant form (no thickness weighted calulation) 221 ! 222 ! ! before scale factors at u- & v-pts (computed from fse3t_b) 223 CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 224 ! 225 DO jk = 1, jpkm1 ! Leap-Frog - Asselin filter and swap: applied on velocity 226 DO jj = 1, jpj ! -------- 257 227 DO ji = 1, jpi 258 228 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) ) … … 267 237 END DO 268 238 ! 269 ELSE 270 ! Temporary filered scale factor at (u/v)-points (will become before scale factor) 271 !----------------------------------------------- 272 ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 273 DO jk = 1, jpkm1 274 DO jj = 1, jpjm1 275 DO ji = 1, jpim1 276 zv_t_ij = zs_t(ji ,jj ) * fse3t_b(ji ,jj ,jk) 277 zv_t_ip1j = zs_t(ji+1,jj ) * fse3t_b(ji+1,jj ,jk) 278 zv_t_ijp1 = zs_t(ji ,jj+1) * fse3t_b(ji ,jj+1,jk) 279 ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 280 ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 281 END DO 282 END DO 283 END DO 284 ! lateral boundary conditions 285 CALL lbc_lnk( ze3u_f, 'U', 1. ) 286 CALL lbc_lnk( ze3v_f, 'V', 1. ) 287 ! Add initial scale factor to scale factor anomaly 288 ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:) 289 ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:) 290 ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 291 ! ----------------------------------- =========================== 292 DO jk = 1, jpkm1 293 DO jj = 1, jpj 294 DO ji = 1, jpim1 239 ELSE ! flux form (thickness weighted calulation) 240 ! 241 CALL dom_vvl_2( kt, ze3u_f, ze3v_f ) ! before scale factors at u- & v-pts (computed from fse3t_b) 242 ! 243 DO jk = 1, jpkm1 ! Leap-Frog - Asselin filter and swap: 244 DO jj = 1, jpj ! applied on thickness weighted velocity 245 DO ji = 1, jpim1 ! --------------------------- 295 246 zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk) 296 247 zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk) … … 300 251 zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk) 301 252 ! 302 zuf = ( zue3n + atfp * ( zue3b - 2.e0* zue3n + zue3a ) ) / ze3u_f(ji,jj,jk)303 zvf = ( zve3n + atfp * ( zve3b - 2.e0* zve3n + zve3a ) ) / ze3v_f(ji,jj,jk)253 zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) 254 zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 304 255 ! 305 ub(ji,jj,jk) = zuf 256 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 306 257 vb(ji,jj,jk) = zvf 307 un(ji,jj,jk) = ua(ji,jj,jk) 258 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 308 259 vn(ji,jj,jk) = va(ji,jj,jk) 309 260 END DO 310 261 END DO 311 262 END DO 312 fse3u_b(:,:, :) = ze3u_f(:,:,:)! e3u_b <-- filtered scale factor313 fse3v_b(:,:, :) = ze3v_f(:,:,:)314 CALL lbc_lnk( ub, 'U', -1. ) 263 fse3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor 264 fse3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) 265 CALL lbc_lnk( ub, 'U', -1. ) ! lateral boundary conditions 315 266 CALL lbc_lnk( vb, 'V', -1. ) 316 267 ENDIF … … 323 274 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 324 275 ! 325 IF( wrk_not_released(2, 1,2,3) ) CALL ctl_stop('dyn_nxt: failed to release workspace arrays')326 !327 276 END SUBROUTINE dyn_nxt 328 277 -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r3094 r3101 118 118 INTEGER :: ji, jj, jk, jn ! dummy loop indices 119 119 INTEGER :: icycle ! local scalar 120 REAL(wp) :: zraur, zcoef, z2dt_e, z2dt_b ! local scalars 121 REAL(wp) :: z1_8, zx1, zy1 ! - - 122 REAL(wp) :: z1_4, zx2, zy2 ! - - 123 REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - - 124 REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - - 120 INTEGER :: ikbu, ikbv ! local scalar 121 REAL(wp) :: zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf ! local scalars 122 REAL(wp) :: z1_8, zx1, zy1 ! - - 123 REAL(wp) :: z1_4, zx2, zy2 ! - - 124 REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - - 125 REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - - 126 REAL(wp) :: ua_btm, va_btm ! - - 125 127 !!---------------------------------------------------------------------- 126 128 … … 146 148 hvr_e (:,:) = hvr (:,:) 147 149 IF( ln_dynvor_een ) THEN 148 ftne(1,:) = 0. e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0150 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 149 151 DO jj = 2, jpj 150 152 DO ji = fs_2, jpi ! vector opt. 151 ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3. 152 ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3. 153 ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3. 154 ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3. 153 ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3._wp 154 ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3._wp 155 ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3._wp 156 ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3._wp 155 157 END DO 156 158 END DO … … 159 161 ENDIF 160 162 161 ! !* Local constant initialization 162 z2dt_b = 2.0 * rdt ! baroclinic time step 163 z1_8 = 0.5 * 0.25 ! coefficient for vorticity estimates 164 z1_4 = 0.5 * 0.5 165 zraur = 1. / rau0 ! 1 / volumic mass 166 ! 167 zhdiv(:,:) = 0.e0 ! barotropic divergence 168 zu_sld = 0.e0 ; zu_asp = 0.e0 ! tides trends (lk_tide=F) 169 zv_sld = 0.e0 ; zv_asp = 0.e0 163 ! !* Local constant initialization 164 z1_2dt_b = 1._wp / ( 2.0_wp * rdt ) ! reciprocal of baroclinic time step 165 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt_b = 1.0_wp / rdt ! reciprocal of baroclinic 166 ! time step (euler timestep) 167 z1_8 = 0.125_wp ! coefficient for vorticity estimates 168 z1_4 = 0.25_wp 169 zraur = 1._wp / rau0 ! 1 / volumic mass 170 ! 171 zhdiv(:,:) = 0._wp ! barotropic divergence 172 zu_sld = 0._wp ; zu_asp = 0._wp ! tides trends (lk_tide=F) 173 zv_sld = 0._wp ; zv_asp = 0._wp 174 175 IF( kt == nit000 .AND. neuler == 0) THEN ! for implicit bottom friction 176 z2dt_bf = rdt 177 ELSE 178 z2dt_bf = 2.0_wp * rdt 179 ENDIF 170 180 171 181 ! ----------------------------------------------------------------------------- … … 175 185 ! !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 176 186 ! ! -------------------------- 177 zua(:,:) = 0. e0 ; zun(:,:) = 0.e0 ; ub_b(:,:) = 0.e0178 zva(:,:) = 0. e0 ; zvn(:,:) = 0.e0 ; vb_b(:,:) = 0.e0187 zua(:,:) = 0._wp ; zun(:,:) = 0._wp ; ub_b(:,:) = 0._wp 188 zva(:,:) = 0._wp ; zvn(:,:) = 0._wp ; vb_b(:,:) = 0._wp 179 189 ! 180 190 DO jk = 1, jpkm1 … … 194 204 ! 195 205 #if defined key_vvl 196 ub_b(ji,jj) = ub_b(ji,jj) + (fse3u_0(ji,jj,jk)*(1.+sshu_b(ji,jj)*muu(ji,jj,jk)))* ub(ji,jj,jk)197 vb_b(ji,jj) = vb_b(ji,jj) + (fse3v_0(ji,jj,jk)*(1.+sshv_b(ji,jj)*muv(ji,jj,jk)))* vb(ji,jj,jk)206 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk) *umask(ji,jj,jk) 207 vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk) *vmask(ji,jj,jk) 198 208 #else 199 209 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) … … 269 279 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 270 280 DO ji = fs_2, fs_jpim1 271 zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)272 zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)273 END DO281 zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 282 zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) 283 END DO 274 284 END DO 275 285 … … 277 287 ! ! Remove barotropic contribution of bottom friction 278 288 ! ! from the barotropic transport trend 279 zcoef = -1. / z2dt_b 289 zcoef = -1._wp * z1_2dt_b 290 291 IF(ln_bfrimp) THEN 292 ! ! Remove the bottom stress trend from 3-D sea surface level gradient 293 ! ! and Coriolis forcing in case of 3D semi-implicit bottom friction 294 DO jj = 2, jpjm1 295 DO ji = fs_2, fs_jpim1 296 ikbu = mbku(ji,jj) 297 ikbv = mbkv(ji,jj) 298 ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu) 299 va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv) 300 301 zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 302 zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 303 END DO 304 END DO 305 306 ELSE 307 280 308 # if defined key_vectopt_loop 281 DO jj = 1, 1282 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)309 DO jj = 1, 1 310 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 283 311 # else 284 DO jj = 2, jpjm1285 DO ji = 2, jpim1312 DO jj = 2, jpjm1 313 DO ji = 2, jpim1 286 314 # endif 287 315 ! Apply stability criteria for bottom friction 288 316 !RBbug for vvl and external mode we may need to use varying fse3 289 317 !!gm Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. 290 zbfru(ji,jj) = MAX( bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef ) 291 zbfrv(ji,jj) = MAX( bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef ) 292 END DO 293 END DO 294 295 IF( lk_vvl ) THEN 296 DO jj = 2, jpjm1 297 DO ji = fs_2, fs_jpim1 ! vector opt. 298 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) & 299 & / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) ) 300 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) & 301 & / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) ) 302 END DO 303 END DO 304 ELSE 305 DO jj = 2, jpjm1 306 DO ji = fs_2, fs_jpim1 ! vector opt. 307 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 308 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 309 END DO 310 END DO 311 ENDIF 312 318 zbfru(ji,jj) = MAX( bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef ) 319 zbfrv(ji,jj) = MAX( bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef ) 320 END DO 321 END DO 322 323 IF( lk_vvl ) THEN 324 DO jj = 2, jpjm1 325 DO ji = fs_2, fs_jpim1 ! vector opt. 326 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) & 327 & / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) ) 328 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) & 329 & / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) ) 330 END DO 331 END DO 332 ELSE 333 DO jj = 2, jpjm1 334 DO ji = fs_2, fs_jpim1 ! vector opt. 335 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 336 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 337 END DO 338 END DO 339 ENDIF 340 END IF ! end (ln_bfrimp) 341 342 313 343 ! !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 314 344 ! ! -------------------------- … … 317 347 ! 318 348 IF( lk_vvl ) THEN 319 ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1. e0- umask(:,:,1) )320 vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1. e0- vmask(:,:,1) )349 ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 350 vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 321 351 ELSE 322 352 ub_b(:,:) = ub_b(:,:) * hur(:,:) … … 354 384 ! set ssh corrections to 0 355 385 ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 356 IF( lp_obc_east ) sshfoe_b(:,:) = 0. e0357 IF( lp_obc_west ) sshfow_b(:,:) = 0. e0358 IF( lp_obc_south ) sshfos_b(:,:) = 0. e0359 IF( lp_obc_north ) sshfon_b(:,:) = 0. e0386 IF( lp_obc_east ) sshfoe_b(:,:) = 0._wp 387 IF( lp_obc_west ) sshfow_b(:,:) = 0._wp 388 IF( lp_obc_south ) sshfos_b(:,:) = 0._wp 389 IF( lp_obc_north ) sshfon_b(:,:) = 0._wp 360 390 #endif 361 391 … … 373 403 ! !* after ssh_e 374 404 ! ! ----------- 375 DO jj = 2, jpjm1 405 DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports 376 406 DO ji = fs_2, fs_jpim1 ! vector opt. 377 407 zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) & … … 385 415 ! ! OBC : zhdiv must be zero behind the open boundary 386 416 !! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 387 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0. e0! east388 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0. e0! west389 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0. e0! north390 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0. e0! south417 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0._wp ! east 418 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0._wp ! west 419 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0._wp ! north 420 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0._wp ! south 391 421 #endif 392 422 #if defined key_bdy … … 402 432 ! !* after barotropic velocities (vorticity scheme dependent) 403 433 ! ! --------------------------- 404 zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) 434 zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! now_e transport 405 435 zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 406 436 ! … … 426 456 zv_cor =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 427 457 ! after velocities with implicit bottom friction 428 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 429 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 430 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 431 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 458 459 IF( ln_bfrimp ) THEN ! implicit bottom friction 460 ! A new method to implement the implicit bottom friction. 461 ! H. Liu 462 ! Sept 2011 463 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 464 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 465 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 466 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 467 ! 468 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 469 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 470 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 471 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 472 ! 473 ELSE 474 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 475 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 476 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 477 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 478 ENDIF 432 479 END DO 433 480 END DO … … 452 499 zv_cor = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj) 453 500 ! after velocities with implicit bottom friction 454 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 455 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 456 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 457 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 501 IF( ln_bfrimp ) THEN 502 ! A new method to implement the implicit bottom friction. 503 ! H. Liu 504 ! Sept 2011 505 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 506 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 507 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 508 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 509 ! 510 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 511 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 512 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 513 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 514 ! 515 ELSE 516 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 517 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 518 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 519 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 520 ENDIF 458 521 END DO 459 522 END DO … … 478 541 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj) 479 542 ! after velocities with implicit bottom friction 480 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 481 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 482 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 483 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 543 IF( ln_bfrimp ) THEN 544 ! A new method to implement the implicit bottom friction. 545 ! H. Liu 546 ! Sept 2011 547 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 548 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 549 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 550 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 551 ! 552 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 553 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 554 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 555 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 556 ! 557 ELSE 558 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 559 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 560 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 561 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 562 ENDIF 484 563 END DO 485 564 END DO 486 565 ! 487 566 ENDIF 488 ! !* domain lateral boundary489 ! ! -----------------------567 ! !* domain lateral boundary 568 ! ! ----------------------- 490 569 491 570 ! OBC open boundaries … … 534 613 DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points 535 614 DO ji = 1, fs_jpim1 ! Vector opt. 536 zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &537 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) &538 & +e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )539 zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &540 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) &541 & +e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )615 zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 616 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & 617 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 618 zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 619 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & 620 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 542 621 END DO 543 622 END DO … … 547 626 hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points 548 627 hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 549 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1. e0- umask(:,:,1) )550 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1. e0- vmask(:,:,1) )628 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 629 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 551 630 ! 552 631 ENDIF … … 567 646 ! 568 647 ! !* Time average ==> after barotropic u, v, ssh 569 zcoef = 1. e0/ ( 2 * nn_baro + 1 )648 zcoef = 1._wp / ( 2 * nn_baro + 1 ) 570 649 zu_sum(:,:) = zcoef * zu_sum (:,:) 571 650 zv_sum(:,:) = zcoef * zv_sum (:,:) … … 573 652 ! !* update the general momentum trend 574 653 DO jk=1,jpkm1 575 ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b576 va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b654 ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b 655 va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b 577 656 END DO 578 657 un_b (:,:) = zu_sum(:,:) … … 608 687 CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 609 688 ELSE 610 un_b (:,:) = 0. e0611 vn_b (:,:) = 0. e0689 un_b (:,:) = 0._wp 690 vn_b (:,:) = 0._wp 612 691 ! vertical sum 613 692 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll … … 630 709 ! Vertically integrated velocity (before) 631 710 IF (neuler/=0) THEN 632 ub_b (:,:) = 0. e0633 vb_b (:,:) = 0. e0711 ub_b (:,:) = 0._wp 712 vb_b (:,:) = 0._wp 634 713 635 714 ! vertical sum … … 649 728 650 729 IF( lk_vvl ) THEN 651 ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1. e0- umask(:,:,1) )652 vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1. e0- vmask(:,:,1) )730 ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 731 vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 653 732 ELSE 654 733 ub_b(:,:) = ub_b(:,:) * hur(:,:) -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r2715 r3101 20 20 USE in_out_manager ! I/O manager 21 21 USE lib_mpp ! MPP library 22 USE zdfbfr ! bottom friction setup 22 23 23 24 IMPLICIT NONE … … 61 62 REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step 62 63 !! 63 INTEGER :: ji, jj, jk ! dummy loop indices 64 REAL(wp) :: z1_p2dt, zcoef, zzwi, zzws, zrhs ! local scalars 64 INTEGER :: ji, jj, jk ! dummy loop indices 65 INTEGER :: ikbum1, ikbvm1 ! local variable 66 REAL(wp) :: z1_p2dt, z2dtf, zcoef, zzwi, zzws, zrhs ! local scalars 67 68 !! * Local variables for implicit bottom friction. H. Liu 69 REAL(wp) :: zbfru, zbfrv 70 REAL(wp) :: zbfr_imp = 0._wp ! toggle (SAVE'd by assignment) 71 !!---------------------------------------------------------------------- 65 72 !!---------------------------------------------------------------------- 66 73 … … 73 80 IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 74 81 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 82 IF(ln_bfrimp) zbfr_imp = 1._wp 75 83 ENDIF 76 84 … … 80 88 81 89 ! 1. Vertical diffusion on u 90 91 ! Vertical diffusion on u&v 82 92 ! --------------------------- 83 93 ! Matrix and second member construction 84 ! bottom boundary condition: both zwi and zws must be masked as avmu can take 85 ! non zero value at the ocean bottom depending on the bottom friction 86 ! used but the bottom velocities have already been updated with the bottom 87 ! friction velocity in dyn_bfr using values from the previous timestep. There 88 ! is no need to include these in the implicit calculation. 89 ! 90 DO jk = 1, jpkm1 ! Matrix 91 DO jj = 2, jpjm1 92 DO ji = fs_2, fs_jpim1 ! vector opt. 94 !! bottom boundary condition: both zwi and zws must be masked as avmu can take 95 !! non zero value at the ocean bottom depending on the bottom friction 96 !! used but the bottom velocities have already been updated with the bottom 97 !! friction velocity in dyn_bfr using values from the previous timestep. There 98 !! is no need to include these in the implicit calculation. 99 100 ! The code has been modified here to implicitly implement bottom 101 ! friction: u(v)mask is not necessary here anymore. 102 ! H. Liu, April 2010. 103 104 ! 1. Vertical diffusion on u 105 DO jj = 2, jpjm1 106 DO ji = fs_2, fs_jpim1 ! vector opt. 107 ikbum1 = mbku(ji,jj) 108 zbfru = bfrua(ji,jj) 109 110 DO jk = 1, ikbum1 93 111 zcoef = - p2dt / fse3u(ji,jj,jk) 94 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 95 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) 96 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 97 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1) 98 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 99 END DO 100 END DO 101 END DO 102 DO jj = 2, jpjm1 ! Surface boudary conditions 103 DO ji = fs_2, fs_jpim1 ! vector opt. 112 zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk ) / fse3uw(ji,jj,jk ) 113 zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 114 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 115 END DO 116 117 ! Surface boundary conditions 104 118 zwi(ji,jj,1) = 0._wp 105 119 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 106 END DO 107 END DO 120 121 ! Bottom boundary conditions ! H. Liu, May, 2010 122 ! !commented out to be consistent with v3.2, h.liu 123 ! z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0_wp * zbfr_imp 124 z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * zbfr_imp 125 zws(ji,jj,ikbum1) = 0._wp 126 zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf 108 127 109 128 ! Matrix inversion starting from the first level … … 121 140 ! The solution (the after velocity) is in ua 122 141 !----------------------------------------------------------------------- 123 ! 124 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 125 DO jj = 2, jpjm1 126 DO ji = fs_2, fs_jpim1 ! vector opt. 142 143 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 144 DO jk = 2, ikbum1 127 145 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 128 146 END DO 129 END DO 130 END DO 131 ! 132 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 133 DO ji = fs_2, fs_jpim1 ! vector opt. 134 ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ( ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 135 & / ( fse3u(ji,jj,1) * rau0 ) ) 136 END DO 137 END DO 138 DO jk = 2, jpkm1 139 DO jj = 2, jpjm1 140 DO ji = fs_2, fs_jpim1 ! vector opt. 147 148 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 149 z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 ) 150 ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj)) 151 DO jk = 2, ikbum1 141 152 zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) ! zrhs=right hand side 142 153 ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 143 154 END DO 144 END DO 145 END DO 146 ! 147 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk == 148 DO ji = fs_2, fs_jpim1 ! vector opt. 149 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 150 END DO 151 END DO 152 DO jk = jpk-2, 1, -1 153 DO jj = 2, jpjm1 154 DO ji = fs_2, fs_jpim1 ! vector opt. 155 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 155 156 157 ! third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 158 ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1) 159 DO jk = ikbum1-1, 1, -1 160 ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 156 161 END DO 157 162 END DO … … 170 175 ! 2. Vertical diffusion on v 171 176 ! --------------------------- 172 ! Matrix and second member construction 173 ! bottom boundary condition: both zwi and zws must be masked as avmv can take 174 ! non zero value at the ocean bottom depending on the bottom friction 175 ! used but the bottom velocities have already been updated with the bottom 176 ! friction velocity in dyn_bfr using values from the previous timestep. There 177 ! is no need to include these in the implicit calculation. 178 ! 179 DO jk = 1, jpkm1 ! Matrix 177 178 DO ji = fs_2, fs_jpim1 ! vector opt. 180 179 DO jj = 2, jpjm1 181 DO ji = fs_2, fs_jpim1 ! vector opt. 180 ikbvm1 = mbkv(ji,jj) 181 zbfrv = bfrva(ji,jj) 182 183 DO jk = 1, ikbvm1 182 184 zcoef = -p2dt / fse3v(ji,jj,jk) 183 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk ) 184 zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk) 185 zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 186 zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1) 187 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 188 END DO 189 END DO 190 END DO 191 DO jj = 2, jpjm1 ! Surface boudary conditions 192 DO ji = fs_2, fs_jpim1 ! vector opt. 185 zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk ) / fse3vw(ji,jj,jk ) 186 zws(ji,jj,jk) = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 187 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 188 END DO 189 190 ! Surface boundary conditions 193 191 zwi(ji,jj,1) = 0._wp 194 192 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 195 END DO 196 END DO 193 194 ! Bottom boundary conditions ! H. Liu, May, 2010 195 ! !commented out to be consistent with v3.2, h.liu 196 ! z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0_wp * zbfr_imp 197 z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * zbfr_imp 198 zws(ji,jj,ikbvm1) = 0._wp 199 zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf 197 200 198 201 ! Matrix inversion … … 206 209 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 207 210 ! 208 ! m is decomposed in the product of an upper and lower triangular matrix 211 ! m is decomposed in the product of an upper and lower triangular 212 ! matrix 209 213 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 210 214 ! The solution (after velocity) is in 2d array va 211 215 !----------------------------------------------------------------------- 212 ! 213 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 214 DO jj = 2, jpjm1 215 DO ji = fs_2, fs_jpim1 ! vector opt. 216 217 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 218 DO jk = 2, ikbvm1 216 219 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 217 220 END DO 218 END DO 219 END DO 220 ! 221 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 222 DO ji = fs_2, fs_jpim1 ! vector opt. 223 va(ji,jj,1) = vb(ji,jj,1) + p2dt * ( va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 224 & / ( fse3v(ji,jj,1) * rau0 ) ) 225 END DO 226 END DO 227 DO jk = 2, jpkm1 228 DO jj = 2, jpjm1 229 DO ji = fs_2, fs_jpim1 ! vector opt. 221 222 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 223 z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 ) 224 va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj)) 225 DO jk = 2, ikbvm1 230 226 zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) ! zrhs=right hand side 231 227 va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 232 228 END DO 233 END DO 234 END DO 235 ! 236 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk == 237 DO ji = fs_2, fs_jpim1 ! vector opt. 238 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 239 END DO 240 END DO 241 DO jk = jpk-2, 1, -1 242 DO jj = 2, jpjm1 243 DO ji = fs_2, fs_jpim1 ! vector opt. 244 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 245 END DO 229 230 ! third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk 231 va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1) 232 233 DO jk = ikbvm1-1, 1, -1 234 va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 235 END DO 236 246 237 END DO 247 238 END DO … … 258 249 IF( wrk_not_released(3, 3) ) CALL ctl_stop('dyn_zdf_imp: failed to release workspace array') 259 250 ! 251 260 252 END SUBROUTINE dyn_zdf_imp 261 253 -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/LBC/lbcnfd.F90
r2715 r3101 236 236 END DO 237 237 END DO 238 CASE ( 'J' ) ! first ice U-V point 239 DO jl =0, ipr2dj 240 pt2d(2,ijpj+jl) = psgn * pt2d(3,ijpj-1+jl) 241 DO ji = 3, jpiglo 242 iju = jpiglo - ji + 3 243 pt2d(ji,ijpj+jl) = psgn * pt2d(iju,ijpj-1-jl) 244 END DO 245 END DO 246 CASE ( 'K' ) ! second ice U-V point 247 DO jl =0, ipr2dj 248 pt2d(2,ijpj+jl) = psgn * pt2d(3,ijpj-1+jl) 249 DO ji = 3, jpiglo 250 iju = jpiglo - ji + 3 251 pt2d(ji,ijpj+jl) = psgn * pt2d(iju,ijpj-1-jl) 252 END DO 253 END DO 238 254 END SELECT 239 255 ! … … 285 301 END DO 286 302 END DO 303 CASE ( 'J' ) ! first ice U-V point 304 pt2d( 2 ,ijpj:ijpj+ipr2dj) = 0.e0 305 DO jl = 0, ipr2dj 306 DO ji = 2 , jpiglo-1 307 ijt = jpiglo - ji + 2 308 pt2d(ji,ijpj+jl)= pt2d(ji,ijpj-1-jl) 309 END DO 310 END DO 311 CASE ( 'K' ) ! second ice U-V point 312 pt2d( 2 ,ijpj:ijpj+ipr2dj) = 0.e0 313 DO jl = 0, ipr2dj 314 DO ji = 2 , jpiglo-1 315 ijt = jpiglo - ji + 2 316 pt2d(ji,ijpj+jl)= pt2d(ijt,ijpj-1-jl) 317 END DO 318 END DO 287 319 END SELECT 288 320 ! … … 298 330 pt2d(:, 1:1-ipr2dj ) = 0.e0 299 331 pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0 332 CASE ( 'J' ) ! first ice U-V point 333 pt2d(:, 1:1-ipr2dj ) = 0.e0 334 pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0 335 CASE ( 'K' ) ! second ice U-V point 336 pt2d(:, 1:1-ipr2dj ) = 0.e0 337 pt2d(:,ijpj:ijpj+ipr2dj) = 0.e0 300 338 END SELECT 301 339 ! -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
r3096 r3101 164 164 REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE, SAVE :: ztab, znorthloc 165 165 REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: znorthgloio 166 REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE, SAVE :: zfoldwk ! Workspace for message transfers avoiding mpi_allgather 166 167 167 168 ! Arrays used in mpp_lbc_north_2d() 168 169 REAL(wp), DIMENSION(:,:) , ALLOCATABLE, SAVE :: ztab_2d, znorthloc_2d 169 170 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: znorthgloio_2d 171 REAL(wp), DIMENSION(:,:) , ALLOCATABLE, SAVE :: zfoldwk_2d ! Workspace for message transfers avoiding mpi_allgather 170 172 171 173 ! Arrays used in mpp_lbc_north_e() … … 173 175 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: znorthgloio_e 174 176 177 ! North fold arrays used to minimise the use of allgather operations. Set in nemo_northcomms (nemogcm) so need to be public 178 INTEGER, PUBLIC, PARAMETER :: jpmaxngh = 8 ! Assumed maximum number of active neighbours 179 INTEGER, PUBLIC, PARAMETER :: jptyps = 5 ! Number of different neighbour lists to be used for northfold exchanges 180 INTEGER, PUBLIC, DIMENSION (jpmaxngh,jptyps) :: isendto 181 INTEGER, PUBLIC, DIMENSION (jptyps) :: nsndto 182 LOGICAL, PUBLIC :: ln_nnogather = .FALSE. ! namelist control of northfold comms 183 LOGICAL, PUBLIC :: l_north_nogather = .FALSE. ! internal control of northfold comms 184 INTEGER, PUBLIC :: ityp 175 185 !!---------------------------------------------------------------------- 176 186 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 203 213 ! 204 214 & ztab(jpiglo,4,jpk) , znorthloc(jpi,4,jpk) , znorthgloio(jpi,4,jpk,jpni) , & 215 & zfoldwk(jpi,4,jpk) , & 205 216 ! 206 217 & ztab_2d(jpiglo,4) , znorthloc_2d(jpi,4) , znorthgloio_2d(jpi,4,jpni) , & 218 & zfoldwk_2d(jpi,4) , & 207 219 ! 208 220 & ztab_e(jpiglo,4+2*jpr2dj) , znorthloc_e(jpi,4+2*jpr2dj) , znorthgloio_e(jpi,4+2*jpr2dj,jpni) , & … … 232 244 LOGICAL :: mpi_was_called 233 245 ! 234 NAMELIST/nammpp/ cn_mpi_send, nn_buffer, jpni, jpnj, jpnij 246 NAMELIST/nammpp/ cn_mpi_send, nn_buffer, jpni, jpnj, jpnij, ln_nnogather 235 247 !!---------------------------------------------------------------------- 236 248 ! … … 269 281 WRITE(ldtxt(ii),*) ' number of local domains jpnij = ',jpnij; ii = ii +1 270 282 END IF 283 284 WRITE(ldtxt(ii),*) ' avoid use of mpi_allgather at the north fold ln_nnogather = ', ln_nnogather ; ii = ii + 1 271 285 272 286 CALL mpi_initialized ( mpi_was_called, code ) … … 441 455 CASE ( -1 ) 442 456 CALL mppsend( 2, t3we(1,1,1,1), imigr, noea, ml_req1 ) 443 CALL mpprecv( 1, t3ew(1,1,1,2), imigr )457 CALL mpprecv( 1, t3ew(1,1,1,2), imigr, noea ) 444 458 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 445 459 CASE ( 0 ) 446 460 CALL mppsend( 1, t3ew(1,1,1,1), imigr, nowe, ml_req1 ) 447 461 CALL mppsend( 2, t3we(1,1,1,1), imigr, noea, ml_req2 ) 448 CALL mpprecv( 1, t3ew(1,1,1,2), imigr )449 CALL mpprecv( 2, t3we(1,1,1,2), imigr )462 CALL mpprecv( 1, t3ew(1,1,1,2), imigr, noea ) 463 CALL mpprecv( 2, t3we(1,1,1,2), imigr, nowe ) 450 464 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 451 465 IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err) 452 466 CASE ( 1 ) 453 467 CALL mppsend( 1, t3ew(1,1,1,1), imigr, nowe, ml_req1 ) 454 CALL mpprecv( 2, t3we(1,1,1,2), imigr )468 CALL mpprecv( 2, t3we(1,1,1,2), imigr, nowe ) 455 469 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 456 470 END SELECT … … 494 508 CASE ( -1 ) 495 509 CALL mppsend( 4, t3sn(1,1,1,1), imigr, nono, ml_req1 ) 496 CALL mpprecv( 3, t3ns(1,1,1,2), imigr )510 CALL mpprecv( 3, t3ns(1,1,1,2), imigr, nono ) 497 511 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 498 512 CASE ( 0 ) 499 513 CALL mppsend( 3, t3ns(1,1,1,1), imigr, noso, ml_req1 ) 500 514 CALL mppsend( 4, t3sn(1,1,1,1), imigr, nono, ml_req2 ) 501 CALL mpprecv( 3, t3ns(1,1,1,2), imigr )502 CALL mpprecv( 4, t3sn(1,1,1,2), imigr )515 CALL mpprecv( 3, t3ns(1,1,1,2), imigr, nono ) 516 CALL mpprecv( 4, t3sn(1,1,1,2), imigr, noso ) 503 517 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 504 518 IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err) 505 519 CASE ( 1 ) 506 520 CALL mppsend( 3, t3ns(1,1,1,1), imigr, noso, ml_req1 ) 507 CALL mpprecv( 4, t3sn(1,1,1,2), imigr )521 CALL mpprecv( 4, t3sn(1,1,1,2), imigr, noso ) 508 522 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 509 523 END SELECT … … 635 649 CASE ( -1 ) 636 650 CALL mppsend( 2, t2we(1,1,1), imigr, noea, ml_req1 ) 637 CALL mpprecv( 1, t2ew(1,1,2), imigr )651 CALL mpprecv( 1, t2ew(1,1,2), imigr, noea ) 638 652 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 639 653 CASE ( 0 ) 640 654 CALL mppsend( 1, t2ew(1,1,1), imigr, nowe, ml_req1 ) 641 655 CALL mppsend( 2, t2we(1,1,1), imigr, noea, ml_req2 ) 642 CALL mpprecv( 1, t2ew(1,1,2), imigr )643 CALL mpprecv( 2, t2we(1,1,2), imigr )656 CALL mpprecv( 1, t2ew(1,1,2), imigr, noea ) 657 CALL mpprecv( 2, t2we(1,1,2), imigr, nowe ) 644 658 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 645 659 IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err) 646 660 CASE ( 1 ) 647 661 CALL mppsend( 1, t2ew(1,1,1), imigr, nowe, ml_req1 ) 648 CALL mpprecv( 2, t2we(1,1,2), imigr )662 CALL mpprecv( 2, t2we(1,1,2), imigr, nowe ) 649 663 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 650 664 END SELECT … … 688 702 CASE ( -1 ) 689 703 CALL mppsend( 4, t2sn(1,1,1), imigr, nono, ml_req1 ) 690 CALL mpprecv( 3, t2ns(1,1,2), imigr )704 CALL mpprecv( 3, t2ns(1,1,2), imigr, nono ) 691 705 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 692 706 CASE ( 0 ) 693 707 CALL mppsend( 3, t2ns(1,1,1), imigr, noso, ml_req1 ) 694 708 CALL mppsend( 4, t2sn(1,1,1), imigr, nono, ml_req2 ) 695 CALL mpprecv( 3, t2ns(1,1,2), imigr )696 CALL mpprecv( 4, t2sn(1,1,2), imigr )709 CALL mpprecv( 3, t2ns(1,1,2), imigr, nono ) 710 CALL mpprecv( 4, t2sn(1,1,2), imigr, noso ) 697 711 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 698 712 IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err) 699 713 CASE ( 1 ) 700 714 CALL mppsend( 3, t2ns(1,1,1), imigr, noso, ml_req1 ) 701 CALL mpprecv( 4, t2sn(1,1,2), imigr )715 CALL mpprecv( 4, t2sn(1,1,2), imigr, noso ) 702 716 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 703 717 END SELECT … … 816 830 CASE ( -1 ) 817 831 CALL mppsend( 2, t4we(1,1,1,1,1), imigr, noea, ml_req1 ) 818 CALL mpprecv( 1, t4ew(1,1,1,1,2), imigr )832 CALL mpprecv( 1, t4ew(1,1,1,1,2), imigr, noea ) 819 833 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 820 834 CASE ( 0 ) 821 835 CALL mppsend( 1, t4ew(1,1,1,1,1), imigr, nowe, ml_req1 ) 822 836 CALL mppsend( 2, t4we(1,1,1,1,1), imigr, noea, ml_req2 ) 823 CALL mpprecv( 1, t4ew(1,1,1,1,2), imigr )824 CALL mpprecv( 2, t4we(1,1,1,1,2), imigr )837 CALL mpprecv( 1, t4ew(1,1,1,1,2), imigr, noea ) 838 CALL mpprecv( 2, t4we(1,1,1,1,2), imigr, nowe ) 825 839 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 826 840 IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err) 827 841 CASE ( 1 ) 828 842 CALL mppsend( 1, t4ew(1,1,1,1,1), imigr, nowe, ml_req1 ) 829 CALL mpprecv( 2, t4we(1,1,1,1,2), imigr )843 CALL mpprecv( 2, t4we(1,1,1,1,2), imigr, nowe ) 830 844 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 831 845 END SELECT … … 875 889 CASE ( -1 ) 876 890 CALL mppsend( 4, t4sn(1,1,1,1,1), imigr, nono, ml_req1 ) 877 CALL mpprecv( 3, t4ns(1,1,1,1,2), imigr )891 CALL mpprecv( 3, t4ns(1,1,1,1,2), imigr, nono ) 878 892 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 879 893 CASE ( 0 ) 880 894 CALL mppsend( 3, t4ns(1,1,1,1,1), imigr, noso, ml_req1 ) 881 895 CALL mppsend( 4, t4sn(1,1,1,1,1), imigr, nono, ml_req2 ) 882 CALL mpprecv( 3, t4ns(1,1,1,1,2), imigr )883 CALL mpprecv( 4, t4sn(1,1,1,1,2), imigr )896 CALL mpprecv( 3, t4ns(1,1,1,1,2), imigr, nono ) 897 CALL mpprecv( 4, t4sn(1,1,1,1,2), imigr, noso ) 884 898 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 885 899 IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err) 886 900 CASE ( 1 ) 887 901 CALL mppsend( 3, t4ns(1,1,1,1,1), imigr, noso, ml_req1 ) 888 CALL mpprecv( 4, t4sn(1,1,1,1,2), imigr )902 CALL mpprecv( 4, t4sn(1,1,1,1,2), imigr, noso ) 889 903 IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err) 890 904 END SELECT … … 1019 1033 CASE ( -1 ) 1020 1034 CALL mppsend( 2, tr2we(1-jpr2dj,1,1), imigr, noea, ml_req1 ) 1021 CALL mpprecv( 1, tr2ew(1-jpr2dj,1,2), imigr )1035 CALL mpprecv( 1, tr2ew(1-jpr2dj,1,2), imigr, noea ) 1022 1036 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1023 1037 CASE ( 0 ) 1024 1038 CALL mppsend( 1, tr2ew(1-jpr2dj,1,1), imigr, nowe, ml_req1 ) 1025 1039 CALL mppsend( 2, tr2we(1-jpr2dj,1,1), imigr, noea, ml_req2 ) 1026 CALL mpprecv( 1, tr2ew(1-jpr2dj,1,2), imigr )1027 CALL mpprecv( 2, tr2we(1-jpr2dj,1,2), imigr )1040 CALL mpprecv( 1, tr2ew(1-jpr2dj,1,2), imigr, noea ) 1041 CALL mpprecv( 2, tr2we(1-jpr2dj,1,2), imigr, nowe ) 1028 1042 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1029 1043 IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err) 1030 1044 CASE ( 1 ) 1031 1045 CALL mppsend( 1, tr2ew(1-jpr2dj,1,1), imigr, nowe, ml_req1 ) 1032 CALL mpprecv( 2, tr2we(1-jpr2dj,1,2), imigr )1046 CALL mpprecv( 2, tr2we(1-jpr2dj,1,2), imigr, nowe ) 1033 1047 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1034 1048 END SELECT … … 1072 1086 CASE ( -1 ) 1073 1087 CALL mppsend( 4, tr2sn(1-jpr2di,1,1), imigr, nono, ml_req1 ) 1074 CALL mpprecv( 3, tr2ns(1-jpr2di,1,2), imigr )1088 CALL mpprecv( 3, tr2ns(1-jpr2di,1,2), imigr, nono ) 1075 1089 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1076 1090 CASE ( 0 ) 1077 1091 CALL mppsend( 3, tr2ns(1-jpr2di,1,1), imigr, noso, ml_req1 ) 1078 1092 CALL mppsend( 4, tr2sn(1-jpr2di,1,1), imigr, nono, ml_req2 ) 1079 CALL mpprecv( 3, tr2ns(1-jpr2di,1,2), imigr )1080 CALL mpprecv( 4, tr2sn(1-jpr2di,1,2), imigr )1093 CALL mpprecv( 3, tr2ns(1-jpr2di,1,2), imigr, nono ) 1094 CALL mpprecv( 4, tr2sn(1-jpr2di,1,2), imigr, noso ) 1081 1095 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1082 1096 IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err) 1083 1097 CASE ( 1 ) 1084 1098 CALL mppsend( 3, tr2ns(1-jpr2di,1,1), imigr, noso, ml_req1 ) 1085 CALL mpprecv( 4, tr2sn(1-jpr2di,1,2), imigr )1099 CALL mpprecv( 4, tr2sn(1-jpr2di,1,2), imigr, noso ) 1086 1100 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err) 1087 1101 END SELECT … … 1138 1152 1139 1153 1140 SUBROUTINE mpprecv( ktyp, pmess, kbytes )1154 SUBROUTINE mpprecv( ktyp, pmess, kbytes, ksource ) 1141 1155 !!---------------------------------------------------------------------- 1142 1156 !! *** routine mpprecv *** … … 1148 1162 INTEGER , INTENT(in ) :: kbytes ! suze of the array pmess 1149 1163 INTEGER , INTENT(in ) :: ktyp ! Tag of the recevied message 1164 INTEGER, OPTIONAL, INTENT(in) :: ksource ! source process number 1150 1165 !! 1151 1166 INTEGER :: istatus(mpi_status_size) 1152 1167 INTEGER :: iflag 1153 !!---------------------------------------------------------------------- 1154 ! 1155 CALL mpi_recv( pmess, kbytes, mpi_double_precision, mpi_any_source, ktyp, mpi_comm_opa, istatus, iflag ) 1168 INTEGER :: use_source 1169 !!---------------------------------------------------------------------- 1170 ! 1171 1172 ! If a specific process number has been passed to the receive call, 1173 ! use that one. Default is to use mpi_any_source 1174 use_source=mpi_any_source 1175 if(present(ksource)) then 1176 use_source=ksource 1177 end if 1178 1179 CALL mpi_recv( pmess, kbytes, mpi_double_precision, use_source, ktyp, mpi_comm_opa, istatus, iflag ) 1156 1180 ! 1157 1181 END SUBROUTINE mpprecv … … 1833 1857 IF( nbondi == -1 ) THEN 1834 1858 CALL mppsend( 2, t2we(1,1,1), imigr, noea, ml_req1 ) 1835 CALL mpprecv( 1, t2ew(1,1,2), imigr )1859 CALL mpprecv( 1, t2ew(1,1,2), imigr, noea ) 1836 1860 IF(l_isend) CALL mpi_wait( ml_req1, ml_stat, ml_err ) 1837 1861 ELSEIF( nbondi == 0 ) THEN 1838 1862 CALL mppsend( 1, t2ew(1,1,1), imigr, nowe, ml_req1 ) 1839 1863 CALL mppsend( 2, t2we(1,1,1), imigr, noea, ml_req2 ) 1840 CALL mpprecv( 1, t2ew(1,1,2), imigr )1841 CALL mpprecv( 2, t2we(1,1,2), imigr )1864 CALL mpprecv( 1, t2ew(1,1,2), imigr, noea ) 1865 CALL mpprecv( 2, t2we(1,1,2), imigr, nowe ) 1842 1866 IF(l_isend) CALL mpi_wait( ml_req1, ml_stat, ml_err ) 1843 1867 IF(l_isend) CALL mpi_wait( ml_req2, ml_stat, ml_err ) 1844 1868 ELSEIF( nbondi == 1 ) THEN 1845 1869 CALL mppsend( 1, t2ew(1,1,1), imigr, nowe, ml_req1 ) 1846 CALL mpprecv( 2, t2we(1,1,2), imigr )1870 CALL mpprecv( 2, t2we(1,1,2), imigr, nowe ) 1847 1871 IF(l_isend) CALL mpi_wait( ml_req1, ml_stat, ml_err ) 1848 1872 ENDIF … … 1879 1903 IF( nbondj == -1 ) THEN 1880 1904 CALL mppsend( 4, t2sn(1,1,1), imigr, nono, ml_req1 ) 1881 CALL mpprecv( 3, t2ns(1,1,2), imigr )1905 CALL mpprecv( 3, t2ns(1,1,2), imigr, nono ) 1882 1906 IF(l_isend) CALL mpi_wait( ml_req1, ml_stat, ml_err ) 1883 1907 ELSEIF( nbondj == 0 ) THEN 1884 1908 CALL mppsend( 3, t2ns(1,1,1), imigr, noso, ml_req1 ) 1885 1909 CALL mppsend( 4, t2sn(1,1,1), imigr, nono, ml_req2 ) 1886 CALL mpprecv( 3, t2ns(1,1,2), imigr )1887 CALL mpprecv( 4, t2sn(1,1,2), imigr )1910 CALL mpprecv( 3, t2ns(1,1,2), imigr, nono ) 1911 CALL mpprecv( 4, t2sn(1,1,2), imigr, noso ) 1888 1912 IF( l_isend ) CALL mpi_wait( ml_req1, ml_stat, ml_err ) 1889 1913 IF( l_isend ) CALL mpi_wait( ml_req2, ml_stat, ml_err ) 1890 1914 ELSEIF( nbondj == 1 ) THEN 1891 1915 CALL mppsend( 3, t2ns(1,1,1), imigr, noso, ml_req1 ) 1892 CALL mpprecv( 4, t2sn(1,1,2), imigr )1916 CALL mpprecv( 4, t2sn(1,1,2), imigr, noso) 1893 1917 IF( l_isend ) CALL mpi_wait( ml_req1, ml_stat, ml_err ) 1894 1918 ENDIF … … 2209 2233 INTEGER :: ierr, itaille, ildi, ilei, iilb 2210 2234 INTEGER :: ijpj, ijpjm1, ij, iproc 2235 INTEGER, DIMENSION (jpmaxngh) :: ml_req_nf ! for mpi_isend when avoiding mpi_allgather 2236 INTEGER :: ml_err ! for mpi_isend when avoiding mpi_allgather 2237 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat ! for mpi_isend when avoiding mpi_allgather 2211 2238 !!---------------------------------------------------------------------- 2212 2239 ! 2213 2240 ijpj = 4 2241 ityp = -1 2214 2242 ijpjm1 = 3 2215 2243 ztab(:,:,:) = 0.e0 … … 2222 2250 ! ! Build in procs of ncomm_north the znorthgloio 2223 2251 itaille = jpi * jpk * ijpj 2224 CALL MPI_ALLGATHER( znorthloc , itaille, MPI_DOUBLE_PRECISION, & 2225 & znorthgloio, itaille, MPI_DOUBLE_PRECISION, ncomm_north, ierr ) 2226 ! 2227 ! ! recover the global north array 2228 DO jr = 1, ndim_rank_north 2229 iproc = nrank_north(jr) + 1 2230 ildi = nldit (iproc) 2231 ilei = nleit (iproc) 2232 iilb = nimppt(iproc) 2233 DO jj = 1, 4 2234 DO ji = ildi, ilei 2235 ztab(ji+iilb-1,jj,:) = znorthgloio(ji,jj,:,jr) 2252 IF ( l_north_nogather ) THEN 2253 ! 2254 ! Avoid the use of mpi_allgather by exchanging only with the processes already identified 2255 ! (in nemo_northcomms) as being involved in this process' northern boundary exchange 2256 ! 2257 DO jj = nlcj-ijpj+1, nlcj ! First put local values into the global array 2258 ij = jj - nlcj + ijpj 2259 DO ji = 1, nlci 2260 ztab(ji+nimpp-1,ij,:) = pt3d(ji,jj,:) 2236 2261 END DO 2237 2262 END DO 2238 END DO 2263 2264 ! 2265 ! Set the exchange type in order to access the correct list of active neighbours 2266 ! 2267 SELECT CASE ( cd_type ) 2268 CASE ( 'T' , 'W' ) 2269 ityp = 1 2270 CASE ( 'U' ) 2271 ityp = 2 2272 CASE ( 'V' ) 2273 ityp = 3 2274 CASE ( 'F' ) 2275 ityp = 4 2276 CASE ( 'I' ) 2277 ityp = 5 2278 CASE DEFAULT 2279 ityp = -1 ! Set a default value for unsupported types which 2280 ! will cause a fallback to the mpi_allgather method 2281 END SELECT 2282 IF ( ityp .gt. 0 ) THEN 2283 2284 DO jr = 1,nsndto(ityp) 2285 CALL mppsend(5, znorthloc, itaille, isendto(jr,ityp), ml_req_nf(jr) ) 2286 END DO 2287 DO jr = 1,nsndto(ityp) 2288 CALL mpprecv(5, zfoldwk, itaille, isendto(jr,ityp)) 2289 iproc = isendto(jr,ityp) + 1 2290 ildi = nldit (iproc) 2291 ilei = nleit (iproc) 2292 iilb = nimppt(iproc) 2293 DO jj = 1, ijpj 2294 DO ji = ildi, ilei 2295 ztab(ji+iilb-1,jj,:) = zfoldwk(ji,jj,:) 2296 END DO 2297 END DO 2298 END DO 2299 IF (l_isend) THEN 2300 DO jr = 1,nsndto(ityp) 2301 CALL mpi_wait(ml_req_nf(jr), ml_stat, ml_err) 2302 END DO 2303 ENDIF 2304 2305 ENDIF 2306 2307 ENDIF 2308 2309 IF ( ityp .lt. 0 ) THEN 2310 CALL MPI_ALLGATHER( znorthloc , itaille, MPI_DOUBLE_PRECISION, & 2311 & znorthgloio, itaille, MPI_DOUBLE_PRECISION, ncomm_north, ierr ) 2312 ! 2313 DO jr = 1, ndim_rank_north ! recover the global north array 2314 iproc = nrank_north(jr) + 1 2315 ildi = nldit (iproc) 2316 ilei = nleit (iproc) 2317 iilb = nimppt(iproc) 2318 DO jj = 1, ijpj 2319 DO ji = ildi, ilei 2320 ztab(ji+iilb-1,jj,:) = znorthgloio(ji,jj,:,jr) 2321 END DO 2322 END DO 2323 END DO 2324 ENDIF 2325 ! 2326 ! The ztab array has been either: 2327 ! a. Fully populated by the mpi_allgather operation or 2328 ! b. Had the active points for this domain and northern neighbours populated 2329 ! by peer to peer exchanges 2330 ! Either way the array may be folded by lbc_nfd and the result for the span of 2331 ! this domain will be identical. 2239 2332 ! 2240 2333 CALL lbc_nfd( ztab, cd_type, psgn ) ! North fold boundary condition … … 2272 2365 INTEGER :: ierr, itaille, ildi, ilei, iilb 2273 2366 INTEGER :: ijpj, ijpjm1, ij, iproc 2367 INTEGER, DIMENSION (jpmaxngh) :: ml_req_nf ! for mpi_isend when avoiding mpi_allgather 2368 INTEGER :: ml_err ! for mpi_isend when avoiding mpi_allgather 2369 INTEGER, DIMENSION(MPI_STATUS_SIZE):: ml_stat ! for mpi_isend when avoiding mpi_allgather 2274 2370 !!---------------------------------------------------------------------- 2275 2371 ! 2276 2372 ijpj = 4 2373 ityp = -1 2277 2374 ijpjm1 = 3 2278 2375 ztab_2d(:,:) = 0.e0 … … 2285 2382 ! ! Build in procs of ncomm_north the znorthgloio_2d 2286 2383 itaille = jpi * ijpj 2287 CALL MPI_ALLGATHER( znorthloc_2d , itaille, MPI_DOUBLE_PRECISION, & 2288 & znorthgloio_2d, itaille, MPI_DOUBLE_PRECISION, ncomm_north, ierr ) 2289 ! 2290 DO jr = 1, ndim_rank_north ! recover the global north array 2291 iproc = nrank_north(jr) + 1 2292 ildi=nldit (iproc) 2293 ilei=nleit (iproc) 2294 iilb=nimppt(iproc) 2295 DO jj = 1, 4 2296 DO ji = ildi, ilei 2297 ztab_2d(ji+iilb-1,jj) = znorthgloio_2d(ji,jj,jr) 2384 IF ( l_north_nogather ) THEN 2385 ! 2386 ! Avoid the use of mpi_allgather by exchanging only with the processes already identified 2387 ! (in nemo_northcomms) as being involved in this process' northern boundary exchange 2388 ! 2389 DO jj = nlcj-ijpj+1, nlcj ! First put local values into the global array 2390 ij = jj - nlcj + ijpj 2391 DO ji = 1, nlci 2392 ztab_2d(ji+nimpp-1,ij) = pt2d(ji,jj) 2298 2393 END DO 2299 2394 END DO 2300 END DO 2395 2396 ! 2397 ! Set the exchange type in order to access the correct list of active neighbours 2398 ! 2399 SELECT CASE ( cd_type ) 2400 CASE ( 'T' , 'W' ) 2401 ityp = 1 2402 CASE ( 'U' ) 2403 ityp = 2 2404 CASE ( 'V' ) 2405 ityp = 3 2406 CASE ( 'F' ) 2407 ityp = 4 2408 CASE ( 'I' ) 2409 ityp = 5 2410 CASE DEFAULT 2411 ityp = -1 ! Set a default value for unsupported types which 2412 ! will cause a fallback to the mpi_allgather method 2413 END SELECT 2414 2415 IF ( ityp .gt. 0 ) THEN 2416 2417 DO jr = 1,nsndto(ityp) 2418 CALL mppsend(5, znorthloc_2d, itaille, isendto(jr,ityp), ml_req_nf(jr) ) 2419 END DO 2420 DO jr = 1,nsndto(ityp) 2421 CALL mpprecv(5, zfoldwk_2d, itaille, isendto(jr,ityp)) 2422 iproc = isendto(jr,ityp) + 1 2423 ildi = nldit (iproc) 2424 ilei = nleit (iproc) 2425 iilb = nimppt(iproc) 2426 DO jj = 1, ijpj 2427 DO ji = ildi, ilei 2428 ztab_2d(ji+iilb-1,jj) = zfoldwk_2d(ji,jj) 2429 END DO 2430 END DO 2431 END DO 2432 IF (l_isend) THEN 2433 DO jr = 1,nsndto(ityp) 2434 CALL mpi_wait(ml_req_nf(jr), ml_stat, ml_err) 2435 END DO 2436 ENDIF 2437 2438 ENDIF 2439 2440 ENDIF 2441 2442 IF ( ityp .lt. 0 ) THEN 2443 CALL MPI_ALLGATHER( znorthloc_2d , itaille, MPI_DOUBLE_PRECISION, & 2444 & znorthgloio_2d, itaille, MPI_DOUBLE_PRECISION, ncomm_north, ierr ) 2445 ! 2446 DO jr = 1, ndim_rank_north ! recover the global north array 2447 iproc = nrank_north(jr) + 1 2448 ildi = nldit (iproc) 2449 ilei = nleit (iproc) 2450 iilb = nimppt(iproc) 2451 DO jj = 1, ijpj 2452 DO ji = ildi, ilei 2453 ztab_2d(ji+iilb-1,jj) = znorthgloio_2d(ji,jj,jr) 2454 END DO 2455 END DO 2456 END DO 2457 ENDIF 2458 ! 2459 ! The ztab array has been either: 2460 ! a. Fully populated by the mpi_allgather operation or 2461 ! b. Had the active points for this domain and northern neighbours populated 2462 ! by peer to peer exchanges 2463 ! Either way the array may be folded by lbc_nfd and the result for the span of 2464 ! this domain will be identical. 2301 2465 ! 2302 2466 CALL lbc_nfd( ztab_2d, cd_type, psgn ) ! North fold boundary condition … … 2499 2663 2500 2664 LOGICAL, PUBLIC, PARAMETER :: lk_mpp = .FALSE. !: mpp flag 2665 LOGICAL, PUBLIC :: ln_nnogather = .FALSE. !: namelist control of northfold comms (needed here in case "key_mpp_mpi" is not used) 2501 2666 INTEGER :: ncomm_ice 2502 2667 !!---------------------------------------------------------------------- -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r2772 r3101 46 46 ! !! Griffies operator 47 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wslp2 !: wslp**2 from Griffies quarter cells 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi_g, triadj_g !: skew flux slopes relative to geopotentials 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi_g, triadj_g !: skew flux slopes relative to geopotentials 49 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi , triadj !: isoneutral slopes relative to model-coordinate 50 50 … … 58 58 59 59 ! Workspace arrays for ldf_slp_grif. These could be replaced by several 3D and 2D workspace 60 ! arrays from the wrk_nemo module with a bit of code re-writing. The 4D workspace 60 ! arrays from the wrk_nemo module with a bit of code re-writing. The 4D workspace 61 61 ! arrays can't be used here because of the zero-indexing of some of the ranks. ARPDBG. 62 62 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: zdzrho , zdyrho, zdxrho ! Horizontal and vertical density gradients … … 93 93 !!---------------------------------------------------------------------- 94 94 !! *** ROUTINE ldf_slp *** 95 !! 95 !! 96 96 !! ** Purpose : Compute the slopes of neutral surface (slope of isopycnal 97 97 !! surfaces referenced locally) (ln_traldf_iso=T). 98 !! 99 !! ** Method : The slope in the i-direction is computed at U- and 100 !! W-points (uslp, wslpi) and the slope in the j-direction is 98 !! 99 !! ** Method : The slope in the i-direction is computed at U- and 100 !! W-points (uslp, wslpi) and the slope in the j-direction is 101 101 !! computed at V- and W-points (vslp, wslpj). 102 102 !! They are bounded by 1/100 over the whole ocean, and within the … … 112 112 !! bottom slope (ln_sco=T) at level jpk in inildf] 113 113 !! 114 !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and j-slopes 114 !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and j-slopes 115 115 !! of now neutral surfaces at u-, w- and v- w-points, resp. 116 116 !!---------------------------------------------------------------------- … … 127 127 INTEGER :: ii0, ii1, iku ! temporary integer 128 128 INTEGER :: ij0, ij1, ikv ! temporary integer 129 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16 129 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16, zcofw ! local scalars 130 130 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 131 131 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - … … 148 148 DO jj = 1, jpjm1 149 149 DO ji = 1, fs_jpim1 ! vector opt. 150 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 151 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 150 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 151 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 152 152 END DO 153 153 END DO 154 154 END DO 155 155 IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level 156 # if defined key_vectopt_loop 156 # if defined key_vectopt_loop 157 157 DO jj = 1, 1 158 158 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) … … 161 161 DO ji = 1, jpim1 162 162 # endif 163 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 164 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 163 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 164 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 165 165 END DO 166 166 END DO … … 181 181 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml 182 182 183 183 184 184 ! I. slopes at u and v point | uslp = d/di( prd ) / d/dz( prd ) 185 185 ! =========================== | vslp = d/dj( prd ) / d/dz( prd ) 186 ! 186 ! 187 187 DO jk = 2, jpkm1 !* Slopes at u and v points 188 188 DO jj = 2, jpjm1 … … 225 225 DO jk = 2, jpkm1 226 226 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 227 DO ji = 2, jpim1 227 DO ji = 2, jpim1 228 228 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 229 229 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & … … 266 266 ! II. slopes at w point | wslpi = mij( d/di( prd ) / d/dz( prd ) 267 267 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 268 ! 268 ! 269 269 DO jk = 2, jpkm1 270 270 DO jj = 2, jpjm1 … … 308 308 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 309 309 DO ji = 2, jpim1 310 zcofw = tmask(ji,jj,jk) * z1_16 310 311 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 311 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &312 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &313 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &314 & + 4.* zwz(ji ,jj ,jk) ) * z1_16 * tmask(ji,jj,jk)312 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 313 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 314 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 315 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 315 316 316 317 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 317 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &318 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &319 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &320 & + 4.* zww(ji ,jj ,jk) ) * z1_16 * tmask(ji,jj,jk)321 END DO 322 END DO 318 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 319 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 320 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 321 & + 4.* zww(ji ,jj ,jk) ) * zcofw 322 END DO 323 END DO 323 324 DO jj = 3, jpj-2 ! other rows 324 325 DO ji = fs_2, fs_jpim1 ! vector opt. 326 zcofw = tmask(ji,jj,jk) * z1_16 325 327 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 326 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &327 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &328 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &329 & + 4.* zwz(ji ,jj ,jk) ) * z1_16 * tmask(ji,jj,jk)328 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 329 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 330 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 331 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 330 332 331 333 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 332 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &333 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &334 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &335 & + 4.* zww(ji ,jj ,jk) ) * z1_16 * tmask(ji,jj,jk)334 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 335 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 336 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 337 & + 4.* zww(ji ,jj ,jk) ) * zcofw 336 338 END DO 337 339 END DO … … 346 348 END DO 347 349 END DO 348 349 ! III. Specific grid points 350 ! =========================== 351 ! 350 351 ! III. Specific grid points 352 ! =========================== 353 ! 352 354 IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN ! ORCA_R4 configuration: horizontal diffusion in specific area 353 355 ! ! Gibraltar Strait … … 368 370 ENDIF 369 371 370 ! IV. Lateral boundary conditions 372 373 ! IV. Lateral boundary conditions 371 374 ! =============================== 372 375 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) … … 382 385 ! 383 386 END SUBROUTINE ldf_slp 384 387 385 388 386 389 SUBROUTINE ldf_slp_grif ( kt ) … … 390 393 !! ** Purpose : Compute the squared slopes of neutral surfaces (slope 391 394 !! of iso-pycnal surfaces referenced locally) (ln_traldf_grif=T) 392 !! at W-points using the Griffies quarter-cells. 393 !! 394 !! ** Method : calculates alpha and beta at T-points 395 !! at W-points using the Griffies quarter-cells. 396 !! 397 !! ** Method : calculates alpha and beta at T-points 395 398 !! 396 399 !! ** Action : - triadi_g, triadj_g T-pts i- and j-slope triads relative to geopot. (used for eiv) … … 399 402 !!---------------------------------------------------------------------- 400 403 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 401 USE oce , ONLY: zdit => ua , zdis => va ! (ua,va) used as workspace 402 USE oce , ONLY: zdjt => ta , zdjs => sa ! (ta,sa) used as workspace 403 USE wrk_nemo, ONLY: zdkt => wrk_3d_2 , zdks => wrk_3d_3 ! 3D workspace 404 USE wrk_nemo, ONLY: zalpha => wrk_3d_4 , zbeta => wrk_3d_5 ! alpha, beta at T points, at depth fsgdept 404 USE oce , ONLY: zalbet => ua ! use ua as workspace 405 405 USE wrk_nemo, ONLY: z1_mlbw => wrk_2d_1 406 ! 407 INTEGER, INTENT( in ) :: kt ! ocean time-step index408 ! 406 !! 407 INTEGER, INTENT( in ) :: kt ! ocean time-step index 408 !! 409 409 INTEGER :: ji, jj, jk, jl, ip, jp, kp ! dummy loop indices 410 INTEGER :: iku, ikv ! local integer 411 REAL(wp) :: zfacti, zfactj, zatempw,zatempu,zatempv ! local scalars 412 REAL(wp) :: zbu, zbv, zbti, zbtj ! - - 413 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_lim2, zti_g_raw, zti_g_lim 414 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_lim2, ztj_g_raw, ztj_g_lim 410 INTEGER :: iku, ikv ! local integer 411 REAL(wp) :: zfacti, zfactj ! local scalars 412 REAL(wp) :: znot_thru_surface ! local scalars 413 REAL(wp) :: zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj 414 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_g_raw, zti_g_lim 415 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim 415 416 REAL(wp) :: zdzrho_raw 417 REAL(wp) :: zbeta0 416 418 !!---------------------------------------------------------------------- 417 419 … … 424 426 !--------------------------------! 425 427 ! 426 CALL eos_alpbet( tsb, zalpha, zbeta ) !== before thermal and haline expension coeff. at T-points ==! 427 ! 428 DO jk = 1, jpkm1 !== before lateral T & S gradients at T-level jk ==! 429 DO jj = 1, jpjm1 430 DO ji = 1, fs_jpim1 ! vector opt. 431 zdit(ji,jj,jk) = ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) * umask(ji,jj,jk) ! i-gradient of T and S at jj 432 zdis(ji,jj,jk) = ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) * umask(ji,jj,jk) 433 zdjt(ji,jj,jk) = ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) * vmask(ji,jj,jk) ! j-gradient of T and S at jj 434 zdjs(ji,jj,jk) = ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) * vmask(ji,jj,jk) 435 END DO 436 END DO 437 END DO 438 IF( ln_zps ) THEN ! partial steps: correction at the last level 428 CALL eos_alpbet( tsb, zalbet, zbeta0 ) !== before local thermal/haline expension ratio at T-points ==! 429 ! 430 DO jl = 0, 1 !== unmasked before density i- j-, k-gradients ==! 431 ! 432 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln) 433 DO jk = 1, jpkm1 ! done each pair of triad 434 DO jj = 1, jpjm1 ! NB: not masked ==> a minimum value is set 435 DO ji = 1, fs_jpim1 ! vector opt. 436 zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! i-gradient of T & S at u-point 437 zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 438 zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! j-gradient of T & S at v-point 439 zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 440 zdxrho_raw = ( - zalbet(ji+ip,jj ,jk) * zdit + zbeta0*zdis ) / e1u(ji,jj) 441 zdyrho_raw = ( - zalbet(ji ,jj+jp,jk) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 442 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 443 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 444 END DO 445 END DO 446 END DO 447 ! 448 IF( ln_zps.and.l_grad_zps ) THEN ! partial steps: correction of i- & j-grad on bottom 439 449 # if defined key_vectopt_loop 440 DO jj = 1, 1441 DO ji = 1, jpij-jpi! vector opt. (forced unrolling)450 DO jj = 1, 1 451 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 442 452 # else 443 DO jj = 1, jpjm1444 DO ji = 1, jpim1453 DO jj = 1, jpjm1 454 DO ji = 1, jpim1 445 455 # endif 446 zdit(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_tem) ! i-gradient of T and S 447 zdis(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_sal) 448 zdjt(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_tem) ! j-gradient of T and S 449 zdjs(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_sal) 450 END DO 451 END DO 452 ENDIF 453 ! 454 zdkt(:,:,1) = 0._wp !== before vertical T & S gradient at w-level ==! 455 zdks(:,:,1) = 0._wp 456 DO jk = 2, jpk 457 zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk) 458 zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk) 459 END DO 460 ! 461 ! 462 DO jl = 0, 1 !== density i-, j-, and k-gradients ==! 463 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln) 464 DO jk = 1, jpkm1 ! done each pair of triad 465 DO jj = 1, jpjm1 ! NB: not masked due to the minimum value set 466 DO ji = 1, fs_jpim1 ! vector opt. 467 zdxrho_raw = ( zalpha(ji+ip,jj ,jk) * zdit(ji,jj,jk) + zbeta(ji+ip,jj ,jk) * zdis(ji,jj,jk) ) / e1u(ji,jj) 468 zdyrho_raw = ( zalpha(ji ,jj+jp,jk) * zdjt(ji,jj,jk) + zbeta(ji ,jj+jp,jk) * zdjs(ji,jj,jk) ) / e2v(ji,jj) 469 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 470 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 456 iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points) 457 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 458 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 459 zdxrho_raw = ( - zalbet(ji+ip,jj ,iku) * zdit + zbeta0*zdis ) / e1u(ji,jj) 460 zdyrho_raw = ( - zalbet(ji ,jj+jp,ikv) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 461 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 462 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 471 463 END DO 472 464 END DO 473 END DO 474 END DO 475 DO kp = 0, 1 !== density i-, j-, and k-gradients ==! 476 DO jk = 1, jpkm1 ! done each pair of triad 477 DO jj = 1, jpj ! NB: not masked due to the minimum value set 478 DO ji = 1, jpi ! vector opt. 479 zdzrho_raw = ( zalpha(ji,jj,jk) * zdkt(ji,jj,jk+kp) + zbeta(ji,jj,jk) * zdks(ji,jj,jk+kp) ) & 480 & / fse3w(ji,jj,jk+kp) 481 zdzrho(ji ,jj ,jk, kp) = - MIN( - repsln, zdzrho_raw ) ! force zdzrho >= repsln 482 END DO 483 END DO 484 END DO 485 END DO 486 ! 487 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==! 465 ENDIF 466 ! 467 END DO 468 469 DO kp = 0, 1 !== unmasked before density i- j-, k-gradients ==! 470 DO jk = 1, jpkm1 ! done each pair of triad 471 DO jj = 1, jpj ! NB: not masked ==> a minimum value is set 472 DO ji = 1, jpi ! vector opt. 473 IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp 474 zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) ) 475 zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) ) 476 ELSE 477 zdkt = 0._wp ! 1st level gradient set to zero 478 zdks = 0._wp 479 ENDIF 480 zdzrho_raw = ( - zalbet(ji ,jj ,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp) 481 zdzrho(ji ,jj ,jk, kp) = - MIN( - repsln, zdzrho_raw ) ! force zdzrho >= repsln 482 END DO 483 END DO 484 END DO 485 END DO 486 ! 487 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==! 488 488 DO ji = 1, jpi 489 489 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth … … 492 492 END DO 493 493 ! 494 ! !== intialisations to zero ==!495 ! 496 wslp2 (:,:,:) = 0._wp 497 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp 494 ! !== intialisations to zero ==! 495 ! 496 wslp2 (:,:,:) = 0._wp ! wslp2 will be cumulated 3D field set to zero 497 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 498 498 triadj_g(:,:,1,:,:) = 0._wp ; triadj_g(:,:,jpk,:,:) = 0._wp 499 !!gm _iso set to zero missing500 triadi (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp! set surface and bottom slope to zero501 triadj (:,:,1,:,:) = 0._wp ; triadj(:,:,jpk,:,:) = 0._wp502 499 !!gm _iso set to zero missing 500 triadi (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 501 triadj (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp 502 503 503 !-------------------------------------! 504 504 ! Triads just below the Mixed Layer ! 505 505 !-------------------------------------! 506 506 ! 507 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base508 DO kp = 0, 1 ! with only the slope-max limit and MASKED507 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 508 DO kp = 0, 1 ! with only the slope-max limit and MASKED 509 509 DO jj = 1, jpjm1 510 510 DO ji = 1, fs_jpim1 511 511 ip = jl ; jp = jl 512 512 jk = MIN( nmln(ji+ip,jj) , mbkt(ji+ip,jj) ) + 1 ! ML level+1 (MIN in case ML depth is the ocean depth) 513 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 513 514 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 514 & +( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj) ) * umask(ji,jj,jk)515 & - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj) ) * umask(ji,jj,jk) 515 516 jk = MIN( nmln(ji,jj+jp) , mbkt(ji,jj+jp) ) + 1 516 517 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 517 & +( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk)518 & - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 518 519 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 519 520 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) … … 527 528 !-------------------------------------! 528 529 ! 529 DO kp = 0, 1 ! k-index of triads530 DO kp = 0, 1 ! k-index of triads 530 531 DO jl = 0, 1 531 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes)532 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes) 532 533 DO jk = 1, jpkm1 534 ! Must mask contribution to slope from dz/dx at constant s for triads jk=1,kp=0 that poke up though ocean surface 535 znot_thru_surface = REAL( 1-1/(jk+kp), wp ) !jk+kp=1,=0.; otherwise=1.0 533 536 DO jj = 1, jpjm1 534 DO ji = 1, fs_jpim1 ! vector opt.537 DO ji = 1, fs_jpim1 ! vector opt. 535 538 ! 536 539 ! Calculate slope relative to geopotentials used for GM skew fluxes 537 ! For s-coordinate, subtract slope at t-points (equivalent to *adding* gradient of depth)540 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 538 541 ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 539 542 ! masked by umask taken at the level of dz(rho) … … 543 546 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 544 547 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 545 zti_coord = ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 546 ztj_coord = ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 547 ! unmasked 548 zti_g_raw = zti_raw + zti_coord ! ref to geopot surfaces 549 ztj_g_raw = ztj_raw + ztj_coord 548 549 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 550 zti_coord = znot_thru_surface * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 551 ztj_coord = znot_thru_surface * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) ! unmasked 552 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 553 ztj_g_raw = ztj_raw - ztj_coord 550 554 zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 551 555 ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 552 556 ! 553 ! Below ML use limited zti_g as is 554 ! Inside ML replace by linearly reducing sx_mlb towards surface 557 ! Below ML use limited zti_g as is & mask 558 ! Inside ML replace by linearly reducing sx_mlb towards surface & mask 555 559 ! 556 560 zfacti = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)), wp ) ! k index of uppermost point(s) of triad is jk+kp-1 557 561 zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1 558 562 ! ! otherwise zfact=0 559 zti_g_lim = zfacti * zti_g_lim &563 zti_g_lim = ( zfacti * zti_g_lim & 560 564 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 561 & * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) 562 ztj_g_lim = zfactj * ztj_g_lim &565 & * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 566 ztj_g_lim = ( zfactj * ztj_g_lim & 563 567 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 564 & * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ! masked565 ! 566 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim * umask(ji,jj,jk+kp)567 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim * vmask(ji,jj,jk+kp)568 & * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 569 ! 570 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim 571 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim 568 572 ! 569 573 ! Get coefficients of isoneutral diffusion tensor … … 574 578 ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 575 579 ! 576 zti_lim = zti_g_lim - zti_coord ! remove the coordinate slope ==> relative to coordinate surfaces 577 ztj_lim = ztj_g_lim - ztj_coord 578 zti_lim2 = zti_lim * zti_lim * umask(ji,jj,jk+kp) ! square of limited slopes ! masked <<== 579 ztj_lim2 = ztj_lim * ztj_lim * vmask(ji,jj,jk+kp) 580 ! 580 zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces 581 ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 582 ! 583 IF( ln_triad_iso ) THEN 584 zti_raw = ( zti_lim*zti_lim ) / zti_raw 585 ztj_raw = ( ztj_lim*ztj_lim ) / ztj_raw 586 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 587 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 588 zti_lim = zfacti * zti_lim & 589 & + ( 1._wp - zfacti ) * zti_raw 590 ztj_lim = zfactj * ztj_lim & 591 & + ( 1._wp - zfactj ) * ztj_raw 592 ENDIF 593 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim 594 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim 595 ! 581 596 zbu = e1u(ji ,jj) * e2u(ji ,jj) * fse3u(ji ,jj,jk ) 582 597 zbv = e1v(ji ,jj) * e2v(ji ,jj) * fse3v(ji ,jj,jk ) … … 584 599 zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 585 600 ! 586 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim2 / zti_raw ! masked 587 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim2 / ztj_raw 588 ! 589 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked 590 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2 ! masked 591 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2 601 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked 602 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * ( zti_g_lim * zti_g_lim ) ! masked 603 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ( ztj_g_lim * ztj_g_lim ) 592 604 END DO 593 605 END DO … … 597 609 ! 598 610 wslp2(:,:,1) = 0._wp ! force the surface wslp to zero 599 611 600 612 CALL lbc_lnk( wslp2, 'W', 1. ) ! lateral boundary confition on wslp2 only ==>>> gm : necessary ? to be checked 601 613 ! … … 610 622 !! *** ROUTINE ldf_slp_mxl *** 611 623 !! 612 !! ** Purpose : Compute the slopes of iso-neutral surface just below 624 !! ** Purpose : Compute the slopes of iso-neutral surface just below 613 625 !! the mixed layer. 614 626 !! … … 619 631 !! 620 632 !! ** Action : uslpml, wslpiml : i- & j-slopes of neutral surfaces 621 !! vslpml, wslpjml just below the mixed layer 633 !! vslpml, wslpjml just below the mixed layer 622 634 !! omlmask : mixed layer mask 623 635 !!---------------------------------------------------------------------- … … 627 639 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_dzr ! z-gradient of density (T-point) 628 640 !! 629 INTEGER :: ji , jj , jk ! dummy loop indices630 INTEGER :: iku, ikv, ik, ikm1 ! local integers641 INTEGER :: ji , jj , jk ! dummy loop indices 642 INTEGER :: iku, ikv, ik, ikm1 ! local integers 631 643 REAL(wp) :: zeps, zm1_g, zm1_2g ! local scalars 632 644 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - … … 644 656 wslpjml(1,:) = 0._wp ; wslpjml(jpi,:) = 0._wp 645 657 ! 646 ! !== surface mixed layer mask !647 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise658 ! !== surface mixed layer mask ! 659 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise 648 660 # if defined key_vectopt_loop 649 661 DO jj = 1, 1 650 DO ji = 1, jpij ! vector opt. (forced unrolling)662 DO ji = 1, jpij ! vector opt. (forced unrolling) 651 663 # else 652 664 DO jj = 1, jpj … … 679 691 DO ji = 2, jpim1 680 692 # endif 681 ! !== Slope at u- & v-points just below the Mixed Layer ==!693 ! !== Slope at u- & v-points just below the Mixed Layer ==! 682 694 ! 683 ! 695 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 684 696 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 685 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 697 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 686 698 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 687 699 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 688 ! 700 ! !- horizontal density gradient at u- & v-points 689 701 zau = p_gru(ji,jj,iku) / e1u(ji,jj) 690 702 zav = p_grv(ji,jj,ikv) / e2v(ji,jj) 691 ! 692 ! 703 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 704 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 693 705 zbu = MIN( zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) ) 694 706 zbv = MIN( zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) ) 695 ! 707 ! !- Slope at u- & v-points (uslpml, vslpml) 696 708 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 697 709 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 698 710 ! 699 ! !== i- & j-slopes at w-points just below the Mixed Layer ==!711 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! 700 712 ! 701 713 ik = MIN( nmln(ji,jj) + 1, jpk ) 702 714 ikm1 = MAX( 1, ik-1 ) 703 ! 715 ! !- vertical density gradient for w-slope (from N^2) 704 716 zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 705 ! 717 ! !- horizontal density i- & j-gradient at w-points 706 718 zci = MAX( umask(ji-1,jj,ik ) + umask(ji,jj,ik ) & 707 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 719 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 708 720 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 709 721 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) … … 712 724 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & 713 725 & + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik) 714 ! 715 ! 726 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 727 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 716 728 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zai ) ) 717 729 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj ) ) 718 ! 730 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 719 731 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 720 732 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 721 733 END DO 722 734 END DO 723 !!gm this lbc_lnk should be useless....735 !!gm this lbc_lnk should be useless.... 724 736 CALL lbc_lnk( uslpml , 'U', -1. ) ; CALL lbc_lnk( vslpml , 'V', -1. ) ! lateral boundary cond. (sign change) 725 737 CALL lbc_lnk( wslpiml, 'W', -1. ) ; CALL lbc_lnk( wslpjml, 'W', -1. ) ! lateral boundary conditions … … 734 746 !! ** Purpose : Initialization for the isopycnal slopes computation 735 747 !! 736 !! ** Method : read the nammbf namelist and check the parameter 737 !! 748 !! ** Method : read the nammbf namelist and check the parameter 749 !! values called by tra_dmp at the first timestep (nit000) 738 750 !!---------------------------------------------------------------------- 739 751 INTEGER :: ji, jj, jk ! dummy loop indices 740 752 INTEGER :: ierr ! local integer 741 753 !!---------------------------------------------------------------------- 742 743 IF(lwp) THEN 754 755 IF(lwp) THEN 744 756 WRITE(numout,*) 745 757 WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing' 746 758 WRITE(numout,*) '~~~~~~~~~~~~' 747 759 ENDIF 748 760 749 761 IF( ln_traldf_grif ) THEN ! Griffies operator : triad of slopes 750 762 ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , wslp2(jpi,jpj,jpk) , STAT=ierr ) … … 755 767 IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 756 768 ! 757 IF( ( ln_traldf_hor .OR. ln_dynldf_hor ) .AND. ln_sco ) &758 CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator in s-coordinate not supported' )759 !760 769 ELSE ! Madec operator : slopes at u-, v-, and w-points 761 770 ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) , & … … 770 779 wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp 771 780 772 !!gm I no longer understand this.....781 !!gm I no longer understand this..... 773 782 IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 774 783 IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' … … 803 812 LOGICAL, PUBLIC, PARAMETER :: lk_ldfslp = .FALSE. !: slopes flag 804 813 CONTAINS 805 SUBROUTINE ldf_slp( kt, prd, pn2 ) 806 INTEGER, INTENT(in) :: kt 814 SUBROUTINE ldf_slp( kt, prd, pn2 ) ! Dummy routine 815 INTEGER, INTENT(in) :: kt 807 816 REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2 808 817 WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1) … … 812 821 WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt 813 822 END SUBROUTINE ldf_slp_grif 814 SUBROUTINE ldf_slp_init ! Dummy routine823 SUBROUTINE ldf_slp_init ! Dummy routine 815 824 END SUBROUTINE ldf_slp_init 816 825 #endif -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
r2715 r3101 67 67 & ln_traldf_level, ln_traldf_hor , ln_traldf_iso, & 68 68 & ln_traldf_grif , ln_traldf_gdia, & 69 & ln_triad_iso , ln_botmix_grif, & 69 70 & rn_aht_0 , rn_ahtb_0 , rn_aeiv_0, & 70 71 & rn_slpmax … … 94 95 WRITE(numout,*) ' maximum isoppycnal slope rn_slpmax = ', rn_slpmax 95 96 WRITE(numout,*) ' + griffies operator internal controls not set via the namelist (experimental): ' 96 WRITE(numout,*) ' calculate triads twice l _triad_iso = ', l_triad_iso97 WRITE(numout,*) ' no Shapiro filter l_no_smooth = ', l_no_smooth97 WRITE(numout,*) ' calculate triads twice ln_triad_iso = ', ln_triad_iso 98 WRITE(numout,*) ' GM -->lat mixing on bottom ln_botmix_grif = ', ln_botmix_grif 98 99 WRITE(numout,*) 99 100 ENDIF -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90
r2715 r3101 32 32 33 33 REAL(wp), PUBLIC :: aht0, ahtb0, aeiv0 !!: OLD namelist names 34 LOGICAL , PUBLIC :: l_triad_iso = .FALSE. !: calculate triads twice 35 LOGICAL , PUBLIC :: l_no_smooth = .FALSE. !: no Shapiro smoothing 34 LOGICAL , PUBLIC :: ln_triad_iso = .FALSE. !: calculate triads twice 35 LOGICAL , PUBLIC :: ln_botmix_grif = .FALSE. !: mixing on bottom 36 LOGICAL , PUBLIC :: l_grad_zps = .FALSE. !: special treatment for Horz Tgradients w partial steps 36 37 37 38 #if defined key_traldf_c3d -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r2715 r3101 3 3 !! *** MODULE eosbn2 *** 4 4 !! Ocean diagnostic variable : equation of state - in situ and potential density 5 !! - Brunt-Vaisala frequency 5 !! - Brunt-Vaisala frequency 6 6 !!============================================================================== 7 7 !! History : OPA ! 1989-03 (O. Marti) Original code … … 27 27 !! eos_insitu_2d : Compute the in situ density for 2d fields 28 28 !! eos_bn2 : Compute the Brunt-Vaisala frequency 29 !! eos_alpbet : calculates the in situ thermal and haline expansion coeff.29 !! eos_alpbet : calculates the in situ thermal/haline expansion ratio 30 30 !! tfreez : Compute the surface freezing temperature 31 31 !! eos_init : set eos parameters (namelist) … … 41 41 PRIVATE 42 42 43 ! !! * Interface 43 ! !! * Interface 44 44 INTERFACE eos 45 45 MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d 46 END INTERFACE 46 END INTERFACE 47 47 INTERFACE bn2 48 48 MODULE PROCEDURE eos_bn2 49 END INTERFACE 49 END INTERFACE 50 50 51 51 PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules … … 61 61 62 62 REAL(wp), PUBLIC :: ralpbet !: alpha / beta ratio 63 63 64 64 !! * Substitutions 65 65 # include "domzgr_substitute.h90" … … 75 75 !!---------------------------------------------------------------------- 76 76 !! *** ROUTINE eos_insitu *** 77 !! 78 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 77 !! 78 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 79 79 !! potential temperature and salinity using an equation of state 80 80 !! defined through the namelist parameter nn_eos. … … 134 134 !CDIR NOVERRCHK 135 135 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 136 ! 136 ! 137 137 DO jk = 1, jpkm1 138 138 DO jj = 1, jpj … … 199 199 !!---------------------------------------------------------------------- 200 200 !! *** ROUTINE eos_insitu_pot *** 201 !! 201 !! 202 202 !! ** Purpose : Compute the in situ density (ratio rho/rau0) and the 203 203 !! potential volumic mass (Kg/m3) from potential temperature and 204 !! salinity fields using an equation of state defined through the 204 !! salinity fields using an equation of state defined through the 205 205 !! namelist parameter nn_eos. 206 206 !! … … 230 230 !! nn_eos = 2 : linear equation of state function of temperature and 231 231 !! salinity 232 !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0 232 !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0 233 233 !! = rn_beta * s - rn_alpha * tn - 1. 234 234 !! rhop(t,s) = rho(t,s) … … 265 265 !CDIR NOVERRCHK 266 266 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 267 ! 267 ! 268 268 DO jk = 1, jpkm1 269 269 DO jj = 1, jpj … … 336 336 !! *** ROUTINE eos_insitu_2d *** 337 337 !! 338 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 338 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 339 339 !! potential temperature and salinity using an equation of state 340 340 !! defined through the namelist parameter nn_eos. * 2D field case … … 374 374 ! ! 2 : salinity [psu] 375 375 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m] 376 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 376 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 377 377 !! 378 378 INTEGER :: ji, jj ! dummy loop indices … … 449 449 DO jj = 1, jpjm1 450 450 DO ji = 1, fs_jpim1 ! vector opt. 451 prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 451 prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 452 452 END DO 453 453 END DO … … 468 468 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the time- 469 469 !! step of the input arguments 470 !! 470 !! 471 471 !! ** Method : 472 472 !! * nn_eos = 0 : UNESCO sea water properties … … 482 482 !! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 483 483 !! The use of potential density to compute N^2 introduces e r r o r 484 !! in the sign of N^2 at great depths. We recommand the use of 484 !! in the sign of N^2 at great depths. We recommand the use of 485 485 !! nn_eos = 0, except for academical studies. 486 486 !! Macro-tasked on horizontal slab (jk-loop) … … 497 497 !! 498 498 INTEGER :: ji, jj, jk ! dummy loop indices 499 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars 499 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars 500 500 #if defined key_zdfddm 501 501 REAL(wp) :: zds ! local scalars … … 504 504 505 505 ! pn2 : interior points only (2=< jk =< jpkm1 ) 506 ! -------------------------- 506 ! -------------------------- 507 507 ! 508 508 SELECT CASE( nn_eos ) … … 542 542 & - 0.121555e-07_wp ) * zh 543 543 ! 544 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 544 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 545 545 & * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 546 546 & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) … … 565 565 & - rn_beta * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) ) ) & 566 566 & / fse3w(:,:,jk) * tmask(:,:,jk) 567 END DO 567 END DO 568 568 #if defined key_zdfddm 569 569 DO jk = 2, jpkm1 ! Rrau = (alpha / beta) (dk[t] / dk[s]) 570 570 DO jj = 1, jpj 571 571 DO ji = 1, jpi 572 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 572 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 573 573 IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 574 574 rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds … … 587 587 588 588 589 SUBROUTINE eos_alpbet( pts, palp h, pbeta)590 !!---------------------------------------------------------------------- 591 !! *** ROUTINE ldf_slp_grif***592 !! 593 !! ** Purpose : Calculates the thermal and haline expansion coefficients at T-points.594 !! 595 !! ** Method : calculates alpha and beta at T-points589 SUBROUTINE eos_alpbet( pts, palpbet, beta0 ) 590 !!---------------------------------------------------------------------- 591 !! *** ROUTINE eos_alpbet *** 592 !! 593 !! ** Purpose : Calculates the in situ thermal/haline expansion ratio at T-points 594 !! 595 !! ** Method : calculates alpha / beta ratio at T-points 596 596 !! * nn_eos = 0 : UNESCO sea water properties 597 !! The brunt-vaisala frequency is computed using the polynomial 598 !! polynomial expression of McDougall (1987): 599 !! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 600 !! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 601 !! computed and used in zdfddm module : 602 !! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 597 !! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial 598 !! polynomial expression of McDougall (1987). 599 !! Scalar beta0 is returned = 1. 603 600 !! * nn_eos = 1 : linear equation of state (temperature only) 604 !! N^2 = grav * rn_alpha * dk[ t ]/e3w 601 !! The ratio is undefined, so we return alpha as palpbet 602 !! Scalar beta0 is returned = 0. 605 603 !! * nn_eos = 2 : linear equation of state (temperature & salinity) 606 !! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 607 !! * nn_eos = 3 : Jackett JAOT 2003 ??? 608 !! 609 !! ** Action : - palph, pbeta : thermal and haline expansion coeff. at T-point 610 !!---------------------------------------------------------------------- 611 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 612 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palph, pbeta ! thermal & haline expansion coeff. 613 ! 604 !! The alpha/beta ratio is returned as ralpbet 605 !! Scalar beta0 is returned = 1. 606 !! 607 !! ** Action : - palpbet : thermal/haline expansion ratio at T-points 608 !! : beta0 : 1. or 0. 609 !!---------------------------------------------------------------------- 610 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 611 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palpbet ! thermal/haline expansion ratio 612 REAL(wp), INTENT( out) :: beta0 ! set = 1 except with case 1 eos, rho=rho(T) 613 !! 614 614 INTEGER :: ji, jj, jk ! dummy loop indices 615 REAL(wp) :: zt, zs, zh ! local scalars 615 REAL(wp) :: zt, zs, zh ! local scalars 616 616 !!---------------------------------------------------------------------- 617 617 ! … … 624 624 zt = pts(ji,jj,jk,jp_tem) ! potential temperature 625 625 zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35) 626 zh = fsdept(ji,jj,jk) ! depth in meters 627 ! 628 pbeta(ji,jj,jk) = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt & 629 & - 0.301985e-05_wp ) * zt & 630 & + 0.785567e-03_wp & 631 & + ( 0.515032e-08_wp * zs & 632 & + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs & 633 & + ( ( 0.121551e-17_wp * zh & 634 & - 0.602281e-15_wp * zs & 635 & - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh & 636 & + 0.408195e-10_wp * zs & 637 & + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt & 638 & - 0.121555e-07_wp ) * zh 639 ! 640 palph(ji,jj,jk) = - pbeta(ji,jj,jk) * & 641 & ((( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & 642 & - 0.203814e-03_wp ) * zt & 643 & + 0.170907e-01_wp ) * zt & 644 & + 0.665157e-01_wp & 645 & + ( - 0.678662e-05_wp * zs & 646 & - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs & 647 & + ( ( - 0.302285e-13_wp * zh & 648 & - 0.251520e-11_wp * zs & 649 & + 0.512857e-12_wp * zt * zt ) * zh & 650 & - 0.164759e-06_wp * zs & 651 & +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt & 652 & + 0.380374e-04_wp ) * zh) 626 zh = fsdept(ji,jj,jk) ! depth in meters 627 ! 628 palpbet(ji,jj,jk) = & 629 & ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & 630 & - 0.203814e-03_wp ) * zt & 631 & + 0.170907e-01_wp ) * zt & 632 & + 0.665157e-01_wp & 633 & + ( - 0.678662e-05_wp * zs & 634 & - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs & 635 & + ( ( - 0.302285e-13_wp * zh & 636 & - 0.251520e-11_wp * zs & 637 & + 0.512857e-12_wp * zt * zt ) * zh & 638 & - 0.164759e-06_wp * zs & 639 & +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt & 640 & + 0.380374e-04_wp ) * zh 653 641 END DO 654 642 END DO 655 643 END DO 656 ! 657 CASE ( 1 ) 658 palph(:,:,:) = - rn_alpha 659 pbeta(:,:,:) = 0._wp 660 ! 661 CASE ( 2 ) 662 palph(:,:,:) = - rn_alpha 663 pbeta(:,:,:) = rn_beta 644 beta0 = 1._wp 645 ! 646 CASE ( 1 ) !== Linear formulation = F( temperature ) ==! 647 palpbet(:,:,:) = rn_alpha 648 beta0 = 0._wp 649 ! 650 CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 651 palpbet(:,:,:) = ralpbet 652 beta0 = 1._wp 664 653 ! 665 654 CASE DEFAULT … … 747 736 748 737 !!====================================================================== 749 END MODULE eosbn2 738 END MODULE eosbn2 -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90
r3094 r3101 4 4 !! Ocean tracers: horizontal component of the lateral tracer mixing trend 5 5 !!====================================================================== 6 !! History : 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) 6 !! History : 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) 7 7 !! ! Griffies operator version adapted from traldf_iso.F90 8 8 !!---------------------------------------------------------------------- … … 11 11 !! 'key_ldfslp' slope of the lateral diffusive direction 12 12 !!---------------------------------------------------------------------- 13 !! tra_ldf_iso_grif : update the tracer trend with the horizontal component 14 !! of the Griffies iso-neutral laplacian operator 13 !! tra_ldf_iso_grif : update the tracer trend with the horizontal component 14 !! of the Griffies iso-neutral laplacian operator 15 15 !!---------------------------------------------------------------------- 16 16 USE oce ! ocean dynamics and active tracers … … 34 34 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: psix_eiv, psiy_eiv !: eiv stream function (diag only) 35 35 REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: ah_wslp2 !: aeiv*w-slope^2 36 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zdkt ! atypic workspace36 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zdkt3d !: vertical tracer gradient at 2 levels 37 37 38 38 !! * Substitutions … … 53 53 !! *** ROUTINE tra_ldf_iso_grif *** 54 54 !! 55 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 56 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 55 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 56 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 57 57 !! add it to the general trend of tracer equation. 58 58 !! 59 !! ** Method : The horizontal component of the lateral diffusive trends 59 !! ** Method : The horizontal component of the lateral diffusive trends 60 60 !! is provided by a 2nd order operator rotated along neural or geopo- 61 61 !! tential surfaces to which an eddy induced advection can be added … … 67 67 !! 68 68 !! 2nd part : horizontal fluxes of the lateral mixing operator 69 !! ======== 69 !! ======== 70 70 !! zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 71 71 !! - aht e2u*uslp dk[ mi(mk(tb)) ] … … 95 95 ! 96 96 INTEGER , INTENT(in ) :: kt ! ocean time-step index 97 INTEGER , INTENT(in ) :: kit000 97 INTEGER , INTENT(in ) :: kit000 ! first time step index 98 98 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 99 99 INTEGER , INTENT(in ) :: kjpt ! number of tracers 100 100 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 101 101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 102 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 102 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 103 103 REAL(wp) , INTENT(in ) :: pahtb0 ! background diffusion coef 104 104 ! … … 109 109 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 110 110 REAL(wp) :: zcoef0, zbtr ! - - 111 !REAL(wp), POINTER, DIMENSION(:,:,:) :: zdkt ! 2D+1 workspace112 111 ! 113 112 REAL(wp) :: zslope_skew, zslope_iso, zslope2, zbu, zbv … … 122 121 CALL ctl_stop('tra_ldf_iso_grif: requested workspace arrays unavailable.') ; RETURN 123 122 ENDIF 124 ! ARP - line below uses 'bounds re-mapping' which is only defined in 125 ! Fortran 2003 and up. We would be OK if code was written to use 126 ! zdkt(:,:,1:2) instead as then wouldn't need to re-map bounds. 127 ! As it is, we make zdkt a module array and allocate it in _alloc(). 128 !zdkt(1:jpi,1:jpj,0:1) => wrk_3d_9(:,:,1:2) 129 130 IF( kt == kit000 ) THEN 123 124 IF( kt == kit000 .AND. .NOT.ALLOCATED(ah_wslp2) ) THEN 131 125 IF(lwp) WRITE(numout,*) 132 126 IF(lwp) WRITE(numout,*) 'tra_ldf_iso_grif : rotated laplacian diffusion operator on ', cdtype 133 IF(lwp) WRITE(numout,*) ' WARNING: STILL UNDER TEST, NOT RECOMMENDED. USE AT YOUR OWN PERIL'134 127 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 135 IF (.not. ALLOCATED(ah_wslp2))THEN 136 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , zdkt(jpi,jpj,0:1), STAT=ierr ) 137 ENDIF 128 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , zdkt3d(jpi,jpj,0:1), STAT=ierr ) 138 129 IF( lk_mpp ) CALL mpp_sum ( ierr ) 139 130 IF( ierr > 0 ) CALL ctl_stop('STOP', 'tra_ldf_iso_grif: unable to allocate arrays') … … 148 139 149 140 !!---------------------------------------------------------------------- 150 !! 0 - calculate ah_wslp2, psix_eiv, psiy_eiv 141 !! 0 - calculate ah_wslp2, psix_eiv, psiy_eiv 151 142 !!---------------------------------------------------------------------- 152 143 153 !!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads144 !!gm Future development: consider using Ah defined at T-points and attached to the 4 t-point triads 154 145 155 146 ah_wslp2(:,:,:) = 0._wp … … 164 155 DO jj = 1, jpjm1 165 156 DO ji = 1, fs_jpim1 157 ze1ur = 1._wp / e1u(ji,jj) 166 158 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 167 159 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 168 zah = fsahtu(ji,jj,jk) ! 160 zah = fsahtu(ji,jj,jk) ! fsaht(ji+ip,jj,jk) 169 161 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 170 zslope2 = zslope_skew - ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+kp) 162 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 163 ! (do this by *adding* gradient of depth) 164 zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+kp) 171 165 zslope2 = zslope2 *zslope2 172 166 ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) & 173 167 & + zah * ( zbu * ze3wr / ( e1t(ji+ip,jj) * e2t(ji+ip,jj) ) ) * zslope2 174 168 IF( ln_traldf_gdia ) THEN 175 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew !fsaeit(ji+ip,jj,jk)*zslope_skew169 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew ! fsaeit(ji+ip,jj,jk)*zslope_skew 176 170 psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp * zaei_slp 177 171 ENDIF … … 187 181 DO jj = 1, jpjm1 188 182 DO ji=1,fs_jpim1 183 ze2vr = 1._wp / e2v(ji,jj) 189 184 ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp) 190 185 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 191 zah = fsaht u(ji,jj,jk) !fsaht(ji,jj+jp,jk)186 zah = fsahtv(ji,jj,jk) ! fsaht(ji,jj+jp,jk) 192 187 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 193 zslope2 = zslope_skew - ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+kp) 188 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 189 ! (do this by *adding* gradient of depth) 190 zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+kp) 194 191 zslope2 = zslope2 * zslope2 195 192 ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) & 196 193 & + zah * ( zbv * ze3wr / ( e1t(ji,jj+jp) * e2t(ji,jj+jp) ) ) * zslope2 197 194 IF( ln_traldf_gdia ) THEN 198 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew !fsaeit(ji,jj+jp,jk)*zslope_skew195 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew 199 196 psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp * zaei_slp 200 197 ENDIF … … 212 209 zftu(:,:,:) = 0._wp 213 210 zftv(:,:,:) = 0._wp 214 ! 211 ! 215 212 DO jk = 1, jpkm1 !== before lateral T & S gradients at T-level jk ==! 216 213 DO jj = 1, jpjm1 … … 221 218 END DO 222 219 END DO 223 IF( ln_zps ) THEN! partial steps: correction at the last level220 IF( ln_zps.and.l_grad_zps ) THEN ! partial steps: correction at the last level 224 221 # if defined key_vectopt_loop 225 222 DO jj = 1, 1 … … 229 226 DO ji = 1, jpim1 230 227 # endif 231 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 232 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 228 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 229 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 233 230 END DO 234 231 END DO … … 242 239 ! 243 240 ! !== Vertical tracer gradient at level jk and jk+1 244 zdkt (:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1)241 zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 245 242 ! 246 ! ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)247 IF( jk == 1 ) THEN ; zdkt (:,:,0) = zdkt(:,:,1)248 ELSE ; zdkt (:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk)243 ! ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1) 244 IF( jk == 1 ) THEN ; zdkt3d(:,:,0) = zdkt3d(:,:,1) 245 ELSE ; zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk) 249 246 ENDIF 250 247 251 IF( .NOT. l_triad_iso ) THEN 252 triadi = triadi_g 253 triadj = triadj_g 254 ENDIF 255 256 DO ip = 0, 1 !== Horizontal & vertical fluxes 257 DO kp = 0, 1 258 DO jj = 1, jpjm1 259 DO ji = 1, fs_jpim1 260 ze1ur = 1._wp / e1u(ji,jj) 261 zdxt = zdit(ji,jj,jk) * ze1ur 262 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 263 zdzt = zdkt(ji+ip,jj,kp) * ze3wr 264 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 265 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 266 267 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 268 zah = fsahtu(ji,jj,jk) !*umask(ji,jj,jk+kp) !fsaht(ji+ip,jj,jk) ===>> ???? 269 zah_slp = zah * zslope_iso 270 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew !fsaeit(ji+ip,jj,jk)*zslope_skew 271 zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 272 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 248 249 IF (ln_botmix_grif) THEN 250 DO ip = 0, 1 !== Horizontal & vertical fluxes 251 DO kp = 0, 1 252 DO jj = 1, jpjm1 253 DO ji = 1, fs_jpim1 254 ze1ur = 1._wp / e1u(ji,jj) 255 zdxt = zdit(ji,jj,jk) * ze1ur 256 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 257 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 258 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 259 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 260 261 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 262 ! ln_botmix_grif is .T. don't mask zah for bottom half cells 263 zah = fsahtu(ji,jj,jk) !*umask(ji,jj,jk+kp) !fsaht(ji+ip,jj,jk) ===>> ???? 264 zah_slp = zah * zslope_iso 265 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew !fsaeit(ji+ip,jj,jk)*zslope_skew 266 zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 267 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 268 END DO 273 269 END DO 274 270 END DO 275 271 END DO 276 END DO 277 278 DO jp = 0, 1 279 DO kp = 0, 1 280 DO jj = 1, jpjm1 281 DO ji = 1, fs_jpim1 282 ze2vr = 1._wp / e2v(ji,jj) 283 zdyt = zdjt(ji,jj,jk) * ze2vr 284 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 285 zdzt = zdkt(ji,jj+jp,kp) * ze3wr 286 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 287 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 288 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 289 zah = fsahtv(ji,jj,jk) !*vmask(ji,jj,jk+kp) !fsaht(ji,jj+jp,jk) 290 zah_slp = zah * zslope_iso 291 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew !fsaeit(ji,jj+jp,jk)*zslope_skew 292 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 293 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 272 273 DO jp = 0, 1 274 DO kp = 0, 1 275 DO jj = 1, jpjm1 276 DO ji = 1, fs_jpim1 277 ze2vr = 1._wp / e2v(ji,jj) 278 zdyt = zdjt(ji,jj,jk) * ze2vr 279 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 280 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 281 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 282 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 283 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 284 ! ln_botmix_grif is .T. don't mask zah for bottom half cells 285 zah = fsahtv(ji,jj,jk) !*vmask(ji,jj,jk+kp) ! fsaht(ji,jj+jp,jk) 286 zah_slp = zah * zslope_iso 287 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew 288 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 289 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 290 END DO 294 291 END DO 295 292 END DO 296 293 END DO 297 END DO 298 299 ! !== divergence and add to the general trend ==! 294 ELSE 295 DO ip = 0, 1 !== Horizontal & vertical fluxes 296 DO kp = 0, 1 297 DO jj = 1, jpjm1 298 DO ji = 1, fs_jpim1 299 ze1ur = 1._wp / e1u(ji,jj) 300 zdxt = zdit(ji,jj,jk) * ze1ur 301 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 302 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 303 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 304 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 305 306 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 307 ! ln_botmix_grif is .F. mask zah for bottom half cells 308 zah = fsahtu(ji,jj,jk) * umask(ji,jj,jk+kp) ! fsaht(ji+ip,jj,jk) ===>> ???? 309 zah_slp = zah * zslope_iso 310 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew ! fsaeit(ji+ip,jj,jk)*zslope_skew 311 zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 312 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 313 END DO 314 END DO 315 END DO 316 END DO 317 318 DO jp = 0, 1 319 DO kp = 0, 1 320 DO jj = 1, jpjm1 321 DO ji = 1, fs_jpim1 322 ze2vr = 1._wp / e2v(ji,jj) 323 zdyt = zdjt(ji,jj,jk) * ze2vr 324 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 325 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 326 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 327 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 328 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 329 ! ln_botmix_grif is .F. mask zah for bottom half cells 330 zah = fsahtv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! fsaht(ji,jj+jp,jk) 331 zah_slp = zah * zslope_iso 332 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew 333 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 334 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 335 END DO 336 END DO 337 END DO 338 END DO 339 END IF 340 ! !== divergence and add to the general trend ==! 300 341 DO jj = 2 , jpjm1 301 DO ji = fs_2, fs_jpim1 342 DO ji = fs_2, fs_jpim1 ! vector opt. 302 343 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 303 344 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zbtr * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) & … … 308 349 END DO 309 350 ! 310 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to the general tracer trend351 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to the general tracer trend 311 352 DO jj = 2, jpjm1 312 DO ji = fs_2, fs_jpim1 353 DO ji = fs_2, fs_jpim1 ! vector opt. 313 354 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) & 314 355 & / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) … … 317 358 END DO 318 359 ! 319 ! ! "Poleward" diffusive heat or salt transports (T-S case only)360 ! ! "Poleward" diffusive heat or salt transports (T-S case only) 320 361 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 321 362 IF( jn == jp_tem) htr_ldf(:) = ptr_vj( zftv(:,:,:) ) ! 3.3 names … … 325 366 #if defined key_diaar5 326 367 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 327 z2d(:,:) = 0._wp 328 zztmp = rau0 * rcp 368 z2d(:,:) = 0._wp 369 zztmp = rau0 * rcp 329 370 DO jk = 1, jpkm1 330 371 DO jj = 2, jpjm1 331 372 DO ji = fs_2, fs_jpim1 ! vector opt. 332 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 373 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 333 374 END DO 334 375 END DO … … 337 378 CALL lbc_lnk( z2d, 'U', -1. ) 338 379 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 339 z2d(:,:) = 0._wp 380 z2d(:,:) = 0._wp 340 381 DO jk = 1, jpkm1 341 382 DO jj = 2, jpjm1 342 383 DO ji = fs_2, fs_jpim1 ! vector opt. 343 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 384 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 344 385 END DO 345 386 END DO … … 347 388 z2d(:,:) = zztmp * z2d(:,:) 348 389 CALL lbc_lnk( z2d, 'V', -1. ) 349 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in i-direction390 CALL iom_put( "vdiff_heattr", z2d ) ! heat transport in j-direction 350 391 END IF 351 392 #endif -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap.F90
r3094 r3101 83 83 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 84 84 ! 85 ALLOCATE( e1ur(jpi,jpj), e2vr(jpi,jpj), STAT=ierr ) 86 IF( lk_mpp ) CALL mpp_sum( ierr ) 87 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'tra_ldf_lap : unable to allocate arrays' ) 88 ! 89 e1ur(:,:) = e2u(:,:) / e1u(:,:) 90 e2vr(:,:) = e1v(:,:) / e2v(:,:) 85 IF( .NOT. ALLOCATED( e1ur ) ) THEN 86 ! This routine may be called for both active and passive tracers. 87 ! Allocate and set saved arrays on first call only. 88 ALLOCATE( e1ur(jpi,jpj), e2vr(jpi,jpj), STAT=ierr ) 89 IF( lk_mpp ) CALL mpp_sum( ierr ) 90 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'tra_ldf_lap : unable to allocate arrays' ) 91 ! 92 e1ur(:,:) = e2u(:,:) / e1u(:,:) 93 e2vr(:,:) = e1v(:,:) / e2v(:,:) 94 ENDIF 91 95 ENDIF 92 96 -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
r2715 r3101 36 36 REAL(wp) :: rn_bfrien = 30._wp ! local factor to enhance coefficient bfri 37 37 LOGICAL :: ln_bfr2d = .false. ! logical switch for 2D enhancement 38 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bfrcoef2d ! 2D bottom drag coefficient38 LOGICAL , PUBLIC :: ln_bfrimp = .false. ! logical switch for implicit bottom friction 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bfrcoef2d ! 2D bottom drag coefficient 40 40 41 41 !! * Substitutions … … 142 142 REAL(wp) :: zfru, zfrv ! - - 143 143 !! 144 NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, ln_bfr2d, rn_bfrien 144 NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, ln_bfr2d, rn_bfrien, ln_bfrimp 145 145 !!---------------------------------------------------------------------- 146 146 … … 156 156 ! ! allocate zdfbfr arrays 157 157 IF( zdf_bfr_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_bfr_init : unable to allocate arrays' ) 158 159 ! ! Make sure ln_zdfexp=.false. when use implicit bfr 160 IF( ln_bfrimp .AND. ln_zdfexp ) THEN 161 IF(lwp) THEN 162 WRITE(numout,*) 163 WRITE(numout,*) 'Implicit bottom friction can only be used when ln_zdfexp=.false.' 164 WRITE(numout,*) ' but you set: ln_bfrimp=.true. and ln_zdfexp=.true.!!!!' 165 WRITE(ctmp1,*) ' bad ln_bfrimp value = .true.' 166 CALL ctl_stop( ctmp1 ) 167 END IF 168 END IF 158 169 159 170 SELECT CASE (nn_bfr) … … 207 218 ! 208 219 END SELECT 220 IF(lwp) WRITE(numout,*) ' implicit bottom friction switch ln_bfrimp = ', ln_bfrimp 209 221 ! 210 222 ! Basic stability check on bottom friction coefficient -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r3094 r3101 69 69 USE c1d ! 1D configuration 70 70 USE step_c1d ! Time stepping loop for the 1D configuration 71 USE dynnept ! simplified form of Neptune effect 71 72 #if defined key_top 72 73 USE trcini ! passive tracer initialisation … … 299 300 CALL dom_init ! Domain 300 301 302 IF( ln_nnogather ) CALL nemo_northcomms ! Initialise the northfold neighbour lists (must be done after the masks are defined) 303 301 304 IF( ln_ctl ) CALL prt_ctl_init ! Print control 302 305 … … 305 308 IF( lk_bdy ) CALL bdy_dta_init ! Open boundaries initialisation of external data arrays 306 309 IF( lk_bdy ) CALL tide_init ! Open boundaries initialisation of tidal harmonic forcing 310 311 CALL flush(numout) 312 CALL dyn_nept_init ! simplified form of Neptune effect 313 CALL flush(numout) 307 314 308 315 CALL istate_init ! ocean initial state (Dynamics and tracers) … … 626 633 END SUBROUTINE factorise 627 634 635 #if defined key_mpp_mpi 636 SUBROUTINE nemo_northcomms 637 !!====================================================================== 638 !! *** ROUTINE nemo_northcomms *** 639 !! nemo_northcomms : Setup for north fold exchanges with explicit peer to peer messaging 640 !!===================================================================== 641 !!---------------------------------------------------------------------- 642 !! 643 !! ** Purpose : Initialization of the northern neighbours lists. 644 !!---------------------------------------------------------------------- 645 !! 1.0 ! 2011-10 (A. C. Coward, NOCS & J. Donners, PRACE) 646 !!---------------------------------------------------------------------- 647 648 INTEGER :: ji, jj, jk, ij, jtyp ! dummy loop indices 649 INTEGER :: ijpj ! number of rows involved in north-fold exchange 650 INTEGER :: northcomms_alloc ! allocate return status 651 REAL(wp), ALLOCATABLE, DIMENSION ( :,: ) :: znnbrs ! workspace 652 LOGICAL, ALLOCATABLE, DIMENSION ( : ) :: lrankset ! workspace 653 654 IF(lwp) WRITE(numout,*) 655 IF(lwp) WRITE(numout,*) 'nemo_northcomms : Initialization of the northern neighbours lists' 656 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 657 658 !!---------------------------------------------------------------------- 659 ALLOCATE( znnbrs(jpi,jpj), stat = northcomms_alloc ) 660 ALLOCATE( lrankset(jpnij), stat = northcomms_alloc ) 661 IF( northcomms_alloc /= 0 ) THEN 662 WRITE(numout,cform_war) 663 WRITE(numout,*) 'northcomms_alloc : failed to allocate arrays' 664 CALL ctl_stop( 'STOP', 'nemo_northcomms : unable to allocate temporary arrays' ) 665 ENDIF 666 nsndto = 0 667 isendto = -1 668 ijpj = 4 669 ! 670 ! This routine has been called because ln_nnogather has been set true ( nammpp ) 671 ! However, these first few exchanges have to use the mpi_allgather method to 672 ! establish the neighbour lists to use in subsequent peer to peer exchanges. 673 ! Consequently, set l_north_nogather to be false here and set it true only after 674 ! the lists have been established. 675 ! 676 l_north_nogather = .FALSE. 677 ! 678 ! Exchange and store ranks on northern rows 679 680 DO jtyp = 1,4 681 682 lrankset = .FALSE. 683 znnbrs = narea 684 SELECT CASE (jtyp) 685 CASE(1) 686 CALL lbc_lnk( znnbrs, 'T', 1. ) ! Type 1: T,W-points 687 CASE(2) 688 CALL lbc_lnk( znnbrs, 'U', 1. ) ! Type 2: U-point 689 CASE(3) 690 CALL lbc_lnk( znnbrs, 'V', 1. ) ! Type 3: V-point 691 CASE(4) 692 CALL lbc_lnk( znnbrs, 'F', 1. ) ! Type 4: F-point 693 END SELECT 694 695 IF ( njmppt(narea) .EQ. MAXVAL( njmppt ) ) THEN 696 DO jj = nlcj-ijpj+1, nlcj 697 ij = jj - nlcj + ijpj 698 DO ji = 1,jpi 699 IF ( INT(znnbrs(ji,jj)) .NE. 0 .AND. INT(znnbrs(ji,jj)) .NE. narea ) & 700 & lrankset(INT(znnbrs(ji,jj))) = .true. 701 END DO 702 END DO 703 704 DO jj = 1,jpnij 705 IF ( lrankset(jj) ) THEN 706 nsndto(jtyp) = nsndto(jtyp) + 1 707 IF ( nsndto(jtyp) .GT. jpmaxngh ) THEN 708 CALL ctl_stop( ' Too many neighbours in nemo_northcomms ', & 709 & ' jpmaxngh will need to be increased ') 710 ENDIF 711 isendto(nsndto(jtyp),jtyp) = jj-1 ! narea converted to MPI rank 712 ENDIF 713 END DO 714 ENDIF 715 716 END DO 717 718 ! 719 ! Type 5: I-point 720 ! 721 ! ICE point exchanges may involve some averaging. The neighbours list is 722 ! built up using two exchanges to ensure that the whole stencil is covered. 723 ! lrankset should not be reset between these 'J' and 'K' point exchanges 724 725 jtyp = 5 726 lrankset = .FALSE. 727 znnbrs = narea 728 CALL lbc_lnk( znnbrs, 'J', 1. ) ! first ice U-V point 729 730 IF ( njmppt(narea) .EQ. MAXVAL( njmppt ) ) THEN 731 DO jj = nlcj-ijpj+1, nlcj 732 ij = jj - nlcj + ijpj 733 DO ji = 1,jpi 734 IF ( INT(znnbrs(ji,jj)) .NE. 0 .AND. INT(znnbrs(ji,jj)) .NE. narea ) & 735 & lrankset(INT(znnbrs(ji,jj))) = .true. 736 END DO 737 END DO 738 ENDIF 739 740 znnbrs = narea 741 CALL lbc_lnk( znnbrs, 'K', 1. ) ! second ice U-V point 742 743 IF ( njmppt(narea) .EQ. MAXVAL( njmppt )) THEN 744 DO jj = nlcj-ijpj+1, nlcj 745 ij = jj - nlcj + ijpj 746 DO ji = 1,jpi 747 IF ( INT(znnbrs(ji,jj)) .NE. 0 .AND. INT(znnbrs(ji,jj)) .NE. narea ) & 748 & lrankset( INT(znnbrs(ji,jj))) = .true. 749 END DO 750 END DO 751 752 DO jj = 1,jpnij 753 IF ( lrankset(jj) ) THEN 754 nsndto(jtyp) = nsndto(jtyp) + 1 755 IF ( nsndto(jtyp) .GT. jpmaxngh ) THEN 756 CALL ctl_stop( ' Too many neighbours in nemo_northcomms ', & 757 & ' jpmaxngh will need to be increased ') 758 ENDIF 759 isendto(nsndto(jtyp),jtyp) = jj-1 ! narea converted to MPI rank 760 ENDIF 761 END DO 762 ! 763 ! For northern row areas, set l_north_nogather so that all subsequent exchanges 764 ! can use peer to peer communications at the north fold 765 ! 766 l_north_nogather = .TRUE. 767 ! 768 ENDIF 769 DEALLOCATE( znnbrs ) 770 DEALLOCATE( lrankset ) 771 772 END SUBROUTINE nemo_northcomms 773 #else 774 SUBROUTINE nemo_northcomms ! Dummy routine 775 WRITE(*,*) 'nemo_northcomms: You should not have seen this print! error?' 776 END SUBROUTINE nemo_northcomms 777 #endif 628 778 !!====================================================================== 629 779 END MODULE nemogcm -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/step.F90
r3094 r3101 36 36 #endif 37 37 USE asminc ! assimilation increments (tra_asm_inc, dyn_asm_inc routines) 38 USE dynnept ! simplified form of Neptune effect 38 39 39 40 IMPLICIT NONE … … 220 221 IF( ln_asmiau .AND. & 221 222 & ln_dyninc ) CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment 223 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! subtract Neptune velocities (simplified) 222 224 CALL dyn_adv( kstp ) ! advection (vector or flux form) 223 225 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 224 226 CALL dyn_ldf( kstp ) ! lateral mixing 227 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! add Neptune velocities (simplified) 225 228 #if defined key_agrif 226 229 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momemtum sponge -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90
r3094 r3101 102 102 zwn(:,:,jpk) = 0.e0 ! no transport trough the bottom 103 103 104 !! add the eiv transport (if necessary)105 IF( lk_traldf_eiv )CALL tra_adv_eiv( kt, nittrc000, zun, zvn, zwn, 'TRC' )104 IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif ) & ! add the eiv transport (if necessary) 105 & CALL tra_adv_eiv( kt, nittrc000, zun, zvn, zwn, 'TRC' ) 106 106 ! 107 107 SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==! -
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90
r3094 r3101 2 2 !!====================================================================== 3 3 !! *** MODULE trcldf *** 4 !! Ocean Passive tracers : lateral diffusive trends 4 !! Ocean Passive tracers : lateral diffusive trends 5 5 !!===================================================================== 6 6 !! History : 9.0 ! 2005-11 (G. Madec) Original code 7 !! NEMO 3.0 ! 2008-01 (C. Ethe, G. Madec) merge TRC-TRA 7 !! NEMO 3.0 ! 2008-01 (C. Ethe, G. Madec) merge TRC-TRA 8 8 !!---------------------------------------------------------------------- 9 9 #if defined key_top … … 23 23 USE traldf_bilap ! lateral mixing (tra_ldf_bilap routine) 24 24 USE traldf_iso ! lateral mixing (tra_ldf_iso routine) 25 USE traldf_iso_grif ! lateral mixing (tra_ldf_iso_grif routine) 25 26 USE traldf_lap ! lateral mixing (tra_ldf_lap routine) 26 27 USE trdmod_oce … … 31 32 PRIVATE 32 33 33 PUBLIC trc_ldf ! called by step.F90 34 PUBLIC trc_ldf ! called by step.F90 34 35 ! !!: ** lateral mixing namelist (nam_trcldf) ** 35 36 INTEGER :: nldf = 0 ! type of lateral diffusion used defined from ln_trcldf_... namlist logicals) … … 39 40 !!---------------------------------------------------------------------- 40 41 !! NEMO/TOP 3.3 , NEMO Consortium (2010) 41 !! $Id$ 42 !! $Id$ 42 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 43 44 !!---------------------------------------------------------------------- … … 48 49 !!---------------------------------------------------------------------- 49 50 !! *** ROUTINE tra_ldf *** 50 !! 51 !! 51 52 !! ** Purpose : compute the lateral ocean tracer physics. 52 53 !! … … 61 62 IF( kt == nittrc000 ) CALL ldf_ctl ! initialisation & control of options 62 63 63 IF( l_trdtrc ) THEN 64 IF( l_trdtrc ) THEN 64 65 ALLOCATE( ztrtrd(jpi,jpj,jpk,jptra) ) ! temporary save of trends 65 66 ztrtrd(:,:,:,:) = tra(:,:,:,:) … … 68 69 SELECT CASE ( nldf ) ! compute lateral mixing trend and add it to the general trend 69 70 CASE ( 0 ) ; CALL tra_ldf_lap ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra ) ! iso-level laplacian 70 CASE ( 1 ) ; CALL tra_ldf_iso ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 ) ! rotated laplacian 71 CASE ( 1 ) ! rotated laplacian 72 IF( ln_traldf_grif ) THEN 73 CALL tra_ldf_iso_grif( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 ) 74 ELSE 75 CALL tra_ldf_iso ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 ) 76 ENDIF 71 77 CASE ( 2 ) ; CALL tra_ldf_bilap ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra ) ! iso-level bilaplacian 72 78 CASE ( 3 ) ; CALL tra_ldf_bilapg( kt, nittrc000, 'TRC', trb, tra, jptra ) ! s-coord. horizontal bilaplacian … … 76 82 WRITE(charout, FMT="('ldf0 ')") ; CALL prt_ctl_trc_info(charout) 77 83 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) 78 CALL tra_ldf_iso ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 ) 84 IF( ln_traldf_grif ) THEN 85 CALL tra_ldf_iso_grif( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 ) 86 ELSE 87 CALL tra_ldf_iso ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 ) 88 ENDIF 79 89 WRITE(charout, FMT="('ldf1 ')") ; CALL prt_ctl_trc_info(charout) 80 90 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) … … 92 102 CALL trd_tra( kt, 'TRC', jn, jptra_trd_ldf, ztrtrd(:,:,:,jn) ) 93 103 END DO 94 DEALLOCATE( ztrtrd ) 104 DEALLOCATE( ztrtrd ) 95 105 ENDIF 96 106 ! ! print mean trends (used for debugging) … … 106 116 !!---------------------------------------------------------------------- 107 117 !! *** ROUTINE ldf_ctl *** 108 !! 118 !! 109 119 !! ** Purpose : Choice of the operator for the lateral tracer diffusion 110 120 !! 111 121 !! ** Method : set nldf from the namtra_ldf logicals 112 !! nldf == -2 No lateral diffusion 122 !! nldf == -2 No lateral diffusion 113 123 !! nldf == -1 ESOPA test: ALL operators are used 114 124 !! nldf == 0 laplacian operator … … 117 127 !! nldf == 3 Rotated bilaplacian 118 128 !!---------------------------------------------------------------------- 119 INTEGER :: ioptio, ierr ! temporary integers 129 INTEGER :: ioptio, ierr ! temporary integers 120 130 !!---------------------------------------------------------------------- 121 131 122 132 ! Define the lateral mixing oparator for tracers 123 133 ! =============================================== 124 134 125 135 ! ! control the input 126 136 ioptio = 0 … … 163 173 ENDIF 164 174 IF ( ln_zps ) THEN ! z-coordinate 165 IF ( ln_trcldf_level ) ierr = 1 ! iso-level not allowed 175 IF ( ln_trcldf_level ) ierr = 1 ! iso-level not allowed 166 176 IF ( ln_trcldf_hor ) nldf = 2 ! horizontal (no rotation) 167 177 IF ( ln_trcldf_iso ) ierr = 2 ! isoneutral ( rotation)
Note: See TracChangeset
for help on using the changeset viewer.