- Timestamp:
- 2020-01-27T15:31:53+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traldf_iso.F90
r12193 r12340 41 41 !! * Substitutions 42 42 # include "vectopt_loop_substitute.h90" 43 # include "do_loop_substitute.h90" 43 44 !!---------------------------------------------------------------------- 44 45 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 145 146 IF( kpass == 1 ) THEN !== first pass only ==! 146 147 ! 147 DO jk = 2, jpkm1 148 DO jj = 2, jpjm1 149 DO ji = fs_2, fs_jpim1 ! vector opt. 150 ! 151 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 152 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 153 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 154 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 155 ! 156 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 157 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 158 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 159 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 160 ! 161 ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 162 & + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 163 END DO 164 END DO 165 END DO 148 DO_3D_00_00( 2, jpkm1 ) 149 ! 150 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 151 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 152 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 153 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 154 ! 155 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 156 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 157 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 158 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 159 ! 160 ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 161 & + zahv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 162 END_3D 166 163 ! 167 164 IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient 168 DO jk = 2, jpkm1 169 DO jj = 2, jpjm1 170 DO ji = fs_2, fs_jpim1 171 akz(ji,jj,jk) = 0.25_wp * ( & 172 & ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & 173 & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & 174 & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & 175 & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) 176 END DO 177 END DO 178 END DO 165 DO_3D_00_00( 2, jpkm1 ) 166 akz(ji,jj,jk) = 0.25_wp * ( & 167 & ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & 168 & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & 169 & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & 170 & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) 171 END_3D 179 172 ! 180 173 IF( ln_traldf_blp ) THEN ! bilaplacian operator 181 DO jk = 2, jpkm1 182 DO jj = 1, jpjm1 183 DO ji = 1, fs_jpim1 184 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 185 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) ) ) 186 END DO 187 END DO 188 END DO 174 DO_3D_10_10( 2, jpkm1 ) 175 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 176 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) ) ) 177 END_3D 189 178 ELSEIF( ln_traldf_lap ) THEN ! laplacian operator 190 DO jk = 2, jpkm1 191 DO jj = 1, jpjm1 192 DO ji = 1, fs_jpim1 193 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 194 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 195 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 196 END DO 197 END DO 198 END DO 179 DO_3D_10_10( 2, jpkm1 ) 180 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 181 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 182 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt 183 END_3D 199 184 ENDIF 200 185 ! … … 217 202 218 203 ! Horizontal tracer gradient 219 DO jk = 1, jpkm1 220 DO jj = 1, jpjm1 221 DO ji = 1, fs_jpim1 ! vector opt. 222 zdit(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 223 zdjt(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 224 END DO 225 END DO 226 END DO 204 DO_3D_10_10( 1, jpkm1 ) 205 zdit(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 206 zdjt(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 207 END_3D 227 208 IF( ln_zps ) THEN ! botton and surface ocean correction of the horizontal gradient 228 DO jj = 1, jpjm1 ! bottom correction (partial bottom cell) 229 DO ji = 1, fs_jpim1 ! vector opt. 230 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 231 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 232 END DO 233 END DO 209 DO_2D_10_10 210 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 211 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 212 END_2D 234 213 IF( ln_isfcav ) THEN ! first wet level beneath a cavity 235 DO jj = 1, jpjm1 236 DO ji = 1, fs_jpim1 ! vector opt. 237 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 238 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 239 END DO 240 END DO 214 DO_2D_10_10 215 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 216 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) 217 END_2D 241 218 ENDIF 242 219 ENDIF … … 254 231 ELSE ; zdkt(:,:) = ( pt(:,:,jk-1,jn) - pt(:,:,jk,jn) ) * wmask(:,:,jk) 255 232 ENDIF 256 DO jj = 1 , jpjm1 !== Horizontal fluxes 257 DO ji = 1, fs_jpim1 ! vector opt. 258 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) 259 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 260 ! 261 zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & 262 & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) 263 ! 264 zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & 265 & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) 266 ! 267 zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 268 zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 269 ! 270 zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & 271 & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & 272 & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) 273 zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & 274 & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & 275 & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) 276 END DO 277 END DO 278 ! 279 DO jj = 2 , jpjm1 !== horizontal divergence and add to pt_rhs 280 DO ji = fs_2, fs_jpim1 ! vector opt. 281 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & 282 & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & 283 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 284 END DO 285 END DO 233 DO_2D_10_10 234 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) 235 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 236 ! 237 zmsku = 1. / MAX( wmask(ji+1,jj,jk ) + wmask(ji,jj,jk+1) & 238 & + wmask(ji+1,jj,jk+1) + wmask(ji,jj,jk ), 1. ) 239 ! 240 zmskv = 1. / MAX( wmask(ji,jj+1,jk ) + wmask(ji,jj,jk+1) & 241 & + wmask(ji,jj+1,jk+1) + wmask(ji,jj,jk ), 1. ) 242 ! 243 zcof1 = - pahu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku 244 zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 245 ! 246 zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & 247 & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & 248 & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) 249 zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & 250 & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & 251 & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) 252 END_2D 253 ! 254 DO_2D_00_00 255 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & 256 & + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & 257 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 258 END_2D 286 259 END DO ! End of slab 287 260 … … 297 270 ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp 298 271 299 DO jk = 2, jpkm1 ! interior (2=<jk=<jpk-1) 300 DO jj = 2, jpjm1 301 DO ji = fs_2, fs_jpim1 ! vector opt. 302 ! 303 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 304 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 305 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 306 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 307 ! 308 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 309 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 310 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 311 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 312 ! 313 zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) !wslpi & j are already w-masked 314 zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 315 ! 316 ztfw(ji,jj,jk) = zcoef3 * ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) & 317 & + zdit(ji-1,jj ,jk-1) + zdit(ji ,jj ,jk) ) & 318 & + zcoef4 * ( zdjt(ji ,jj ,jk-1) + zdjt(ji ,jj-1,jk) & 319 & + zdjt(ji ,jj-1,jk-1) + zdjt(ji ,jj ,jk) ) 320 END DO 321 END DO 322 END DO 272 DO_3D_00_00( 2, jpkm1 ) 273 ! 274 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & 275 & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) 276 zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & 277 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 278 ! 279 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 280 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 281 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 282 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 283 ! 284 zcoef3 = - zahu_w * e2t(ji,jj) * zmsku * wslpi (ji,jj,jk) !wslpi & j are already w-masked 285 zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 286 ! 287 ztfw(ji,jj,jk) = zcoef3 * ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) & 288 & + zdit(ji-1,jj ,jk-1) + zdit(ji ,jj ,jk) ) & 289 & + zcoef4 * ( zdjt(ji ,jj ,jk-1) + zdjt(ji ,jj-1,jk) & 290 & + zdjt(ji ,jj-1,jk-1) + zdjt(ji ,jj ,jk) ) 291 END_3D 323 292 ! !== add the vertical 33 flux ==! 324 293 IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz 325 DO jk = 2, jpkm1 326 DO jj = 2, jpjm1 327 DO ji = fs_2, fs_jpim1 328 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) & 329 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 330 & * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 331 END DO 332 END DO 333 END DO 294 DO_3D_00_00( 2, jpkm1 ) 295 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) & 296 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 297 & * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 298 END_3D 334 299 ! 335 300 ELSE ! bilaplacian 336 301 SELECT CASE( kpass ) 337 302 CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2 338 DO jk = 2, jpkm1 339 DO jj = 2, jpjm1 340 DO ji = fs_2, fs_jpim1 341 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) & 342 & + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & 343 & * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 344 END DO 345 END DO 346 END DO 303 DO_3D_00_00( 2, jpkm1 ) 304 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) & 305 & + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & 306 & * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 307 END_3D 347 308 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt and pt2 gradients, resp. 348 DO jk = 2, jpkm1 349 DO jj = 2, jpjm1 350 DO ji = fs_2, fs_jpim1 351 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) & 352 & * ( ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) ) & 353 & + akz(ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) ) ) 354 END DO 355 END DO 356 END DO 309 DO_3D_00_00( 2, jpkm1 ) 310 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) & 311 & * ( ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) ) & 312 & + akz(ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) ) ) 313 END_3D 357 314 END SELECT 358 315 ENDIF 359 316 ! 360 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pt_rhs ==! 361 DO jj = 2, jpjm1 362 DO ji = fs_2, fs_jpim1 ! vector opt. 363 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) & 364 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 365 END DO 366 END DO 367 END DO 317 DO_3D_00_00( 1, jpkm1 ) 318 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) & 319 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 320 END_3D 368 321 ! 369 322 IF( ( kpass == 1 .AND. ln_traldf_lap ) .OR. & !== first pass only ( laplacian) ==!
Note: See TracChangeset
for help on using the changeset viewer.