Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
- Property svn:eol-style deleted
r1662 r2528 4 4 !! Ocean dynamics: vertical component(s) of the momentum mixing trend 5 5 !!============================================================================== 6 !! History : ! 90-10 (B. Blanke) Original code 7 !! ! 97-05 (G. Madec) vertical component of isopycnal 8 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 6 !! History : OPA ! 1990-10 (B. Blanke) Original code 7 !! 8.0 ! 1997-05 (G. Madec) vertical component of isopycnal 8 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 9 10 !!---------------------------------------------------------------------- 10 11 … … 13 14 !! sion using a implicit time-stepping. 14 15 !!---------------------------------------------------------------------- 15 !! OPA 9.0 , LOCEAN-IPSL (2005)16 !! $Id$17 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)18 !!----------------------------------------------------------------------19 !! * Modules used20 16 USE oce ! ocean dynamics and tracers 21 17 USE dom_oce ! ocean space and time domain … … 28 24 PRIVATE 29 25 30 !! * Routine accessibility 31 PUBLIC dyn_zdf_imp ! called by step.F90 26 PUBLIC dyn_zdf_imp ! called by step.F90 32 27 33 28 !! * Substitutions … … 35 30 # include "vectopt_loop_substitute.h90" 36 31 !!---------------------------------------------------------------------- 37 !! OPA 9.0 , LOCEAN-IPSL (2005)32 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 38 33 !! $Id$ 39 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 40 !!---------------------------------------------------------------------- 41 34 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 35 !!---------------------------------------------------------------------- 42 36 CONTAINS 43 44 37 45 38 SUBROUTINE dyn_zdf_imp( kt, p2dt ) … … 54 47 !! dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) ) 55 48 !! backward time stepping 56 !! Surface boundary conditions: wind stress input 49 !! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) 57 50 !! Bottom boundary conditions : bottom stress (cf zdfbfr.F) 58 51 !! Add this trend to the general trend ua : 59 52 !! ua = ua + dz( avmu dz(u) ) 60 53 !! 61 !! ** Action : - Update (ua,va) arrays with the after vertical diffusive 62 !! mixing trend. 54 !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend. 63 55 !!--------------------------------------------------------------------- 64 !! * Modules used 65 USE oce, ONLY : zwd => ta, & ! use ta as workspace 66 zws => sa ! use sa as workspace 67 68 !! * Arguments 69 INTEGER , INTENT( in ) :: kt ! ocean time-step index 70 REAL(wp), INTENT( in ) :: p2dt ! vertical profile of tracer time-step 71 72 !! * Local declarations 73 INTEGER :: ji, jj, jk ! dummy loop indices 74 REAL(wp) :: zrau0r, z2dtf, zcoef, zzws, zrhs ! temporary scalars 75 REAL(wp) :: zzwi ! temporary scalars 76 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi ! temporary workspace arrays 56 USE oce, ONLY : zwd => ta ! use ta as workspace 57 USE oce, ONLY : zws => sa ! use sa as workspace 58 !! 59 INTEGER , INTENT( in ) :: kt ! ocean time-step index 60 REAL(wp), INTENT( in ) :: p2dt ! vertical profile of tracer time-step 61 !! 62 INTEGER :: ji, jj, jk ! dummy loop indices 63 REAL(wp) :: z1_p2dt, zcoef ! temporary scalars 64 REAL(wp) :: zzwi, zzws, zrhs ! temporary scalars 65 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi ! 3D workspace 77 66 !!---------------------------------------------------------------------- 78 67 … … 85 74 ! 0. Local constant initialization 86 75 ! -------------------------------- 87 z rau0r = 1. / rau0 ! inverse of the reference density76 z1_p2dt = 1._wp / p2dt ! inverse of the timestep 88 77 89 78 ! 1. Vertical diffusion on u … … 95 84 ! friction velocity in dyn_bfr using values from the previous timestep. There 96 85 ! is no need to include these in the implicit calculation. 97 DO jk = 1, jpkm1 86 ! 87 DO jk = 1, jpkm1 ! Matrix 98 88 DO jj = 2, jpjm1 99 89 DO ji = fs_2, fs_jpim1 ! vector opt. 100 90 zcoef = - p2dt / fse3u(ji,jj,jk) 101 zzwi = zcoef * avmu(ji,jj,jk ) / fse3uw(ji,jj,jk ) 102 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) 103 zzws = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 104 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1) 105 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws 106 END DO 107 END DO 108 END DO 109 110 ! Surface boudary conditions 111 DO jj = 2, jpjm1 112 DO ji = fs_2, fs_jpim1 ! vector opt. 113 zwi(ji,jj,1) = 0. 114 zwd(ji,jj,1) = 1. - zws(ji,jj,1) 91 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 92 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) 93 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 94 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1) 95 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 96 END DO 97 END DO 98 END DO 99 DO jj = 2, jpjm1 ! Surface boudary conditions 100 DO ji = fs_2, fs_jpim1 ! vector opt. 101 zwi(ji,jj,1) = 0._wp 102 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 115 103 END DO 116 104 END DO … … 130 118 ! The solution (the after velocity) is in ua 131 119 !----------------------------------------------------------------------- 132 133 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 134 DO jk = 2, jpkm1 120 ! 121 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 135 122 DO jj = 2, jpjm1 136 123 DO ji = fs_2, fs_jpim1 ! vector opt. … … 139 126 END DO 140 127 END DO 141 142 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 143 DO jj = 2, jpjm1 144 DO ji = fs_2, fs_jpim1 ! vector opt. 145 !!! change les resultats (derniers digit, pas significativement + rapide 1* de moins) 146 !!! ua(ji,jj,1) = ub(ji,jj,1) & 147 !!! + p2dt * ( ua(ji,jj,1) + utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) ) 148 z2dtf = p2dt / ( fse3u(ji,jj,1)*rau0 ) 149 ua(ji,jj,1) = ub(ji,jj,1) & 150 + p2dt * ua(ji,jj,1) + z2dtf * utau(ji,jj) 128 ! 129 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 130 DO ji = fs_2, fs_jpim1 ! vector opt. 131 ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ( ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 132 & / ( fse3u(ji,jj,1) * rau0 ) ) 151 133 END DO 152 134 END DO … … 159 141 END DO 160 142 END DO 161 162 ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 163 DO jj = 2, jpjm1 143 ! 144 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk == 164 145 DO ji = fs_2, fs_jpim1 ! vector opt. 165 146 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 169 150 DO jj = 2, jpjm1 170 151 DO ji = fs_2, fs_jpim1 ! vector opt. 171 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)152 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 172 153 END DO 173 154 END DO … … 178 159 DO jj = 2, jpjm1 179 160 DO ji = fs_2, fs_jpim1 ! vector opt. 180 ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) /p2dt161 ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt 181 162 END DO 182 163 END DO … … 192 173 ! friction velocity in dyn_bfr using values from the previous timestep. There 193 174 ! is no need to include these in the implicit calculation. 194 DO jk = 1, jpkm1 175 ! 176 DO jk = 1, jpkm1 ! Matrix 195 177 DO jj = 2, jpjm1 196 178 DO ji = fs_2, fs_jpim1 ! vector opt. 197 179 zcoef = -p2dt / fse3v(ji,jj,jk) 198 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk )180 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk ) 199 181 zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk) 200 zzws = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)182 zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 201 183 zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1) 202 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws 203 END DO 204 END DO 205 END DO 206 207 ! Surface boudary conditions 208 DO jj = 2, jpjm1 209 DO ji = fs_2, fs_jpim1 ! vector opt. 210 zwi(ji,jj,1) = 0.e0 211 zwd(ji,jj,1) = 1. - zws(ji,jj,1) 184 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 185 END DO 186 END DO 187 END DO 188 DO jj = 2, jpjm1 ! Surface boudary conditions 189 DO ji = fs_2, fs_jpim1 ! vector opt. 190 zwi(ji,jj,1) = 0._wp 191 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 212 192 END DO 213 193 END DO … … 223 203 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 224 204 ! 225 ! m is decomposed in the product of an upper and lower triangular 226 ! matrix 205 ! m is decomposed in the product of an upper and lower triangular matrix 227 206 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 228 207 ! The solution (after velocity) is in 2d array va 229 208 !----------------------------------------------------------------------- 230 231 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 232 DO jk = 2, jpkm1 209 ! 210 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 233 211 DO jj = 2, jpjm1 234 212 DO ji = fs_2, fs_jpim1 ! vector opt. … … 237 215 END DO 238 216 END DO 239 240 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 241 DO jj = 2, jpjm1 242 DO ji = fs_2, fs_jpim1 ! vector opt. 243 !!! change les resultats (derniers digit, pas significativement + rapide 1* de moins) 244 !!! va(ji,jj,1) = vb(ji,jj,1) & 245 !!! + p2dt * ( va(ji,jj,1) + vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) ) 246 z2dtf = p2dt / ( fse3v(ji,jj,1)*rau0 ) 247 va(ji,jj,1) = vb(ji,jj,1) & 248 + p2dt * va(ji,jj,1) + z2dtf * vtau(ji,jj) 217 ! 218 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 219 DO ji = fs_2, fs_jpim1 ! vector opt. 220 va(ji,jj,1) = vb(ji,jj,1) + p2dt * ( va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 221 & / ( fse3v(ji,jj,1) * rau0 ) ) 249 222 END DO 250 223 END DO … … 257 230 END DO 258 231 END DO 259 260 ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk 261 DO jj = 2, jpjm1 232 ! 233 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk == 262 234 DO ji = fs_2, fs_jpim1 ! vector opt. 263 235 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 267 239 DO jj = 2, jpjm1 268 240 DO ji = fs_2, fs_jpim1 ! vector opt. 269 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)241 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 270 242 END DO 271 243 END DO … … 276 248 DO jj = 2, jpjm1 277 249 DO ji = fs_2, fs_jpim1 ! vector opt. 278 va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) /p2dt279 END DO 280 END DO 281 END DO 282 250 va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt 251 END DO 252 END DO 253 END DO 254 ! 283 255 END SUBROUTINE dyn_zdf_imp 284 256
Note: See TracChangeset
for help on using the changeset viewer.
