- Timestamp:
- 2016-11-28T17:04:10+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_triad.F90
r5836 r7351 11 11 !! tra_ldf_triad : update the tracer trend with the iso-neutral laplacian triad-operator 12 12 !!---------------------------------------------------------------------- 13 USE oce 14 USE dom_oce 15 USE phycst 16 USE trc_oce 17 USE zdf_oce 18 USE ldftra 19 USE ldfslp 20 USE traldf_iso 21 USE diaptr 22 USE zpshde 13 USE oce ! ocean dynamics and active tracers 14 USE dom_oce ! ocean space and time domain 15 USE phycst ! physical constants 16 USE trc_oce ! share passive tracers/Ocean variables 17 USE zdf_oce ! ocean vertical physics 18 USE ldftra ! lateral physics: eddy diffusivity 19 USE ldfslp ! lateral physics: iso-neutral slopes 20 USE traldf_iso ! lateral diffusion (Madec operator) (tra_ldf_iso routine) 21 USE diaptr ! poleward transport diagnostics 22 USE zpshde ! partial step: hor. derivative (zps_hde routine) 23 23 ! 24 USE in_out_manager 25 USE iom 26 USE lbclnk 27 USE lib_mpp 28 USE wrk_nemo 29 USE timing 24 USE in_out_manager ! I/O manager 25 USE iom ! I/O library 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 27 USE lib_mpp ! MPP library 28 USE wrk_nemo ! Memory Allocation 29 USE timing ! Timing 30 30 31 31 IMPLICIT NONE … … 37 37 38 38 !! * Substitutions 39 # include "domzgr_substitute.h90"40 39 # include "vectopt_loop_substitute.h90" 41 40 !!---------------------------------------------------------------------- … … 113 112 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 114 113 ENDIF 115 !116 114 ! ! set time step size (Euler/Leapfrog) 117 IF( neuler == 0 .AND. kt == kit000 ) THEN ; z2dt = rdt tra(1)! at nit000 (Euler)118 ELSE ; z2dt = 2.* rdt tra(1)! (Leapfrog)115 IF( neuler == 0 .AND. kt == kit000 ) THEN ; z2dt = rdt ! at nit000 (Euler) 116 ELSE ; z2dt = 2.* rdt ! (Leapfrog) 119 117 ENDIF 120 118 z1_2dt = 1._wp / z2dt … … 123 121 ELSE ; zsign = -1._wp 124 122 ENDIF 125 123 ! 126 124 !!---------------------------------------------------------------------- 127 125 !! 0 - calculate ah_wslp2, akz, and optionally zpsi_uw, zpsi_vw … … 142 140 DO jj = 1, jpjm1 143 141 DO ji = 1, fs_jpim1 144 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)145 zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk)142 ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 143 zbu = e1e2u(ji,jj) * e3u_n(ji,jj,jk) 146 144 zah = 0.25_wp * pahu(ji,jj,jk) 147 145 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 148 146 ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 149 zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)147 zslope2 = zslope_skew + ( gdept_n(ji+1,jj,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 150 148 zslope2 = zslope2 *zslope2 151 149 ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2 152 150 akz (ji+ip,jj,jk+kp) = akz (ji+ip,jj,jk+kp) + zah * r1_e1u(ji,jj) & 153 151 & * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 154 !152 ! 155 153 IF( ln_ldfeiv_dia ) zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 156 154 & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * zslope_skew … … 166 164 DO jj = 1, jpjm1 167 165 DO ji = 1, fs_jpim1 168 ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp)169 zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk)166 ze3wr = 1.0_wp / e3w_n(ji,jj+jp,jk+kp) 167 zbv = e1e2v(ji,jj) * e3v_n(ji,jj,jk) 170 168 zah = 0.25_wp * pahv(ji,jj,jk) 171 169 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 172 170 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 173 171 ! (do this by *adding* gradient of depth) 174 zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)172 zslope2 = zslope_skew + ( gdept_n(ji,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) 175 173 zslope2 = zslope2 * zslope2 176 174 ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2 … … 193 191 DO ji = 1, fs_jpim1 194 192 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 195 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) ) )193 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) 196 194 END DO 197 195 END DO … … 201 199 DO jj = 1, jpjm1 202 200 DO ji = 1, fs_jpim1 203 ze3w_2 = fse3w(ji,jj,jk) * fse3w(ji,jj,jk)201 ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 204 202 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 205 203 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt … … 250 248 ENDIF 251 249 ENDIF 252 250 ! 253 251 !!---------------------------------------------------------------------- 254 252 !! II - horizontal trend (full) … … 256 254 ! 257 255 DO jk = 1, jpkm1 258 !259 256 ! !== Vertical tracer gradient at level jk and jk+1 260 257 zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) … … 274 271 ze1ur = r1_e1u(ji,jj) 275 272 zdxt = zdit(ji,jj,jk) * ze1ur 276 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)273 ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 277 274 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 278 275 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 279 276 zslope_iso = triadi (ji+ip,jj,jk,1-ip,kp) 280 281 zbu = 0.25_wp * e1e2u(ji,jj) * fse3u(ji,jj,jk)277 ! 278 zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk) 282 279 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked.... 283 280 zah = pahu(ji,jj,jk) … … 290 287 END DO 291 288 END DO 292 289 ! 293 290 DO jp = 0, 1 294 291 DO kp = 0, 1 … … 297 294 ze2vr = r1_e2v(ji,jj) 298 295 zdyt = zdjt(ji,jj,jk) * ze2vr 299 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)296 ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp) 300 297 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 301 298 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 302 299 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 303 zbv = 0.25_wp * e1e2v(ji,jj) * fse3v(ji,jj,jk)300 zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk) 304 301 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahv is masked... 305 302 zah = pahv(ji,jj,jk) … … 312 309 END DO 313 310 END DO 314 311 ! 315 312 ELSE 316 313 ! 317 314 DO ip = 0, 1 !== Horizontal & vertical fluxes 318 315 DO kp = 0, 1 … … 321 318 ze1ur = r1_e1u(ji,jj) 322 319 zdxt = zdit(ji,jj,jk) * ze1ur 323 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp)320 ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 324 321 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 325 322 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 326 323 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 327 328 zbu = 0.25_wp * e1e2u(ji,jj) * fse3u(ji,jj,jk)324 ! 325 zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk) 329 326 ! ln_botmix_triad is .F. mask zah for bottom half cells 330 327 zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp) ! pahu(ji+ip,jj,jk) ===>> ???? 331 328 zah_slp = zah * zslope_iso 332 IF( ln_ldfeiv ) zaei_slp = aeiu(ji,jj,jk) * zslope_skew ! fsaeit(ji+ip,jj,jk)*zslope_skew329 IF( ln_ldfeiv ) zaei_slp = aeiu(ji,jj,jk) * zslope_skew ! aeit(ji+ip,jj,jk)*zslope_skew 333 330 zftu(ji ,jj,jk ) = zftu(ji ,jj,jk ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 334 331 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr … … 337 334 END DO 338 335 END DO 339 336 ! 340 337 DO jp = 0, 1 341 338 DO kp = 0, 1 … … 344 341 ze2vr = r1_e2v(ji,jj) 345 342 zdyt = zdjt(ji,jj,jk) * ze2vr 346 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp)343 ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp) 347 344 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 348 345 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 349 346 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 350 zbv = 0.25_wp * e1e2v(ji,jj) * fse3v(ji,jj,jk)347 zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk) 351 348 ! ln_botmix_triad is .F. mask zah for bottom half cells 352 349 zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! pahv(ji,jj+jp,jk) ???? 353 350 zah_slp = zah * zslope_iso 354 IF( ln_ldfeiv ) zaei_slp = aeiv(ji,jj,jk) * zslope_skew ! fsaeit(ji,jj+jp,jk)*zslope_skew351 IF( ln_ldfeiv ) zaei_slp = aeiv(ji,jj,jk) * zslope_skew ! aeit(ji,jj+jp,jk)*zslope_skew 355 352 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 356 353 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr … … 365 362 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) & 366 363 & + zftv(ji,jj-1,jk) - zftv(ji,jj,jk) ) & 367 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )364 & / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 368 365 END DO 369 366 END DO … … 376 373 DO jj = 1, jpjm1 377 374 DO ji = fs_2, fs_jpim1 378 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) &375 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) & 379 376 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 380 377 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) … … 388 385 DO jj = 1, jpjm1 389 386 DO ji = fs_2, fs_jpim1 390 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) &387 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) & 391 388 & * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 392 389 END DO … … 397 394 DO jj = 1, jpjm1 398 395 DO ji = fs_2, fs_jpim1 399 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / fse3w(ji,jj,jk) * tmask(ji,jj,jk) &396 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) & 400 397 & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & 401 398 & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) … … 410 407 DO ji = fs_2, fs_jpim1 ! vector opt. 411 408 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) & 412 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )409 & / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 413 410 END DO 414 411 END DO
Note: See TracChangeset
for help on using the changeset viewer.
