- Timestamp:
- 2013-11-20T17:28:04+01:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r3625 r4292 16 16 USE oce ! ocean dynamics and tracers 17 17 USE dom_oce ! ocean space and time domain 18 USE domvvl ! variable volume 18 19 USE sbc_oce ! surface boundary condition: ocean 19 20 USE zdf_oce ! ocean vertical physics … … 24 25 USE wrk_nemo ! Memory Allocation 25 26 USE timing ! Timing 27 USE dynadv ! dynamics: vector invariant versus flux form 28 USE dynspg_oce, ONLY: lk_dynspg_ts, ua_b, va_b 29 USE dynspg_ts 26 30 27 31 IMPLICIT NONE … … 29 33 30 34 PUBLIC dyn_zdf_imp ! called by step.F90 35 36 REAL(wp) :: r_vvl ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise 31 37 32 38 !! * Substitutions … … 64 70 INTEGER :: ikbu, ikbv ! local integers 65 71 REAL(wp) :: z1_p2dt, zcoef, zzwi, zzws, zrhs ! local scalars 72 REAL(wp) :: ze3ua, ze3va 66 73 !!---------------------------------------------------------------------- 67 74 68 75 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwd, zws 69 REAL(wp), POINTER, DIMENSION(:,:) :: zavmu, zavmv70 76 !!---------------------------------------------------------------------- 71 77 ! … … 73 79 ! 74 80 CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws ) 75 CALL wrk_alloc( jpi,jpj, zavmu, zavmv )76 81 ! 77 82 IF( kt == nit000 ) THEN … … 79 84 IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 80 85 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 86 ! 87 IF( lk_vvl ) THEN ; r_vvl = 1._wp ! Variable volume indicator 88 ELSE ; r_vvl = 0._wp 89 ENDIF 81 90 ENDIF 82 91 … … 94 103 IF( ln_bfrimp ) THEN 95 104 # if defined key_vectopt_loop 96 DO jj = 1, 197 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)105 DO jj = 1, 1 106 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 98 107 # else 99 DO jj = 2, jpjm1100 DO ji = 2, jpim1108 DO jj = 2, jpjm1 109 DO ji = 2, jpim1 101 110 # endif 102 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 103 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 104 zavmu(ji,jj) = avmu(ji,jj,ikbu+1) 105 zavmv(ji,jj) = avmv(ji,jj,ikbv+1) 106 avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 107 avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 108 END DO 109 END DO 110 ENDIF 111 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 112 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 113 avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 114 avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 115 END DO 116 END DO 117 ENDIF 118 119 #if defined key_dynspg_ts 120 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 121 DO jk = 1, jpkm1 122 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 123 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 124 END DO 125 ELSE ! applied on thickness weighted velocity 126 DO jk = 1, jpkm1 127 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 128 & + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 129 & / fse3u_a(:,:,jk) * umask(:,:,jk) 130 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 131 & + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 132 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 133 END DO 134 ENDIF 135 136 IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN 137 ! remove barotropic velocities: 138 DO jk = 1, jpkm1 139 ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk) 140 va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk) 141 ENDDO 142 ! Add bottom stress due to barotropic component only: 143 DO jj = 2, jpjm1 144 DO ji = fs_2, fs_jpim1 ! vector opt. 145 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 146 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 147 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl * fse3u_a(ji,jj,ikbu) 148 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl * fse3v_a(ji,jj,ikbv) 149 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 150 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 151 END DO 152 END DO 153 ENDIF 154 #endif 111 155 112 156 ! 2. Vertical diffusion on u … … 119 163 DO jj = 2, jpjm1 120 164 DO ji = fs_2, fs_jpim1 ! vector opt. 121 zcoef = - p2dt / fse3u(ji,jj,jk) 165 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl * fse3u_a(ji,jj,jk) ! after scale factor at T-point 166 zcoef = - p2dt / ze3ua 122 167 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 123 168 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) … … 128 173 END DO 129 174 END DO 130 DO jj = 2, jpjm1 ! Surface bou dary conditions175 DO jj = 2, jpjm1 ! Surface boundary conditions 131 176 DO ji = fs_2, fs_jpim1 ! vector opt. 132 177 zwi(ji,jj,1) = 0._wp … … 160 205 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 161 206 DO ji = fs_2, fs_jpim1 ! vector opt. 162 ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ( ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 163 & * r1_rau0 / fse3u(ji,jj,1) ) 207 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1) 208 #if defined key_dynspg_ts 209 ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 210 & / ( ze3ua * rau0 ) 211 #else 212 ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 213 & / ( fse3u(ji,jj,1) * rau0 ) ) 214 #endif 164 215 END DO 165 216 END DO … … 167 218 DO jj = 2, jpjm1 168 219 DO ji = fs_2, fs_jpim1 ! vector opt. 169 zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) ! zrhs=right hand side 220 #if defined key_dynspg_ts 221 zrhs = ua(ji,jj,jk) ! zrhs=right hand side 222 #else 223 zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) 224 #endif 170 225 ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 171 226 END DO … … 186 241 END DO 187 242 243 #if ! defined key_dynspg_ts 188 244 ! Normalization to obtain the general momentum trend ua 189 245 DO jk = 1, jpkm1 … … 194 250 END DO 195 251 END DO 196 252 #endif 197 253 198 254 ! 3. Vertical diffusion on v … … 205 261 DO jj = 2, jpjm1 206 262 DO ji = fs_2, fs_jpim1 ! vector opt. 207 zcoef = -p2dt / fse3v(ji,jj,jk) 263 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk) + r_vvl * fse3v_a(ji,jj,jk) ! after scale factor at T-point 264 zcoef = - p2dt / ze3va 208 265 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk ) 209 266 zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk) … … 214 271 END DO 215 272 END DO 216 DO jj = 2, jpjm1 ! Surface bou dary conditions273 DO jj = 2, jpjm1 ! Surface boundary conditions 217 274 DO ji = fs_2, fs_jpim1 ! vector opt. 218 275 zwi(ji,jj,1) = 0._wp … … 246 303 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 247 304 DO ji = fs_2, fs_jpim1 ! vector opt. 248 va(ji,jj,1) = vb(ji,jj,1) + p2dt * ( va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 249 & * r1_rau0 / fse3v(ji,jj,1) ) 305 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl * fse3v_a(ji,jj,1) 306 #if defined key_dynspg_ts 307 va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 308 & / ( ze3va * rau0 ) 309 #else 310 va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 311 & / ( fse3v(ji,jj,1) * rau0 ) ) 312 #endif 250 313 END DO 251 314 END DO … … 253 316 DO jj = 2, jpjm1 254 317 DO ji = fs_2, fs_jpim1 ! vector opt. 255 zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) ! zrhs=right hand side 318 #if defined key_dynspg_ts 319 zrhs = va(ji,jj,jk) ! zrhs=right hand side 320 #else 321 zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) 322 #endif 256 323 va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 257 324 END DO … … 259 326 END DO 260 327 ! 261 DO jj = 2, jpjm1 !== th rid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==328 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk == 262 329 DO ji = fs_2, fs_jpim1 ! vector opt. 263 330 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 273 340 274 341 ! Normalization to obtain the general momentum trend va 342 #if ! defined key_dynspg_ts 275 343 DO jk = 1, jpkm1 276 344 DO jj = 2, jpjm1 … … 280 348 END DO 281 349 END DO 282 350 #endif 351 352 ! J. Chanut: Lines below are useless ? 283 353 !! restore bottom layer avmu(v) 284 354 IF( ln_bfrimp ) THEN … … 292 362 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 293 363 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 294 avmu(ji,jj,ikbu+1) = zavmu(ji,jj)295 avmv(ji,jj,ikbv+1) = zavmv(ji,jj)364 avmu(ji,jj,ikbu+1) = 0.e0 365 avmv(ji,jj,ikbv+1) = 0.e0 296 366 END DO 297 367 END DO … … 299 369 ! 300 370 CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 301 CALL wrk_dealloc( jpi,jpj, zavmu, zavmv)302 371 ! 303 372 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_imp')
Note: See TracChangeset
for help on using the changeset viewer.
