Changeset 3072 for branches/2011/dev_NOC_2011_MERGE
- Timestamp:
- 2011-11-09T18:07:32+01:00 (12 years ago)
- Location:
- branches/2011/dev_NOC_2011_MERGE
- Files:
-
- 11 edited
-
DOC/TexFiles/Chapters/Chap_ZDF.tex (modified) (6 diffs)
-
DOC/TexFiles/Namelist/nambfr (modified) (1 diff)
-
NEMOGCM/CONFIG/GYRE/EXP00/namelist (modified) (1 diff)
-
NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist (modified) (1 diff)
-
NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist (modified) (1 diff)
-
NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist (modified) (1 diff)
-
NEMOGCM/CONFIG/POMME/EXP00/namelist (modified) (1 diff)
-
NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90 (modified) (2 diffs)
-
NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (22 diffs)
-
NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 (modified) (8 diffs)
-
NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Chapters/Chap_ZDF.tex
r2541 r3072 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_2011_MERGE/DOC/TexFiles/Namelist/nambfr
r2540 r3072 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_2011_MERGE/NEMOGCM/CONFIG/GYRE/EXP00/namelist
r3009 r3072 417 417 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 418 418 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d = .true.) 419 ln_bfrimp = .false. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 419 420 / 420 421 !----------------------------------------------------------------------- -
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist
r2735 r3072 417 417 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 418 418 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d=T) 419 ln_bfrimp = .false. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 419 420 / 420 421 !----------------------------------------------------------------------- -
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
r3044 r3072 417 417 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 418 418 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d=T) 419 ln_bfrimp = .false. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 419 420 / 420 421 !----------------------------------------------------------------------- -
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist
r3009 r3072 417 417 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 418 418 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d = .true.) 419 ln_bfrimp = .false. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 419 420 / 420 421 !----------------------------------------------------------------------- -
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/POMME/EXP00/namelist
r2986 r3072 417 417 ln_bfr2d = .false. ! horizontal variation of the bottom friction coef (read a 2D mask file ) 418 418 rn_bfrien = 50. ! local multiplying factor of bfr (ln_bfr2d=T) 419 ln_bfrimp = .false. ! implicit bottom friction (requires ln_zdfexp = .false. if true) 419 420 / 420 421 !----------------------------------------------------------------------- -
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
r2715 r3072 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_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r3008 r3072 119 119 INTEGER :: ji, jj, jk, jn ! dummy loop indices 120 120 INTEGER :: icycle ! local scalar 121 REAL(wp) :: zraur, zcoef, z2dt_e, z1_2dt_b ! 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 ! - - 121 INTEGER :: ikbu, ikbv ! local scalar 122 REAL(wp) :: zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf ! local scalars 123 REAL(wp) :: z1_8, zx1, zy1 ! - - 124 REAL(wp) :: z1_4, zx2, zy2 ! - - 125 REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - - 126 REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - - 127 REAL(wp) :: ua_btm, va_btm ! - - 126 128 !!---------------------------------------------------------------------- 127 129 … … 147 149 hvr_e (:,:) = hvr (:,:) 148 150 IF( ln_dynvor_een ) THEN 149 ftne(1,:) = 0. e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0151 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 150 152 DO jj = 2, jpj 151 153 DO ji = fs_2, jpi ! vector opt. 152 ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3. 153 ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3. 154 ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3. 155 ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3. 154 ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3._wp 155 ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3._wp 156 ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3._wp 157 ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3._wp 156 158 END DO 157 159 END DO … … 160 162 ENDIF 161 163 162 ! !* Local constant initialization 163 z1_2dt_b = 1._wp / ( 2.0_wp * rdt ) ! baroclinic time step 164 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt_b = 1.0_wp / rdt ! baroclinic time step (starting with euler timestep) 165 z1_8 = 0.5 * 0.25 ! coefficient for vorticity estimates 166 z1_4 = 0.5 * 0.5 167 zraur = 1. / rau0 ! 1 / volumic mass 168 ! 169 zhdiv(:,:) = 0.e0 ! barotropic divergence 170 zu_sld = 0.e0 ; zu_asp = 0.e0 ! tides trends (lk_tide=F) 171 zv_sld = 0.e0 ; zv_asp = 0.e0 164 ! !* Local constant initialization 165 z1_2dt_b = 1._wp / ( 2.0_wp * rdt ) ! reciprocal of baroclinic time step 166 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt_b = 1.0_wp / rdt ! reciprocal of baroclinic 167 ! time step (euler timestep) 168 z1_8 = 0.125_wp ! coefficient for vorticity estimates 169 z1_4 = 0.25_wp 170 zraur = 1._wp / rau0 ! 1 / volumic mass 171 ! 172 zhdiv(:,:) = 0._wp ! barotropic divergence 173 zu_sld = 0._wp ; zu_asp = 0._wp ! tides trends (lk_tide=F) 174 zv_sld = 0._wp ; zv_asp = 0._wp 175 176 IF( kt == nit000 .AND. neuler == 0) THEN ! for implicit bottom friction 177 z2dt_bf = rdt 178 ELSE 179 z2dt_bf = 2.0_wp * rdt 180 ENDIF 172 181 173 182 ! ----------------------------------------------------------------------------- … … 177 186 ! !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 178 187 ! ! -------------------------- 179 zua(:,:) = 0. e0 ; zun(:,:) = 0.e0 ; ub_b(:,:) = 0.e0180 zva(:,:) = 0. e0 ; zvn(:,:) = 0.e0 ; vb_b(:,:) = 0.e0188 zua(:,:) = 0._wp ; zun(:,:) = 0._wp ; ub_b(:,:) = 0._wp 189 zva(:,:) = 0._wp ; zvn(:,:) = 0._wp ; vb_b(:,:) = 0._wp 181 190 ! 182 191 DO jk = 1, jpkm1 … … 196 205 ! 197 206 #if defined key_vvl 198 ! 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)199 ! 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)200 207 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk) *umask(ji,jj,jk) 201 208 vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk) *vmask(ji,jj,jk) … … 273 280 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 274 281 DO ji = fs_2, fs_jpim1 275 zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)276 zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)277 END DO282 zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 283 zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) 284 END DO 278 285 END DO 279 286 … … 282 289 ! ! from the barotropic transport trend 283 290 zcoef = -1._wp * z1_2dt_b 291 292 IF(ln_bfrimp) THEN 293 ! ! Remove the bottom stress trend from 3-D sea surface level gradient 294 ! ! and Coriolis forcing in case of 3D semi-implicit bottom friction 295 DO jj = 2, jpjm1 296 DO ji = fs_2, fs_jpim1 297 ikbu = mbku(ji,jj) 298 ikbv = mbkv(ji,jj) 299 ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu) 300 va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv) 301 302 zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 303 zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 304 END DO 305 END DO 306 307 ELSE 308 284 309 # if defined key_vectopt_loop 285 DO jj = 1, 1286 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)310 DO jj = 1, 1 311 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 287 312 # else 288 DO jj = 2, jpjm1289 DO ji = 2, jpim1313 DO jj = 2, jpjm1 314 DO ji = 2, jpim1 290 315 # endif 291 316 ! Apply stability criteria for bottom friction 292 317 !RBbug for vvl and external mode we may need to use varying fse3 293 318 !!gm Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. 294 zbfru(ji,jj) = MAX( bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef ) 295 zbfrv(ji,jj) = MAX( bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef ) 296 END DO 297 END DO 298 299 IF( lk_vvl ) THEN 300 DO jj = 2, jpjm1 301 DO ji = fs_2, fs_jpim1 ! vector opt. 302 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) & 303 & / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) ) 304 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) & 305 & / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) ) 306 END DO 307 END DO 308 ELSE 309 DO jj = 2, jpjm1 310 DO ji = fs_2, fs_jpim1 ! vector opt. 311 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 312 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 313 END DO 314 END DO 315 ENDIF 316 319 zbfru(ji,jj) = MAX( bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef ) 320 zbfrv(ji,jj) = MAX( bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef ) 321 END DO 322 END DO 323 324 IF( lk_vvl ) THEN 325 DO jj = 2, jpjm1 326 DO ji = fs_2, fs_jpim1 ! vector opt. 327 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) & 328 & / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) ) 329 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) & 330 & / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) ) 331 END DO 332 END DO 333 ELSE 334 DO jj = 2, jpjm1 335 DO ji = fs_2, fs_jpim1 ! vector opt. 336 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 337 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 338 END DO 339 END DO 340 ENDIF 341 END IF ! end (ln_bfrimp) 342 343 317 344 ! !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 318 345 ! ! -------------------------- … … 321 348 ! 322 349 IF( lk_vvl ) THEN 323 ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1. e0- umask(:,:,1) )324 vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1. e0- vmask(:,:,1) )350 ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 351 vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 325 352 ELSE 326 353 ub_b(:,:) = ub_b(:,:) * hur(:,:) … … 358 385 ! set ssh corrections to 0 359 386 ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 360 IF( lp_obc_east ) sshfoe_b(:,:) = 0. e0361 IF( lp_obc_west ) sshfow_b(:,:) = 0. e0362 IF( lp_obc_south ) sshfos_b(:,:) = 0. e0363 IF( lp_obc_north ) sshfon_b(:,:) = 0. e0387 IF( lp_obc_east ) sshfoe_b(:,:) = 0._wp 388 IF( lp_obc_west ) sshfow_b(:,:) = 0._wp 389 IF( lp_obc_south ) sshfos_b(:,:) = 0._wp 390 IF( lp_obc_north ) sshfon_b(:,:) = 0._wp 364 391 #endif 365 392 … … 377 404 ! !* after ssh_e 378 405 ! ! ----------- 379 DO jj = 2, jpjm1 406 DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports 380 407 DO ji = fs_2, fs_jpim1 ! vector opt. 381 408 zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) & … … 389 416 ! ! OBC : zhdiv must be zero behind the open boundary 390 417 !! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 391 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0. e0! east392 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0. e0! west393 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0. e0! north394 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0. e0! south418 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0._wp ! east 419 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0._wp ! west 420 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0._wp ! north 421 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0._wp ! south 395 422 #endif 396 423 #if defined key_bdy … … 406 433 ! !* after barotropic velocities (vorticity scheme dependent) 407 434 ! ! --------------------------- 408 zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) 435 zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! now_e transport 409 436 zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 410 437 ! … … 430 457 zv_cor =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 431 458 ! after velocities with implicit bottom friction 432 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) & 433 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 434 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) & 435 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 459 460 IF( ln_bfrimp ) THEN ! implicit bottom friction 461 ! A new method to implement the implicit bottom friction. 462 ! H. Liu 463 ! Sept 2011 464 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 465 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 466 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 467 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 468 ! 469 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 470 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 471 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 472 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 473 ! 474 ELSE 475 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) & 476 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 477 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) & 478 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 479 ENDIF 436 480 END DO 437 481 END DO … … 456 500 zv_cor = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj) 457 501 ! after velocities with implicit bottom friction 458 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) & 459 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 460 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) & 461 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 502 IF( ln_bfrimp ) THEN 503 ! A new method to implement the implicit bottom friction. 504 ! H. Liu 505 ! Sept 2011 506 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 507 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 508 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 509 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 510 ! 511 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 512 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 513 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 514 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 515 ! 516 ELSE 517 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) & 518 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 519 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) & 520 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 521 ENDIF 462 522 END DO 463 523 END DO … … 482 542 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj) 483 543 ! after velocities with implicit bottom friction 484 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) & 485 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 486 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) & 487 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 544 IF( ln_bfrimp ) THEN 545 ! A new method to implement the implicit bottom friction. 546 ! H. Liu 547 ! Sept 2011 548 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 549 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 550 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 551 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 552 ! 553 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 554 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 555 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 556 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 557 ! 558 ELSE 559 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) & 560 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 561 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) & 562 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 563 ENDIF 488 564 END DO 489 565 END DO … … 492 568 ! !* domain lateral boundary 493 569 ! ! ----------------------- 494 ! 495 ! 496 ! 570 ! ! Flather's boundary condition for the barotropic loop : 571 ! ! - Update sea surface height on each open boundary 572 ! ! - Correct the velocity 497 573 498 574 IF( lk_obc ) CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e ) … … 529 605 DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points 530 606 DO ji = 1, fs_jpim1 ! Vector opt. 531 zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &532 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) &533 & +e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )534 zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &535 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) &536 & +e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )607 zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 608 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & 609 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 610 zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 611 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & 612 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 537 613 END DO 538 614 END DO … … 542 618 hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points 543 619 hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 544 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1. e0- umask(:,:,1) )545 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1. e0- vmask(:,:,1) )620 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 621 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 546 622 ! 547 623 ENDIF … … 562 638 ! 563 639 ! !* Time average ==> after barotropic u, v, ssh 564 zcoef = 1. e0/ ( 2 * nn_baro + 1 )640 zcoef = 1._wp / ( 2 * nn_baro + 1 ) 565 641 zu_sum(:,:) = zcoef * zu_sum (:,:) 566 642 zv_sum(:,:) = zcoef * zv_sum (:,:) … … 603 679 CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 604 680 ELSE 605 un_b (:,:) = 0. e0606 vn_b (:,:) = 0. e0681 un_b (:,:) = 0._wp 682 vn_b (:,:) = 0._wp 607 683 ! vertical sum 608 684 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll … … 625 701 ! Vertically integrated velocity (before) 626 702 IF (neuler/=0) THEN 627 ub_b (:,:) = 0. e0628 vb_b (:,:) = 0. e0703 ub_b (:,:) = 0._wp 704 vb_b (:,:) = 0._wp 629 705 630 706 ! vertical sum … … 644 720 645 721 IF( lk_vvl ) THEN 646 ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1. e0- umask(:,:,1) )647 vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1. e0- vmask(:,:,1) )722 ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 723 vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 648 724 ELSE 649 725 ub_b(:,:) = ub_b(:,:) * hur(:,:) -
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r2715 r3072 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_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
r2715 r3072 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
Note: See TracChangeset
for help on using the changeset viewer.
