- Timestamp:
- 2010-10-04T15:53:42+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r1662 r2148 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 , LOCEAN-IPSL (2010) 38 33 !! $Id$ 39 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt34 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 40 35 !!---------------------------------------------------------------------- 41 36 42 37 CONTAINS 43 44 38 45 39 SUBROUTINE dyn_zdf_imp( kt, p2dt ) … … 54 48 !! dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) ) 55 49 !! backward time stepping 56 !! Surface boundary conditions: wind stress input 50 !! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) 57 51 !! Bottom boundary conditions : bottom stress (cf zdfbfr.F) 58 52 !! Add this trend to the general trend ua : 59 53 !! ua = ua + dz( avmu dz(u) ) 60 54 !! 61 !! ** Action : - Update (ua,va) arrays with the after vertical diffusive 62 !! mixing trend. 55 !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend. 63 56 !!--------------------------------------------------------------------- 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 57 USE oce, ONLY : zwd => ta ! use ta as workspace 58 USE oce, ONLY : zws => sa ! use sa as workspace 59 !! 60 INTEGER , INTENT( in ) :: kt ! ocean time-step index 61 REAL(wp), INTENT( in ) :: p2dt ! vertical profile of tracer time-step 62 !! 63 INTEGER :: ji, jj, jk ! dummy loop indices 64 REAL(wp) :: zrau0r, zcoef ! temporary scalars 65 REAL(wp) :: zzwi, zzws, zrhs ! temporary scalars 66 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi ! 3D workspace 77 67 !!---------------------------------------------------------------------- 78 68 … … 95 85 ! friction velocity in dyn_bfr using values from the previous timestep. There 96 86 ! is no need to include these in the implicit calculation. 97 DO jk = 1, jpkm1 87 ! 88 DO jk = 1, jpkm1 ! Matrix 98 89 DO jj = 2, jpjm1 99 90 DO ji = fs_2, fs_jpim1 ! vector opt. 100 91 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)92 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 93 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) 94 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 95 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1) 105 96 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws 106 97 END DO 107 98 END DO 108 99 END DO 109 110 ! Surface boudary conditions 111 DO jj = 2, jpjm1 100 DO jj = 2, jpjm1 ! Surface boudary conditions 112 101 DO ji = fs_2, fs_jpim1 ! vector opt. 113 102 zwi(ji,jj,1) = 0. … … 130 119 ! The solution (the after velocity) is in ua 131 120 !----------------------------------------------------------------------- 132 133 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 134 DO jk = 2, jpkm1 121 ! 122 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 135 123 DO jj = 2, jpjm1 136 124 DO ji = fs_2, fs_jpim1 ! vector opt. … … 139 127 END DO 140 128 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) 129 ! 130 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 131 DO ji = fs_2, fs_jpim1 ! vector opt. 132 ua(ji,jj,1) = ub(ji,jj,1) & 133 & + p2dt * ( ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 ) ) 151 134 END DO 152 135 END DO … … 159 142 END DO 160 143 END DO 161 162 ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 163 DO jj = 2, jpjm1 144 ! 145 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk == 164 146 DO ji = fs_2, fs_jpim1 ! vector opt. 165 147 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 173 155 END DO 174 156 END DO 175 176 ! Normalization to obtain the general momentum trend ua 177 DO jk = 1, jpkm1 157 ! 158 DO jk = 1, jpkm1 !== Normalization to obtain the general momentum trend ua == 178 159 DO jj = 2, jpjm1 179 160 DO ji = fs_2, fs_jpim1 ! vector opt. … … 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 184 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws … … 204 186 END DO 205 187 END DO 206 207 ! Surface boudary conditions 208 DO jj = 2, jpjm1 188 DO jj = 2, jpjm1 ! Surface boudary conditions 209 189 DO ji = fs_2, fs_jpim1 ! vector opt. 210 190 zwi(ji,jj,1) = 0.e0 … … 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) & 221 & + p2dt * ( va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( 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) … … 271 243 END DO 272 244 END DO 273 274 ! Normalization to obtain the general momentum trend va 275 DO jk = 1, jpkm1 245 ! 246 DO jk = 1, jpkm1 !== Normalization to obtain the general momentum trend va == 276 247 DO jj = 2, jpjm1 277 248 DO ji = fs_2, fs_jpim1 ! vector opt. … … 280 251 END DO 281 252 END DO 282 253 ! 283 254 END SUBROUTINE dyn_zdf_imp 284 255
Note: See TracChangeset
for help on using the changeset viewer.
