Changeset 12377 for NEMO/trunk/src/OCE/LDF
- Timestamp:
- 2020-02-12T15:39:06+01:00 (4 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 5 edited
-
. (modified) (1 prop)
-
src/OCE/LDF/ldfc1d_c2d.F90 (modified) (3 diffs)
-
src/OCE/LDF/ldfdyn.F90 (modified) (9 diffs)
-
src/OCE/LDF/ldfslp.F90 (modified) (23 diffs)
-
src/OCE/LDF/ldftra.F90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEAD ext/AGRIF5 ^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL
-
- Property svn:externals
-
NEMO/trunk/src/OCE/LDF/ldfc1d_c2d.F90
r10425 r12377 30 30 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 31 31 32 !! * Substitutions 33 # include "do_loop_substitute.h90" 32 34 !!---------------------------------------------------------------------- 33 35 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 78 80 pah1(:,:,jk) = pahs1(:,:) * ( zratio + zc * ( 1._wp + TANH( - ( gdept_0(:,:,jk) - zh ) * zw) ) ) 79 81 END DO 80 DO jk = jpkm1, 1, -1 ! pah2 at F-point (zdep2 is an approximation in zps-coord.) 81 DO jj = 1, jpjm1 82 DO ji = 1, jpim1 83 zdep2 = ( gdept_0(ji,jj+1,jk) + gdept_0(ji+1,jj+1,jk) & 84 & + gdept_0(ji,jj ,jk) + gdept_0(ji+1,jj ,jk) ) * r1_4 85 pah2(ji,jj,jk) = pahs2(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) 86 END DO 87 END DO 88 END DO 82 DO_3DS_10_10( jpkm1, 1, -1 ) 83 zdep2 = ( gdept_0(ji,jj+1,jk) + gdept_0(ji+1,jj+1,jk) & 84 & + gdept_0(ji,jj ,jk) + gdept_0(ji+1,jj ,jk) ) * r1_4 85 pah2(ji,jj,jk) = pahs2(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) 86 END_3D 89 87 CALL lbc_lnk( 'ldfc1d_c2d', pah2, 'F', 1. ) ! Lateral boundary conditions 90 88 ! 91 89 CASE( 'TRA' ) ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.) 92 DO jk = jpkm1, 1, -1 93 DO jj = 1, jpjm1 94 DO ji = 1, jpim1 95 zdep1 = ( gdept_0(ji,jj,jk) + gdept_0(ji+1,jj,jk) ) * 0.5_wp 96 zdep2 = ( gdept_0(ji,jj,jk) + gdept_0(ji,jj+1,jk) ) * 0.5_wp 97 pah1(ji,jj,jk) = pahs1(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) ) ) 98 pah2(ji,jj,jk) = pahs2(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) 99 END DO 100 END DO 101 END DO 90 DO_3DS_10_10( jpkm1, 1, -1 ) 91 zdep1 = ( gdept_0(ji,jj,jk) + gdept_0(ji+1,jj,jk) ) * 0.5_wp 92 zdep2 = ( gdept_0(ji,jj,jk) + gdept_0(ji,jj+1,jk) ) * 0.5_wp 93 pah1(ji,jj,jk) = pahs1(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) ) ) 94 pah2(ji,jj,jk) = pahs2(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) 95 END_3D 102 96 ! Lateral boundary conditions 103 97 CALL lbc_lnk_multi( 'ldfc1d_c2d', pah1, 'U', 1. , pah2, 'V', 1. ) … … 141 135 ! 142 136 CASE( 'DYN' ) ! T- and F-points 143 DO jj = 1, jpj 144 DO ji = 1, jpi 145 pah1(ji,jj,1) = pUfac * MAX( e1t(ji,jj) , e2t(ji,jj) )**knn 146 pah2(ji,jj,1) = pUfac * MAX( e1f(ji,jj) , e2f(ji,jj) )**knn 147 END DO 148 END DO 137 DO_2D_11_11 138 pah1(ji,jj,1) = pUfac * MAX( e1t(ji,jj) , e2t(ji,jj) )**knn 139 pah2(ji,jj,1) = pUfac * MAX( e1f(ji,jj) , e2f(ji,jj) )**knn 140 END_2D 149 141 CASE( 'TRA' ) ! U- and V-points 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 pah1(ji,jj,1) = pUfac * MAX( e1u(ji,jj), e2u(ji,jj) )**knn 153 pah2(ji,jj,1) = pUfac * MAX( e1v(ji,jj), e2v(ji,jj) )**knn 154 END DO 155 END DO 142 DO_2D_11_11 143 pah1(ji,jj,1) = pUfac * MAX( e1u(ji,jj), e2u(ji,jj) )**knn 144 pah2(ji,jj,1) = pUfac * MAX( e1v(ji,jj), e2v(ji,jj) )**knn 145 END_2D 156 146 CASE DEFAULT ! error 157 147 CALL ctl_stop( 'ldf_c2d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' ) -
NEMO/trunk/src/OCE/LDF/ldfdyn.F90
r12276 r12377 73 73 74 74 !! * Substitutions 75 # include " vectopt_loop_substitute.h90"75 # include "do_loop_substitute.h90" 76 76 !!---------------------------------------------------------------------- 77 77 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 115 115 !!---------------------------------------------------------------------- 116 116 ! 117 REWIND( numnam_ref )118 117 READ ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901) 119 118 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in reference namelist' ) 120 119 121 REWIND( numnam_cfg )122 120 READ ( numnam_cfg, namdyn_ldf, IOSTAT = ios, ERR = 902 ) 123 121 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_ldf in configuration namelist' ) … … 313 311 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays') 314 312 ! 315 DO jj = 1, jpj ! Set local gridscale values 316 DO ji = 1, jpi 317 esqt(ji,jj) = ( 2._wp * e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 318 esqf(ji,jj) = ( 2._wp * e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2 319 END DO 320 END DO 313 DO_2D_11_11 314 esqt(ji,jj) = ( 2._wp * e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 315 esqf(ji,jj) = ( 2._wp * e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2 316 END_2D 321 317 ! 322 318 CASE DEFAULT … … 339 335 340 336 341 SUBROUTINE ldf_dyn( kt )337 SUBROUTINE ldf_dyn( kt, Kbb ) 342 338 !!---------------------------------------------------------------------- 343 339 !! *** ROUTINE ldf_dyn *** … … 357 353 !!---------------------------------------------------------------------- 358 354 INTEGER, INTENT(in) :: kt ! time step index 355 INTEGER, INTENT(in) :: Kbb ! ocean time level indices 359 356 ! 360 357 INTEGER :: ji, jj, jk ! dummy loop indices … … 371 368 IF( ln_dynldf_lap ) THEN ! laplacian operator : |u| e /12 = |u/144| e 372 369 DO jk = 1, jpkm1 373 DO jj = 2, jpjm1 374 DO ji = fs_2, fs_jpim1 375 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 376 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 377 zemax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 378 ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax * tmask(ji,jj,jk) ! 288= 12*12 * 2 379 END DO 380 END DO 381 DO jj = 1, jpjm1 382 DO ji = 1, fs_jpim1 383 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 384 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 385 zemax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 386 ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax * fmask(ji,jj,jk) ! 288= 12*12 * 2 387 END DO 388 END DO 370 DO_2D_00_00 371 zu2pv2_ij = uu(ji ,jj ,jk,Kbb) * uu(ji ,jj ,jk,Kbb) + vv(ji ,jj ,jk,Kbb) * vv(ji ,jj ,jk,Kbb) 372 zu2pv2_ij_m1 = uu(ji-1,jj ,jk,Kbb) * uu(ji-1,jj ,jk,Kbb) + vv(ji ,jj-1,jk,Kbb) * vv(ji ,jj-1,jk,Kbb) 373 zemax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 374 ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax * tmask(ji,jj,jk) ! 288= 12*12 * 2 375 END_2D 376 DO_2D_10_10 377 zu2pv2_ij_p1 = uu(ji ,jj+1,jk, Kbb) * uu(ji ,jj+1,jk, Kbb) + vv(ji+1,jj ,jk, Kbb) * vv(ji+1,jj ,jk, Kbb) 378 zu2pv2_ij = uu(ji ,jj ,jk, Kbb) * uu(ji ,jj ,jk, Kbb) + vv(ji ,jj ,jk, Kbb) * vv(ji ,jj ,jk, Kbb) 379 zemax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 380 ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax * fmask(ji,jj,jk) ! 288= 12*12 * 2 381 END_2D 389 382 END DO 390 383 ELSEIF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e 391 384 DO jk = 1, jpkm1 392 DO jj = 2, jpjm1 393 DO ji = fs_2, fs_jpim1 394 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 395 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 396 zemax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 397 ahmt(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax ) * zemax * tmask(ji,jj,jk) 398 END DO 399 END DO 400 DO jj = 1, jpjm1 401 DO ji = 1, fs_jpim1 402 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 403 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 404 zemax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 405 ahmf(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax ) * zemax * fmask(ji,jj,jk) 406 END DO 407 END DO 385 DO_2D_00_00 386 zu2pv2_ij = uu(ji ,jj ,jk,Kbb) * uu(ji ,jj ,jk,Kbb) + vv(ji ,jj ,jk,Kbb) * vv(ji ,jj ,jk,Kbb) 387 zu2pv2_ij_m1 = uu(ji-1,jj ,jk,Kbb) * uu(ji-1,jj ,jk,Kbb) + vv(ji ,jj-1,jk,Kbb) * vv(ji ,jj-1,jk,Kbb) 388 zemax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 389 ahmt(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zemax ) * zemax * tmask(ji,jj,jk) 390 END_2D 391 DO_2D_10_10 392 zu2pv2_ij_p1 = uu(ji ,jj+1,jk, Kbb) * uu(ji ,jj+1,jk, Kbb) + vv(ji+1,jj ,jk, Kbb) * vv(ji+1,jj ,jk, Kbb) 393 zu2pv2_ij = uu(ji ,jj ,jk, Kbb) * uu(ji ,jj ,jk, Kbb) + vv(ji ,jj ,jk, Kbb) * vv(ji ,jj ,jk, Kbb) 394 zemax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 395 ahmf(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zemax ) * zemax * fmask(ji,jj,jk) 396 END_2D 408 397 END DO 409 398 ENDIF … … 423 412 DO jk = 1, jpkm1 424 413 ! 425 DO jj = 2, jpjm1426 DO ji = 2, jpim1427 zdb = ( ub(ji,jj,jk) * r1_e2u(ji,jj) - ub(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) * r1_e1t(ji,jj) * e2t(ji,jj)&428 & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) * r1_e2t(ji,jj) * e1t(ji,jj)429 dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk)430 END DO431 END DO414 DO_2D_00_00 415 zdb = ( uu(ji,jj,jk,Kbb) * r1_e2u(ji,jj) - uu(ji-1,jj,jk,Kbb) * r1_e2u(ji-1,jj) ) & 416 & * r1_e1t(ji,jj) * e2t(ji,jj) & 417 & - ( vv(ji,jj,jk,Kbb) * r1_e1v(ji,jj) - vv(ji,jj-1,jk,Kbb) * r1_e1v(ji,jj-1) ) & 418 & * r1_e2t(ji,jj) * e1t(ji,jj) 419 dtensq(ji,jj,jk) = zdb * zdb * tmask(ji,jj,jk) 420 END_2D 432 421 ! 433 DO jj = 1, jpjm1434 DO ji = 1, jpim1435 zdb = ( ub(ji,jj+1,jk) * r1_e1u(ji,jj+1) - ub(ji,jj,jk) * r1_e1u(ji,jj) ) * r1_e2f(ji,jj) * e1f(ji,jj)&436 & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) * r1_e1f(ji,jj) * e2f(ji,jj)437 dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk)438 END DO439 END DO422 DO_2D_10_10 423 zdb = ( uu(ji,jj+1,jk,Kbb) * r1_e1u(ji,jj+1) - uu(ji,jj,jk,Kbb) * r1_e1u(ji,jj) ) & 424 & * r1_e2f(ji,jj) * e1f(ji,jj) & 425 & + ( vv(ji+1,jj,jk,Kbb) * r1_e2v(ji+1,jj) - vv(ji,jj,jk,Kbb) * r1_e2v(ji,jj) ) & 426 & * r1_e1f(ji,jj) * e2f(ji,jj) 427 dshesq(ji,jj,jk) = zdb * zdb * fmask(ji,jj,jk) 428 END_2D 440 429 ! 441 430 END DO … … 445 434 DO jk = 1, jpkm1 446 435 ! 447 DO jj = 2, jpjm1 ! T-point value 448 DO ji = fs_2, fs_jpim1 449 ! 450 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 451 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 452 ! 453 zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2 454 ahmt(ji,jj,jk) = zdelta * SQRT( dtensq(ji ,jj,jk) + & 455 & r1_4 * ( dshesq(ji ,jj,jk) + dshesq(ji ,jj-1,jk) + & 456 & dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) ) 457 ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 458 ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 459 ! 460 END DO 461 END DO 436 DO_2D_00_00 437 ! 438 zu2pv2_ij = uu(ji ,jj ,jk,Kbb) * uu(ji ,jj ,jk,Kbb) + vv(ji ,jj ,jk,Kbb) * vv(ji ,jj ,jk,Kbb) 439 zu2pv2_ij_m1 = uu(ji-1,jj ,jk,Kbb) * uu(ji-1,jj ,jk,Kbb) + vv(ji ,jj-1,jk,Kbb) * vv(ji ,jj-1,jk,Kbb) 440 ! 441 zdelta = zcmsmag * esqt(ji,jj) ! L^2 * (C_smag/pi)^2 442 ahmt(ji,jj,jk) = zdelta * SQRT( dtensq(ji ,jj,jk) + & 443 & r1_4 * ( dshesq(ji ,jj,jk) + dshesq(ji ,jj-1,jk) + & 444 & dshesq(ji-1,jj,jk) + dshesq(ji-1,jj-1,jk) ) ) 445 ahmt(ji,jj,jk) = MAX( ahmt(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 446 ahmt(ji,jj,jk) = MIN( ahmt(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 447 ! 448 END_2D 462 449 ! 463 DO jj = 1, jpjm1 ! F-point value 464 DO ji = 1, fs_jpim1 465 ! 466 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 467 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 468 ! 469 zdelta = zcmsmag * esqf(ji,jj) ! L^2 * (C_smag/pi)^2 470 ahmf(ji,jj,jk) = zdelta * SQRT( dshesq(ji ,jj,jk) + & 471 & r1_4 * ( dtensq(ji ,jj,jk) + dtensq(ji ,jj+1,jk) + & 472 & dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) ) 473 ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 474 ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 475 ! 476 END DO 477 END DO 450 DO_2D_10_10 451 ! 452 zu2pv2_ij_p1 = uu(ji ,jj+1,jk, kbb) * uu(ji ,jj+1,jk, kbb) + vv(ji+1,jj ,jk, kbb) * vv(ji+1,jj ,jk, kbb) 453 zu2pv2_ij = uu(ji ,jj ,jk, kbb) * uu(ji ,jj ,jk, kbb) + vv(ji ,jj ,jk, kbb) * vv(ji ,jj ,jk, kbb) 454 ! 455 zdelta = zcmsmag * esqf(ji,jj) ! L^2 * (C_smag/pi)^2 456 ahmf(ji,jj,jk) = zdelta * SQRT( dshesq(ji ,jj,jk) + & 457 & r1_4 * ( dtensq(ji ,jj,jk) + dtensq(ji ,jj+1,jk) + & 458 & dtensq(ji+1,jj,jk) + dtensq(ji+1,jj+1,jk) ) ) 459 ahmf(ji,jj,jk) = MAX( ahmf(ji,jj,jk), SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * zdelta * zstabf_lo ) ) ! Impose lower limit == minfac * |U|L/2 460 ahmf(ji,jj,jk) = MIN( ahmf(ji,jj,jk), zdelta * zstabf_up ) ! Impose upper limit == maxfac * L^2/(4*2dt) 461 ! 462 END_2D 478 463 ! 479 464 END DO … … 486 471 ! ! effective default limits are 1/12 |U|L^3 < B_hm < 1//(32*2dt) L^4 487 472 DO jk = 1, jpkm1 488 DO jj = 2, jpjm1 489 DO ji = fs_2, fs_jpim1 490 ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 491 END DO 492 END DO 493 DO jj = 1, jpjm1 494 DO ji = 1, fs_jpim1 495 ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 496 END DO 497 END DO 473 DO_2D_00_00 474 ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 475 END_2D 476 DO_2D_10_10 477 ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 478 END_2D 498 479 END DO 499 480 ! -
NEMO/trunk/src/OCE/LDF/ldfslp.F90
r12198 r12377 21 21 !!---------------------------------------------------------------------- 22 22 USE oce ! ocean dynamics and tracers 23 USE isf_oce ! ice shelf 23 24 USE dom_oce ! ocean space and time domain 24 25 ! USE ldfdyn ! lateral diffusion: eddy viscosity coef. … … 73 74 74 75 !! * Substitutions 75 # include " vectopt_loop_substitute.h90"76 # include "do_loop_substitute.h90" 76 77 !!---------------------------------------------------------------------- 77 78 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 81 82 CONTAINS 82 83 83 SUBROUTINE ldf_slp( kt, prd, pn2 )84 SUBROUTINE ldf_slp( kt, prd, pn2, Kbb, Kmm ) 84 85 !!---------------------------------------------------------------------- 85 86 !! *** ROUTINE ldf_slp *** … … 107 108 !!---------------------------------------------------------------------- 108 109 INTEGER , INTENT(in) :: kt ! ocean time-step index 110 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 109 111 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: prd ! in situ density 110 112 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: pn2 ! Brunt-Vaisala frequency (locally ref.) … … 134 136 zwz(:,:,:) = 0._wp 135 137 ! 136 DO jk = 1, jpk !== i- & j-gradient of density ==! 137 DO jj = 1, jpjm1 138 DO ji = 1, fs_jpim1 ! vector opt. 139 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 140 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 141 END DO 142 END DO 143 END DO 138 DO_3D_10_10( 1, jpk ) 139 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 140 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 141 END_3D 144 142 IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level 145 DO jj = 1, jpjm1 146 DO ji = 1, jpim1 147 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 148 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 149 END DO 150 END DO 143 DO_2D_10_10 144 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 145 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 146 END_2D 151 147 ENDIF 152 148 IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the bottom ocean level 153 DO jj = 1, jpjm1 154 DO ji = 1, jpim1 155 IF( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 156 IF( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 157 END DO 158 END DO 149 DO_2D_10_10 150 IF( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 151 IF( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 152 END_2D 159 153 ENDIF 160 154 ! … … 171 165 ! 172 166 ! !== Slopes just below the mixed layer ==! 173 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml167 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr, Kmm ) ! output: uslpml, vslpml, wslpiml, wslpjml 174 168 175 169 … … 178 172 ! 179 173 IF ( ln_isfcav ) THEN 180 DO jj = 2, jpjm1 181 DO ji = fs_2, fs_jpim1 ! vector opt. 182 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt (ji,jj), hmlpt (ji+1,jj ), 5._wp) & 183 & - MAX(risfdep(ji,jj), risfdep(ji+1,jj ) ) ) 184 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt (ji,jj), hmlpt (ji ,jj+1), 5._wp) & 185 & - MAX(risfdep(ji,jj), risfdep(ji ,jj+1) ) ) 186 END DO 187 END DO 174 DO_2D_00_00 175 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt (ji,jj), hmlpt (ji+1,jj ), 5._wp) & 176 & - MAX(risfdep(ji,jj), risfdep(ji+1,jj ) ) ) 177 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt (ji,jj), hmlpt (ji ,jj+1), 5._wp) & 178 & - MAX(risfdep(ji,jj), risfdep(ji ,jj+1) ) ) 179 END_2D 188 180 ELSE 189 DO jj = 2, jpjm1 190 DO ji = fs_2, fs_jpim1 ! vector opt. 191 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) 192 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) 193 END DO 194 END DO 181 DO_2D_00_00 182 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) 183 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) 184 END_2D 195 185 END IF 196 186 197 DO jk = 2, jpkm1 !* Slopes at u and v points 198 DO jj = 2, jpjm1 199 DO ji = fs_2, fs_jpim1 ! vector opt. 200 ! ! horizontal and vertical density gradient at u- and v-points 201 zau = zgru(ji,jj,jk) * r1_e1u(ji,jj) 202 zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj) 203 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 204 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 205 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 206 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 207 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,jk)* ABS( zau ) ) 208 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,jk)* ABS( zav ) ) 209 ! ! Fred Dupont: add a correction for bottom partial steps: 210 ! ! max slope = 1/2 * e3 / e1 211 IF (ln_zps .AND. jk==mbku(ji,jj)) & 212 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , - 2._wp * e1u(ji,jj) / e3u_n(ji,jj,jk)* ABS( zau ) ) 213 IF (ln_zps .AND. jk==mbkv(ji,jj)) & 214 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , - 2._wp * e2v(ji,jj) / e3v_n(ji,jj,jk)* ABS( zav ) ) 215 ! ! uslp and vslp output in zwz and zww, resp. 216 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 217 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 218 ! thickness of water column between surface and level k at u/v point 219 zdepu = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji+1,jj,jk) ) & 220 - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj) ) - e3u_n(ji,jj,miku(ji,jj)) ) 221 zdepv = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji,jj+1,jk) ) & 222 - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) - e3v_n(ji,jj,mikv(ji,jj)) ) 223 ! 224 zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps ) & 225 & + zfi * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk) 226 zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps ) & 227 & + zfj * zdepv * zslpml_hmlpv(ji,jj) ) * vmask(ji,jj,jk) 187 DO_3D_00_00( 2, jpkm1 ) 188 ! ! horizontal and vertical density gradient at u- and v-points 189 zau = zgru(ji,jj,jk) * r1_e1u(ji,jj) 190 zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj) 191 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 192 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 193 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 194 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 195 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u(ji,jj,jk,Kmm)* ABS( zau ) ) 196 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v(ji,jj,jk,Kmm)* ABS( zav ) ) 197 ! ! Fred Dupont: add a correction for bottom partial steps: 198 ! ! max slope = 1/2 * e3 / e1 199 IF (ln_zps .AND. jk==mbku(ji,jj)) & 200 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , - 2._wp * e1u(ji,jj) / e3u(ji,jj,jk,Kmm)* ABS( zau ) ) 201 IF (ln_zps .AND. jk==mbkv(ji,jj)) & 202 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , - 2._wp * e2v(ji,jj) / e3v(ji,jj,jk,Kmm)* ABS( zav ) ) 203 ! ! uslp and vslp output in zwz and zww, resp. 204 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 205 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 206 ! thickness of water column between surface and level k at u/v point 207 zdepu = 0.5_wp * ( ( gdept (ji,jj,jk,Kmm) + gdept (ji+1,jj,jk,Kmm) ) & 208 - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj) ) - e3u(ji,jj,miku(ji,jj),Kmm) ) 209 zdepv = 0.5_wp * ( ( gdept (ji,jj,jk,Kmm) + gdept (ji,jj+1,jk,Kmm) ) & 210 - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) - e3v(ji,jj,mikv(ji,jj),Kmm) ) 211 ! 212 zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps ) & 213 & + zfi * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk) 214 zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps ) & 215 & + zfj * zdepv * zslpml_hmlpv(ji,jj) ) * vmask(ji,jj,jk) 228 216 !!gm modif to suppress omlmask.... (as in Griffies case) 229 217 ! ! ! jk must be >= ML level for zf=1. otherwise zf=0. 230 218 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 231 219 ! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 232 ! zci = 0.5 * ( gdept _n(ji+1,jj,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )233 ! zcj = 0.5 * ( gdept _n(ji,jj+1,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )220 ! zci = 0.5 * ( gdept(ji+1,jj,jk,Kmm)+gdept(ji,jj,jk,Kmm) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 221 ! zcj = 0.5 * ( gdept(ji,jj+1,jk,Kmm)+gdept(ji,jj,jk,Kmm) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 234 222 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 235 223 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 236 224 !!gm end modif 237 END DO 238 END DO 239 END DO 225 END_3D 240 226 CALL lbc_lnk_multi( 'ldfslp', zwz, 'U', -1., zww, 'V', -1. ) ! lateral boundary conditions 241 227 ! 242 228 ! !* horizontal Shapiro filter 243 229 DO jk = 2, jpkm1 244 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 245 DO ji = 2, jpim1 246 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 247 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 248 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 249 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 250 & + 4.* zwz(ji ,jj ,jk) ) 251 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 252 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 253 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 254 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 255 & + 4.* zww(ji,jj ,jk) ) 256 END DO 257 END DO 230 DO_2D_00_00 231 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 232 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 233 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 234 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 235 & + 4.* zwz(ji ,jj ,jk) ) 236 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 237 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 238 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 239 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 240 & + 4.* zww(ji,jj ,jk) ) 241 END_2D 258 242 DO jj = 3, jpj-2 ! other rows 259 DO ji = fs_2, fs_jpim1 ! vector opt.243 DO ji = 2, jpim1 ! vector opt. 260 244 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 261 245 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & … … 271 255 END DO 272 256 ! !* decrease along coastal boundaries 273 DO jj = 2, jpjm1 274 DO ji = fs_2, fs_jpim1 ! vector opt. 275 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 276 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 277 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 278 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 279 END DO 280 END DO 257 DO_2D_00_00 258 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 259 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 260 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 261 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 262 END_2D 281 263 END DO 282 264 … … 285 267 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 286 268 ! 287 DO jk = 2, jpkm1 288 DO jj = 2, jpjm1 289 DO ji = fs_2, fs_jpim1 ! vector opt. 290 ! !* Local vertical density gradient evaluated from N^2 291 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 292 ! !* Slopes at w point 293 ! ! i- & j-gradient of density at w-points 294 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 295 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 296 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 297 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 298 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 299 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * wmask (ji,jj,jk) 300 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 301 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * wmask (ji,jj,jk) 302 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 303 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 304 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zai ) ) 305 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zaj ) ) 306 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 307 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 308 zck = ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj) - gdepw_n(ji,jj,mikt(ji,jj)), 10._wp ) 309 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 310 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk) 269 DO_3D_00_00( 2, jpkm1 ) 270 ! !* Local vertical density gradient evaluated from N^2 271 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 272 ! !* Slopes at w point 273 ! ! i- & j-gradient of density at w-points 274 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 275 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 276 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 277 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 278 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 279 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * wmask (ji,jj,jk) 280 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 281 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * wmask (ji,jj,jk) 282 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 283 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 284 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w(ji,jj,jk,Kmm)* ABS( zai ) ) 285 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w(ji,jj,jk,Kmm)* ABS( zaj ) ) 286 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 287 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 288 zck = ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) ) / MAX( hmlp(ji,jj) - gdepw(ji,jj,mikt(ji,jj),Kmm), 10._wp ) 289 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 290 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk) 311 291 312 292 !!gm modif to suppress omlmask.... (as in Griffies operator) … … 317 297 ! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 318 298 !!gm end modif 319 END DO 320 END DO 321 END DO 299 END_3D 322 300 CALL lbc_lnk_multi( 'ldfslp', zwz, 'T', -1., zww, 'T', -1. ) ! lateral boundary conditions 323 301 ! 324 302 ! !* horizontal Shapiro filter 325 303 DO jk = 2, jpkm1 326 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 327 DO ji = 2, jpim1 328 zcofw = wmask(ji,jj,jk) * z1_16 329 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 330 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 331 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 332 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 333 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 334 335 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 336 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 337 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 338 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 339 & + 4.* zww(ji ,jj ,jk) ) * zcofw 340 END DO 341 END DO 304 DO_2D_00_00 305 zcofw = wmask(ji,jj,jk) * z1_16 306 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 307 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 308 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 309 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 310 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 311 312 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 313 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 314 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 315 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 316 & + 4.* zww(ji ,jj ,jk) ) * zcofw 317 END_2D 342 318 DO jj = 3, jpj-2 ! other rows 343 DO ji = fs_2, fs_jpim1 ! vector opt.319 DO ji = 2, jpim1 ! vector opt. 344 320 zcofw = wmask(ji,jj,jk) * z1_16 345 321 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & … … 357 333 END DO 358 334 ! !* decrease in vicinity of topography 359 DO jj = 2, jpjm1 360 DO ji = fs_2, fs_jpim1 ! vector opt. 361 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 362 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 363 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 364 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 365 END DO 366 END DO 335 DO_2D_00_00 336 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 337 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 338 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 339 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 340 END_2D 367 341 END DO 368 342 … … 371 345 CALL lbc_lnk_multi( 'ldfslp', uslp , 'U', -1. , vslp , 'V', -1. , wslpi, 'W', -1., wslpj, 'W', -1. ) 372 346 373 IF( ln_ctl) THEN347 IF(sn_cfctl%l_prtctl) THEN 374 348 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 375 349 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) … … 381 355 382 356 383 SUBROUTINE ldf_slp_triad ( kt )357 SUBROUTINE ldf_slp_triad ( kt, Kbb, Kmm ) 384 358 !!---------------------------------------------------------------------- 385 359 !! *** ROUTINE ldf_slp_triad *** … … 396 370 !!---------------------------------------------------------------------- 397 371 INTEGER, INTENT( in ) :: kt ! ocean time-step index 372 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 398 373 !! 399 374 INTEGER :: ji, jj, jk, jl, ip, jp, kp ! dummy loop indices … … 421 396 ! 422 397 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln) 423 DO jk = 1, jpkm1 ! done each pair of triad 424 DO jj = 1, jpjm1 ! NB: not masked ==> a minimum value is set 425 DO ji = 1, fs_jpim1 ! vector opt. 426 zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! i-gradient of T & S at u-point 427 zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 428 zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! j-gradient of T & S at v-point 429 zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 430 zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj) 431 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj) 432 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 433 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 434 END DO 435 END DO 436 END DO 398 DO_3D_10_10( 1, jpkm1 ) 399 zdit = ( ts(ji+1,jj,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) ) ! i-gradient of T & S at u-point 400 zdis = ( ts(ji+1,jj,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) ) 401 zdjt = ( ts(ji,jj+1,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) ) ! j-gradient of T & S at v-point 402 zdjs = ( ts(ji,jj+1,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) ) 403 zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj) 404 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj) 405 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 406 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 407 END_3D 437 408 ! 438 409 IF( ln_zps .AND. l_grad_zps ) THEN ! partial steps: correction of i- & j-grad on bottom 439 DO jj = 1, jpjm1 440 DO ji = 1, jpim1 441 iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points) 442 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 443 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 444 zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj) 445 zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj) 446 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 447 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 448 END DO 449 END DO 410 DO_2D_10_10 411 iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points) 412 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 413 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 414 zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj) 415 zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj) 416 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 417 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 418 END_2D 450 419 ENDIF 451 420 ! … … 453 422 454 423 DO kp = 0, 1 !== unmasked before density i- j-, k-gradients ==! 455 DO jk = 1, jpkm1 ! done each pair of triad 456 DO jj = 1, jpj ! NB: not masked ==> a minimum value is set 457 DO ji = 1, jpi ! vector opt. 458 IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp 459 zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) ) 460 zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) ) 461 ELSE 462 zdkt = 0._wp ! 1st level gradient set to zero 463 zdks = 0._wp 464 ENDIF 465 zdzrho_raw = ( - rab_b(ji,jj,jk ,jp_tem) * zdkt & 466 & + rab_b(ji,jj,jk ,jp_sal) * zdks & 467 & ) / e3w_n(ji,jj,jk+kp) 468 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln 469 END DO 470 END DO 471 END DO 424 DO_3D_11_11( 1, jpkm1 ) 425 IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp 426 zdkt = ( ts(ji,jj,jk+kp-1,jp_tem,Kbb) - ts(ji,jj,jk+kp,jp_tem,Kbb) ) 427 zdks = ( ts(ji,jj,jk+kp-1,jp_sal,Kbb) - ts(ji,jj,jk+kp,jp_sal,Kbb) ) 428 ELSE 429 zdkt = 0._wp ! 1st level gradient set to zero 430 zdks = 0._wp 431 ENDIF 432 zdzrho_raw = ( - rab_b(ji,jj,jk ,jp_tem) * zdkt & 433 & + rab_b(ji,jj,jk ,jp_sal) * zdks & 434 & ) / e3w(ji,jj,jk+kp,Kmm) 435 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln 436 END_3D 472 437 END DO 473 438 ! 474 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==! 475 DO ji = 1, jpi 476 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth 477 z1_mlbw(ji,jj) = 1._wp / gdepw_n(ji,jj,jk) 478 END DO 479 END DO 439 DO_2D_11_11 440 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth 441 z1_mlbw(ji,jj) = 1._wp / gdepw(ji,jj,jk,Kmm) 442 END_2D 480 443 ! 481 444 ! !== intialisations to zero ==! … … 494 457 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 495 458 DO kp = 0, 1 ! with only the slope-max limit and MASKED 496 DO jj = 1, jpjm1 497 DO ji = 1, fs_jpim1 498 ip = jl ; jp = jl 499 ! 500 jk = nmln(ji+ip,jj) + 1 501 IF( jk > mbkt(ji+ip,jj) ) THEN ! ML reaches bottom 502 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 503 ELSE 504 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 505 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 506 & - ( gdept_n(ji+1,jj,jk-kp) - gdept_n(ji,jj,jk-kp) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk) 507 ze3_e1 = e3w_n(ji+ip,jj,jk-kp) * r1_e1u(ji,jj) 508 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1 , ABS( zti_g_raw ) ), zti_g_raw ) 509 ENDIF 510 ! 511 jk = nmln(ji,jj+jp) + 1 512 IF( jk > mbkt(ji,jj+jp) ) THEN !ML reaches bottom 513 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp 514 ELSE 515 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 516 & - ( gdept_n(ji,jj+1,jk-kp) - gdept_n(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 517 ze3_e2 = e3w_n(ji,jj+jp,jk-kp) / e2v(ji,jj) 518 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2 , ABS( ztj_g_raw ) ), ztj_g_raw ) 519 ENDIF 520 END DO 521 END DO 459 DO_2D_10_10 460 ip = jl ; jp = jl 461 ! 462 jk = nmln(ji+ip,jj) + 1 463 IF( jk > mbkt(ji+ip,jj) ) THEN ! ML reaches bottom 464 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 465 ELSE 466 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 467 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 468 & - ( gdept(ji+1,jj,jk-kp,Kmm) - gdept(ji,jj,jk-kp,Kmm) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk) 469 ze3_e1 = e3w(ji+ip,jj,jk-kp,Kmm) * r1_e1u(ji,jj) 470 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1 , ABS( zti_g_raw ) ), zti_g_raw ) 471 ENDIF 472 ! 473 jk = nmln(ji,jj+jp) + 1 474 IF( jk > mbkt(ji,jj+jp) ) THEN !ML reaches bottom 475 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp 476 ELSE 477 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 478 & - ( gdept(ji,jj+1,jk-kp,Kmm) - gdept(ji,jj,jk-kp,Kmm) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 479 ze3_e2 = e3w(ji,jj+jp,jk-kp,Kmm) / e2v(ji,jj) 480 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2 , ABS( ztj_g_raw ) ), ztj_g_raw ) 481 ENDIF 482 END_2D 522 483 END DO 523 484 END DO … … 533 494 ! Must mask contribution to slope from dz/dx at constant s for triads jk=1,kp=0 that poke up though ocean surface 534 495 znot_thru_surface = REAL( 1-1/(jk+kp), wp ) !jk+kp=1,=0.; otherwise=1.0 535 DO jj = 1, jpjm1 536 DO ji = 1, fs_jpim1 ! vector opt. 537 ! 538 ! Calculate slope relative to geopotentials used for GM skew fluxes 539 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 540 ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 541 ! masked by umask taken at the level of dz(rho) 542 ! 543 ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti) 544 ! 545 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 546 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 547 ! 548 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 549 zti_coord = znot_thru_surface * ( gdept_n(ji+1,jj ,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) 550 ztj_coord = znot_thru_surface * ( gdept_n(ji ,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) ! unmasked 551 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 552 ztj_g_raw = ztj_raw - ztj_coord 553 ! additional limit required in bilaplacian case 554 ze3_e1 = e3w_n(ji+ip,jj ,jk+kp) * r1_e1u(ji,jj) 555 ze3_e2 = e3w_n(ji ,jj+jp,jk+kp) * r1_e2v(ji,jj) 556 ! NB: hard coded factor 5 (can be a namelist parameter...) 557 zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) 558 ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw ) 559 ! 560 ! Below ML use limited zti_g as is & mask 561 ! Inside ML replace by linearly reducing sx_mlb towards surface & mask 562 ! 563 zfacti = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)), wp ) ! k index of uppermost point(s) of triad is jk+kp-1 564 zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1 565 ! ! otherwise zfact=0 566 zti_g_lim = ( zfacti * zti_g_lim & 567 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 568 & * gdepw_n(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 569 ztj_g_lim = ( zfactj * ztj_g_lim & 570 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 571 & * gdepw_n(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 572 ! 573 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim 574 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim 575 ! 576 ! Get coefficients of isoneutral diffusion tensor 577 ! 1. Utilise gradients *relative* to s-coordinate, so add t-point slopes (*subtract* depth gradients) 578 ! 2. We require that isoneutral diffusion gives no vertical buoyancy flux 579 ! i.e. 33 term = (real slope* 31, 13 terms) 580 ! To do this, retain limited sx**2 in vertical flux, but divide by real slope for 13/31 terms 581 ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 582 ! 583 zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces 584 ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 585 ! 586 IF( ln_triad_iso ) THEN 587 zti_raw = zti_lim*zti_lim / zti_raw 588 ztj_raw = ztj_lim*ztj_lim / ztj_raw 589 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 590 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 591 zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw 592 ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw 593 ENDIF 594 ! ! switching triad scheme 595 zisw = (1._wp - rn_sw_triad ) + rn_sw_triad & 596 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) ) ) 597 zjsw = (1._wp - rn_sw_triad ) + rn_sw_triad & 598 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) ) ) 599 ! 600 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim * zisw 601 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 602 ! 603 zbu = e1e2u(ji ,jj ) * e3u_n(ji ,jj ,jk ) 604 zbv = e1e2v(ji ,jj ) * e3v_n(ji ,jj ,jk ) 605 zbti = e1e2t(ji+ip,jj ) * e3w_n(ji+ip,jj ,jk+kp) 606 zbtj = e1e2t(ji ,jj+jp) * e3w_n(ji ,jj+jp,jk+kp) 607 ! 608 wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked 609 wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim 610 END DO 611 END DO 496 DO_2D_10_10 497 ! 498 ! Calculate slope relative to geopotentials used for GM skew fluxes 499 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 500 ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 501 ! masked by umask taken at the level of dz(rho) 502 ! 503 ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti) 504 ! 505 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 506 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 507 ! 508 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 509 zti_coord = znot_thru_surface * ( gdept(ji+1,jj ,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) 510 ztj_coord = znot_thru_surface * ( gdept(ji ,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) ! unmasked 511 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 512 ztj_g_raw = ztj_raw - ztj_coord 513 ! additional limit required in bilaplacian case 514 ze3_e1 = e3w(ji+ip,jj ,jk+kp,Kmm) * r1_e1u(ji,jj) 515 ze3_e2 = e3w(ji ,jj+jp,jk+kp,Kmm) * r1_e2v(ji,jj) 516 ! NB: hard coded factor 5 (can be a namelist parameter...) 517 zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) 518 ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw ) 519 ! 520 ! Below ML use limited zti_g as is & mask 521 ! Inside ML replace by linearly reducing sx_mlb towards surface & mask 522 ! 523 zfacti = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)), wp ) ! k index of uppermost point(s) of triad is jk+kp-1 524 zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1 525 ! ! otherwise zfact=0 526 zti_g_lim = ( zfacti * zti_g_lim & 527 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 528 & * gdepw(ji+ip,jj,jk+kp,Kmm) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 529 ztj_g_lim = ( zfactj * ztj_g_lim & 530 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 531 & * gdepw(ji,jj+jp,jk+kp,Kmm) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 532 ! 533 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim 534 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim 535 ! 536 ! Get coefficients of isoneutral diffusion tensor 537 ! 1. Utilise gradients *relative* to s-coordinate, so add t-point slopes (*subtract* depth gradients) 538 ! 2. We require that isoneutral diffusion gives no vertical buoyancy flux 539 ! i.e. 33 term = (real slope* 31, 13 terms) 540 ! To do this, retain limited sx**2 in vertical flux, but divide by real slope for 13/31 terms 541 ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 542 ! 543 zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces 544 ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 545 ! 546 IF( ln_triad_iso ) THEN 547 zti_raw = zti_lim*zti_lim / zti_raw 548 ztj_raw = ztj_lim*ztj_lim / ztj_raw 549 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 550 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 551 zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw 552 ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw 553 ENDIF 554 ! ! switching triad scheme 555 zisw = (1._wp - rn_sw_triad ) + rn_sw_triad & 556 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) ) ) 557 zjsw = (1._wp - rn_sw_triad ) + rn_sw_triad & 558 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) ) ) 559 ! 560 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim * zisw 561 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 562 ! 563 zbu = e1e2u(ji ,jj ) * e3u(ji ,jj ,jk ,Kmm) 564 zbv = e1e2v(ji ,jj ) * e3v(ji ,jj ,jk ,Kmm) 565 zbti = e1e2t(ji+ip,jj ) * e3w(ji+ip,jj ,jk+kp,Kmm) 566 zbtj = e1e2t(ji ,jj+jp) * e3w(ji ,jj+jp,jk+kp,Kmm) 567 ! 568 wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked 569 wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim 570 END_2D 612 571 END DO 613 572 END DO … … 623 582 624 583 625 SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr )584 SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr, Kmm ) 626 585 !!---------------------------------------------------------------------- 627 586 !! *** ROUTINE ldf_slp_mxl *** … … 643 602 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_gru, p_grv ! i- & j-gradient of density (u- & v-pts) 644 603 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_dzr ! z-gradient of density (T-point) 604 INTEGER , INTENT(in) :: Kmm ! ocean time level indices 645 605 !! 646 606 INTEGER :: ji , jj , jk ! dummy loop indices … … 663 623 ! 664 624 ! !== surface mixed layer mask ! 665 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise 666 DO jj = 1, jpj 667 DO ji = 1, jpi 668 ik = nmln(ji,jj) - 1 669 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp 670 ELSE ; omlmask(ji,jj,jk) = 0._wp 671 ENDIF 672 END DO 673 END DO 674 END DO 625 DO_3D_11_11( 1, jpk ) 626 ik = nmln(ji,jj) - 1 627 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp 628 ELSE ; omlmask(ji,jj,jk) = 0._wp 629 ENDIF 630 END_3D 675 631 676 632 … … 685 641 !----------------------------------------------------------------------- 686 642 ! 687 DO jj = 2, jpjm1 688 DO ji = 2, jpim1 689 ! !== Slope at u- & v-points just below the Mixed Layer ==! 690 ! 691 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 692 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 693 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 694 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 695 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 696 ! !- horizontal density gradient at u- & v-points 697 zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj) 698 zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj) 699 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 700 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 701 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,iku)* ABS( zau ) ) 702 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,ikv)* ABS( zav ) ) 703 ! !- Slope at u- & v-points (uslpml, vslpml) 704 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 705 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 706 ! 707 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! 708 ! 709 ik = MIN( nmln(ji,jj) + 1, jpk ) 710 ikm1 = MAX( 1, ik-1 ) 711 ! !- vertical density gradient for w-slope (from N^2) 712 zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 713 ! !- horizontal density i- & j-gradient at w-points 714 zci = MAX( umask(ji-1,jj,ik ) + umask(ji,jj,ik ) & 715 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 716 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 717 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) 718 zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) & 719 & + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1 ) ) / zci * tmask(ji,jj,ik) 720 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & 721 & + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik) 722 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 723 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 724 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zai ) ) 725 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zaj ) ) 726 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 727 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 728 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 729 END DO 730 END DO 643 DO_2D_00_00 644 ! !== Slope at u- & v-points just below the Mixed Layer ==! 645 ! 646 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 647 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 648 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 649 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 650 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 651 ! !- horizontal density gradient at u- & v-points 652 zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj) 653 zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj) 654 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 655 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 656 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u(ji,jj,iku,Kmm)* ABS( zau ) ) 657 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v(ji,jj,ikv,Kmm)* ABS( zav ) ) 658 ! !- Slope at u- & v-points (uslpml, vslpml) 659 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 660 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 661 ! 662 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! 663 ! 664 ik = MIN( nmln(ji,jj) + 1, jpk ) 665 ikm1 = MAX( 1, ik-1 ) 666 ! !- vertical density gradient for w-slope (from N^2) 667 zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 668 ! !- horizontal density i- & j-gradient at w-points 669 zci = MAX( umask(ji-1,jj,ik ) + umask(ji,jj,ik ) & 670 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 671 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 672 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) 673 zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) & 674 & + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1 ) ) / zci * tmask(ji,jj,ik) 675 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & 676 & + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik) 677 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 678 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 679 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w(ji,jj,ik,Kmm)* ABS( zai ) ) 680 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w(ji,jj,ik,Kmm)* ABS( zaj ) ) 681 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 682 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 683 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 684 END_2D 731 685 !!gm this lbc_lnk should be useless.... 732 686 CALL lbc_lnk_multi( 'ldfslp', uslpml , 'U', -1. , vslpml , 'V', -1. , wslpiml, 'W', -1. , wslpjml, 'W', -1. ) … … 790 744 ! DO jk = 1, jpk 791 745 ! DO jj = 2, jpjm1 792 ! DO ji = fs_2, fs_jpim1 ! vector opt.793 ! uslp (ji,jj,jk) = - ( gdept _n(ji+1,jj,jk) - gdept_n(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)794 ! vslp (ji,jj,jk) = - ( gdept _n(ji,jj+1,jk) - gdept_n(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)795 ! wslpi(ji,jj,jk) = - ( gdepw _n(ji+1,jj,jk) - gdepw_n(ji-1,jj,jk) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5796 ! wslpj(ji,jj,jk) = - ( gdepw _n(ji,jj+1,jk) - gdepw_n(ji,jj-1,jk) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5746 ! DO ji = 2, jpim1 ! vector opt. 747 ! uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kmm) - gdept(ji ,jj ,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 748 ! vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kmm) - gdept(ji ,jj ,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 749 ! wslpi(ji,jj,jk) = - ( gdepw(ji+1,jj,jk,Kmm) - gdepw(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5 750 ! wslpj(ji,jj,jk) = - ( gdepw(ji,jj+1,jk,Kmm) - gdepw(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5 797 751 ! END DO 798 752 ! END DO -
NEMO/trunk/src/OCE/LDF/ldftra.F90
r12296 r12377 94 94 95 95 !! * Substitutions 96 # include " vectopt_loop_substitute.h90"96 # include "do_loop_substitute.h90" 97 97 !!---------------------------------------------------------------------- 98 98 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 152 152 ! ================================= 153 153 ! 154 REWIND( numnam_ref )155 154 READ ( numnam_ref, namtra_ldf, IOSTAT = ios, ERR = 901) 156 155 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldf in reference namelist' ) 157 158 REWIND( numnam_cfg )159 156 READ ( numnam_cfg, namtra_ldf, IOSTAT = ios, ERR = 902 ) 160 157 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_ldf in configuration namelist' ) … … 383 380 384 381 385 SUBROUTINE ldf_tra( kt )382 SUBROUTINE ldf_tra( kt, Kbb, Kmm ) 386 383 !!---------------------------------------------------------------------- 387 384 !! *** ROUTINE ldf_tra *** … … 404 401 !!---------------------------------------------------------------------- 405 402 INTEGER, INTENT(in) :: kt ! time step 403 INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices 406 404 ! 407 405 INTEGER :: ji, jj, jk ! dummy loop indices … … 412 410 ! ! =F(growth rate of baroclinic instability) 413 411 ! ! max value aeiv_0 ; decreased to 0 within 20N-20S 414 CALL ldf_eiv( kt, aei0, aeiu, aeiv )412 CALL ldf_eiv( kt, aei0, aeiu, aeiv, Kmm ) 415 413 ENDIF 416 414 ! … … 425 423 ahtv(:,:,1) = aeiv(:,:,1) 426 424 ELSE ! compute aht. 427 CALL ldf_eiv( kt, aht0, ahtu, ahtv )425 CALL ldf_eiv( kt, aht0, ahtu, ahtv, Kmm ) 428 426 ENDIF 429 427 ! … … 431 429 zaht_min = 0.2_wp * aht0 ! minimum value for aht 432 430 zDaht = aht0 - zaht_min 433 DO jj = 1, jpj 434 DO ji = 1, jpi 435 !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg) 436 !! ==>>> The Coriolis value is identical for t- & u_points, and for v- and f-points 437 zaht = ( 1._wp - MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * zDaht 438 zahf = ( 1._wp - MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * zDaht 439 ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) ! min value zaht_min 440 ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zahf ) ! increase within 20S-20N 441 END DO 442 END DO 431 DO_2D_11_11 432 !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg) 433 !! ==>>> The Coriolis value is identical for t- & u_points, and for v- and f-points 434 zaht = ( 1._wp - MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * zDaht 435 zahf = ( 1._wp - MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * zDaht 436 ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) ! min value zaht_min 437 ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zahf ) ! increase within 20S-20N 438 END_2D 443 439 DO jk = 1, jpkm1 ! deeper value = surface value + mask for all levels 444 440 ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk) … … 449 445 IF( ln_traldf_lap ) THEN ! laplacian operator |u| e /12 450 446 DO jk = 1, jpkm1 451 ahtu(:,:,jk) = ABS( u b(:,:,jk) ) * e1u(:,:) * r1_12 ! n.b. ub,vbare masked452 ahtv(:,:,jk) = ABS( v b(:,:,jk) ) * e2v(:,:) * r1_12447 ahtu(:,:,jk) = ABS( uu(:,:,jk,Kbb) ) * e1u(:,:) * r1_12 ! n.b. uu,vv are masked 448 ahtv(:,:,jk) = ABS( vv(:,:,jk,Kbb) ) * e2v(:,:) * r1_12 453 449 END DO 454 450 ELSEIF( ln_traldf_blp ) THEN ! bilaplacian operator sqrt( |u| e^3 /12 ) = sqrt( |u| e /12 ) * e 455 451 DO jk = 1, jpkm1 456 ahtu(:,:,jk) = SQRT( ABS( u b(:,:,jk) ) * e1u(:,:) * r1_12 ) * e1u(:,:)457 ahtv(:,:,jk) = SQRT( ABS( v b(:,:,jk) ) * e2v(:,:) * r1_12 ) * e2v(:,:)452 ahtu(:,:,jk) = SQRT( ABS( uu(:,:,jk,Kbb) ) * e1u(:,:) * r1_12 ) * e1u(:,:) 453 ahtv(:,:,jk) = SQRT( ABS( vv(:,:,jk,Kbb) ) * e2v(:,:) * r1_12 ) * e2v(:,:) 458 454 END DO 459 455 ENDIF … … 511 507 ENDIF 512 508 ! 513 REWIND( numnam_ref )514 509 READ ( numnam_ref, namtra_eiv, IOSTAT = ios, ERR = 901) 515 510 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_eiv in reference namelist' ) 516 511 ! 517 REWIND( numnam_cfg )518 512 READ ( numnam_cfg, namtra_eiv, IOSTAT = ios, ERR = 902 ) 519 513 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_eiv in configuration namelist' ) … … 626 620 627 621 628 SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv )622 SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv, Kmm ) 629 623 !!---------------------------------------------------------------------- 630 624 !! *** ROUTINE ldf_eiv *** … … 638 632 !!---------------------------------------------------------------------- 639 633 INTEGER , INTENT(in ) :: kt ! ocean time-step index 634 INTEGER , INTENT(in ) :: Kmm ! ocean time level indices 640 635 REAL(wp) , INTENT(inout) :: paei0 ! max value [m2/s] 641 636 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: paeiu, paeiv ! eiv coefficient [m2/s] … … 652 647 ! ! Compute lateral diffusive coefficient at T-point 653 648 IF( ln_traldf_triad ) THEN 654 DO jk = 1, jpk 655 DO jj = 2, jpjm1 656 DO ji = 2, jpim1 657 ! Take the max of N^2 and zero then take the vertical sum 658 ! of the square root of the resulting N^2 ( required to compute 659 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 660 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 661 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk) 662 ! Compute elements required for the inverse time scale of baroclinic 663 ! eddies using the isopycnal slopes calculated in ldfslp.F : 664 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 665 ze3w = e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 666 zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w 667 zhw(ji,jj) = zhw(ji,jj) + ze3w 668 END DO 669 END DO 670 END DO 649 DO_3D_00_00( 1, jpk ) 650 ! Take the max of N^2 and zero then take the vertical sum 651 ! of the square root of the resulting N^2 ( required to compute 652 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 653 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 654 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm) 655 ! Compute elements required for the inverse time scale of baroclinic 656 ! eddies using the isopycnal slopes calculated in ldfslp.F : 657 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 658 ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 659 zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w 660 zhw(ji,jj) = zhw(ji,jj) + ze3w 661 END_3D 671 662 ELSE 672 DO jk = 1, jpk 673 DO jj = 2, jpjm1 674 DO ji = 2, jpim1 675 ! Take the max of N^2 and zero then take the vertical sum 676 ! of the square root of the resulting N^2 ( required to compute 677 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 678 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 679 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk) 680 ! Compute elements required for the inverse time scale of baroclinic 681 ! eddies using the isopycnal slopes calculated in ldfslp.F : 682 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 683 ze3w = e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 684 zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 685 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w 686 zhw(ji,jj) = zhw(ji,jj) + ze3w 687 END DO 688 END DO 689 END DO 690 ENDIF 691 692 DO jj = 2, jpjm1 693 DO ji = fs_2, fs_jpim1 ! vector opt. 694 zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 695 ! Rossby radius at w-point taken betwenn 2 km and 40km 696 zRo(ji,jj) = MAX( 2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 ) ) 697 ! Compute aeiw by multiplying Ro^2 and T^-1 698 zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 699 END DO 700 END DO 663 DO_3D_00_00( 1, jpk ) 664 ! Take the max of N^2 and zero then take the vertical sum 665 ! of the square root of the resulting N^2 ( required to compute 666 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 667 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 668 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm) 669 ! Compute elements required for the inverse time scale of baroclinic 670 ! eddies using the isopycnal slopes calculated in ldfslp.F : 671 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 672 ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 673 zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 674 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w 675 zhw(ji,jj) = zhw(ji,jj) + ze3w 676 END_3D 677 ENDIF 678 679 DO_2D_00_00 680 zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 681 ! Rossby radius at w-point taken betwenn 2 km and 40km 682 zRo(ji,jj) = MAX( 2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 ) ) 683 ! Compute aeiw by multiplying Ro^2 and T^-1 684 zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 685 END_2D 701 686 702 687 ! !== Bound on eiv coeff. ==! 703 688 z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) ) 704 DO jj = 2, jpjm1 705 DO ji = fs_2, fs_jpim1 ! vector opt. 706 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease 707 zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0 708 END DO 709 END DO 689 DO_2D_00_00 690 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease 691 zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0 692 END_2D 710 693 CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1. ) ! lateral boundary condition 711 694 ! 712 DO jj = 2, jpjm1 !== aei at u- and v-points ==! 713 DO ji = fs_2, fs_jpim1 ! vector opt. 714 paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj ) ) * umask(ji,jj,1) 715 paeiv(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji ,jj+1) ) * vmask(ji,jj,1) 716 END DO 717 END DO 695 DO_2D_00_00 696 paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj ) ) * umask(ji,jj,1) 697 paeiv(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji ,jj+1) ) * vmask(ji,jj,1) 698 END_2D 718 699 CALL lbc_lnk_multi( 'ldftra', paeiu(:,:,1), 'U', 1. , paeiv(:,:,1), 'V', 1. ) ! lateral boundary condition 719 700 … … 726 707 727 708 728 SUBROUTINE ldf_eiv_trp( kt, kit000, pu n, pvn, pwn, cdtype)709 SUBROUTINE ldf_eiv_trp( kt, kit000, pu, pv, pw, cdtype, Kmm, Krhs ) 729 710 !!---------------------------------------------------------------------- 730 711 !! *** ROUTINE ldf_eiv_trp *** … … 742 723 !! velocity and heat transport (call ldf_eiv_dia) 743 724 !! 744 !! ** Action : pun, pvn increased by the eiv transport 745 !!---------------------------------------------------------------------- 746 INTEGER , INTENT(in ) :: kt ! ocean time-step index 747 INTEGER , INTENT(in ) :: kit000 ! first time step index 748 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 749 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pun ! in : 3 ocean transport components [m3/s] 750 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pvn ! out: 3 ocean transport components [m3/s] 751 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pwn ! increased by the eiv [m3/s] 725 !! ** Action : pu, pv increased by the eiv transport 726 !!---------------------------------------------------------------------- 727 INTEGER , INTENT(in ) :: kt ! ocean time-step index 728 INTEGER , INTENT(in ) :: kit000 ! first time step index 729 INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices 730 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 731 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components [m3/s] 732 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pv ! out: 3 ocean transport components [m3/s] 733 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the eiv [m3/s] 752 734 !! 753 735 INTEGER :: ji, jj, jk ! dummy loop indices … … 767 749 zpsi_uw(:,:,jpk) = 0._wp ; zpsi_vw(:,:,jpk) = 0._wp 768 750 ! 769 DO jk = 2, jpkm1 770 DO jj = 1, jpjm1 771 DO ji = 1, fs_jpim1 ! vector opt. 772 zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk ) + wslpi(ji+1,jj,jk) ) & 773 & * ( aeiu (ji,jj,jk-1) + aeiu (ji ,jj,jk) ) * wumask(ji,jj,jk) 774 zpsi_vw(ji,jj,jk) = - r1_4 * e1v(ji,jj) * ( wslpj(ji,jj,jk ) + wslpj(ji,jj+1,jk) ) & 775 & * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj ,jk) ) * wvmask(ji,jj,jk) 776 END DO 777 END DO 778 END DO 779 ! 780 DO jk = 1, jpkm1 781 DO jj = 1, jpjm1 782 DO ji = 1, fs_jpim1 ! vector opt. 783 pun(ji,jj,jk) = pun(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 784 pvn(ji,jj,jk) = pvn(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 785 END DO 786 END DO 787 END DO 788 DO jk = 1, jpkm1 789 DO jj = 2, jpjm1 790 DO ji = fs_2, fs_jpim1 ! vector opt. 791 pwn(ji,jj,jk) = pwn(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj ,jk) & 792 & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji ,jj-1,jk) ) 793 END DO 794 END DO 795 END DO 751 DO_3D_10_10( 2, jpkm1 ) 752 zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk ) + wslpi(ji+1,jj,jk) ) & 753 & * ( aeiu (ji,jj,jk-1) + aeiu (ji ,jj,jk) ) * wumask(ji,jj,jk) 754 zpsi_vw(ji,jj,jk) = - r1_4 * e1v(ji,jj) * ( wslpj(ji,jj,jk ) + wslpj(ji,jj+1,jk) ) & 755 & * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj ,jk) ) * wvmask(ji,jj,jk) 756 END_3D 757 ! 758 DO_3D_10_10( 1, jpkm1 ) 759 pu(ji,jj,jk) = pu(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 760 pv(ji,jj,jk) = pv(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 761 END_3D 762 DO_3D_00_00( 1, jpkm1 ) 763 pw(ji,jj,jk) = pw(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj ,jk) & 764 & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji ,jj-1,jk) ) 765 END_3D 796 766 ! 797 767 ! ! diagnose the eddy induced velocity and associated heat transport 798 IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) CALL ldf_eiv_dia( zpsi_uw, zpsi_vw )768 IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 799 769 ! 800 770 END SUBROUTINE ldf_eiv_trp 801 771 802 772 803 SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw )773 SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw, Kmm ) 804 774 !!---------------------------------------------------------------------- 805 775 !! *** ROUTINE ldf_eiv_dia *** … … 812 782 !!---------------------------------------------------------------------- 813 783 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: psi_uw, psi_vw ! streamfunction [m3/s] 784 INTEGER , INTENT(in ) :: Kmm ! ocean time level indices 814 785 ! 815 786 INTEGER :: ji, jj, jk ! dummy loop indices … … 832 803 ! 833 804 DO jk = 1, jpkm1 ! e2u e3u u_eiv = -dk[psi_uw] 834 zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u _n(:,:,jk) )805 zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u(:,:,jk,Kmm) ) 835 806 END DO 836 807 CALL iom_put( "uoce_eiv", zw3d ) 837 808 ! 838 809 DO jk = 1, jpkm1 ! e1v e3v v_eiv = -dk[psi_vw] 839 zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v _n(:,:,jk) )810 zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v(:,:,jk,Kmm) ) 840 811 END DO 841 812 CALL iom_put( "voce_eiv", zw3d ) 842 813 ! 843 DO jk = 1, jpkm1 ! e1 e2 w_eiv = dk[psix] + dk[psix] 844 DO jj = 2, jpjm1 845 DO ji = fs_2, fs_jpim1 ! vector opt. 846 zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk) - psi_vw(ji ,jj-1,jk) & 847 & + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj ,jk) ) / e1e2t(ji,jj) 848 END DO 849 END DO 850 END DO 814 DO_3D_00_00( 1, jpkm1 ) 815 zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk) - psi_vw(ji ,jj-1,jk) & 816 & + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj ,jk) ) / e1e2t(ji,jj) 817 END_3D 851 818 CALL lbc_lnk( 'ldftra', zw3d, 'T', 1. ) ! lateral boundary condition 852 819 CALL iom_put( "woce_eiv", zw3d ) … … 872 839 zw2d(:,:) = 0._wp 873 840 zw3d(:,:,:) = 0._wp 874 DO jk = 1, jpkm1 875 DO jj = 2, jpjm1 876 DO ji = fs_2, fs_jpim1 ! vector opt. 877 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) & 878 & * ( tsn (ji,jj,jk,jp_tem) + tsn (ji+1,jj,jk,jp_tem) ) 879 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 880 END DO 881 END DO 882 END DO 841 DO_3D_00_00( 1, jpkm1 ) 842 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1) - psi_uw(ji ,jj,jk) ) & 843 & * ( ts (ji,jj,jk,jp_tem,Kmm) + ts (ji+1,jj,jk,jp_tem,Kmm) ) 844 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 845 END_3D 883 846 CALL lbc_lnk( 'ldftra', zw2d, 'U', -1. ) 884 847 CALL lbc_lnk( 'ldftra', zw3d, 'U', -1. ) … … 897 860 zw2d(:,:) = 0._wp 898 861 zw3d(:,:,:) = 0._wp 899 DO jk = 1, jpkm1 900 DO jj = 2, jpjm1 901 DO ji = fs_2, fs_jpim1 ! vector opt. 902 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) & 903 & * ( tsn (ji,jj,jk,jp_tem) + tsn (ji,jj+1,jk,jp_tem) ) 904 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 905 END DO 906 END DO 907 END DO 862 DO_3D_00_00( 1, jpkm1 ) 863 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj ,jk) ) & 864 & * ( ts (ji,jj,jk,jp_tem,Kmm) + ts (ji,jj+1,jk,jp_tem,Kmm) ) 865 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 866 END_3D 908 867 CALL lbc_lnk( 'ldftra', zw2d, 'V', -1. ) 909 868 CALL iom_put( "veiv_heattr", zztmp * zw2d ) ! heat transport in j-direction 910 869 CALL iom_put( "veiv_heattr", zztmp * zw3d ) ! heat transport in j-direction 911 870 ! 912 IF( ln_diaptr) CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )871 IF( iom_use( 'sophteiv' ) ) CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d ) 913 872 ! 914 873 zztmp = 0.5_wp * 0.5 … … 916 875 zw2d(:,:) = 0._wp 917 876 zw3d(:,:,:) = 0._wp 918 DO jk = 1, jpkm1 919 DO jj = 2, jpjm1 920 DO ji = fs_2, fs_jpim1 ! vector opt. 921 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) & 922 & * ( tsn (ji,jj,jk,jp_sal) + tsn (ji+1,jj,jk,jp_sal) ) 923 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 924 END DO 925 END DO 926 END DO 877 DO_3D_00_00( 1, jpkm1 ) 878 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1) - psi_uw(ji ,jj,jk) ) & 879 & * ( ts (ji,jj,jk,jp_sal,Kmm) + ts (ji+1,jj,jk,jp_sal,Kmm) ) 880 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 881 END_3D 927 882 CALL lbc_lnk( 'ldftra', zw2d, 'U', -1. ) 928 883 CALL lbc_lnk( 'ldftra', zw3d, 'U', -1. ) 929 884 CALL iom_put( "ueiv_salttr", zztmp * zw2d ) ! salt transport in i-direction 930 CALL iom_put( "ueiv_salttr3d", zztmp * zw3d ) 885 CALL iom_put( "ueiv_salttr3d", zztmp * zw3d ) ! salt transport in i-direction 931 886 ENDIF 932 887 zw2d(:,:) = 0._wp 933 888 zw3d(:,:,:) = 0._wp 934 DO jk = 1, jpkm1 935 DO jj = 2, jpjm1 936 DO ji = fs_2, fs_jpim1 ! vector opt. 937 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) & 938 & * ( tsn (ji,jj,jk,jp_sal) + tsn (ji,jj+1,jk,jp_sal) ) 939 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 940 END DO 941 END DO 942 END DO 889 DO_3D_00_00( 1, jpkm1 ) 890 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj ,jk) ) & 891 & * ( ts (ji,jj,jk,jp_sal,Kmm) + ts (ji,jj+1,jk,jp_sal,Kmm) ) 892 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 893 END_3D 943 894 CALL lbc_lnk( 'ldftra', zw2d, 'V', -1. ) 944 895 CALL iom_put( "veiv_salttr", zztmp * zw2d ) ! salt transport in j-direction 945 896 CALL iom_put( "veiv_salttr", zztmp * zw3d ) ! salt transport in j-direction 946 897 ! 947 IF( ln_diaptr) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )898 IF( iom_use( 'sopsteiv' ) ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d ) 948 899 ! 949 900 !
Note: See TracChangeset
for help on using the changeset viewer.
