- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- 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) (24 diffs)
-
src/OCE/LDF/ldftra.F90 (modified) (19 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/LDF/ldfc1d_c2d.F90
r10425 r12928 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/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/LDF/ldfdyn.F90
r12178 r12928 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 ) ! Namelist namdyn_ldf in reference namelist : Lateral physics118 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 ) ! Namelist namdyn_ldf in configuration namelist : Lateral physics122 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 … … 417 406 ! 418 407 zcmsmag = (rn_csmc/rpi)**2 ! (C_smag/pi)^2 419 zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 4._wp * zcmsmag )! lower limit stability factor scaling420 zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * r dt ) ! upper limit stability factor scaling408 zstabf_lo = rn_minfac * rn_minfac / ( 2._wp * 12._wp * 12._wp * zcmsmag ) ! lower limit stability factor scaling 409 zstabf_up = rn_maxfac / ( 4._wp * zcmsmag * 2._wp * rn_Dt ) ! upper limit stability factor scaling 421 410 IF( ln_dynldf_blp ) zstabf_lo = ( 16._wp / 9._wp ) * zstabf_lo ! provide |U|L^3/12 lower limit instead 422 411 ! ! of |U|L^3/16 in blp case 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/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/LDF/ldfslp.F90
r10425 r12928 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 ! ! uslp and vslp output in zwz and zww, resp. 210 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 211 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 212 ! thickness of water column between surface and level k at u/v point 213 zdepu = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji+1,jj,jk) ) & 214 - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj) ) - e3u_n(ji,jj,miku(ji,jj)) ) 215 zdepv = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji,jj+1,jk) ) & 216 - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) - e3v_n(ji,jj,mikv(ji,jj)) ) 217 ! 218 zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps ) & 219 & + zfi * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk) 220 zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps ) & 221 & + 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) 222 216 !!gm modif to suppress omlmask.... (as in Griffies case) 223 217 ! ! ! jk must be >= ML level for zf=1. otherwise zf=0. 224 218 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 225 219 ! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 226 ! zci = 0.5 * ( gdept _n(ji+1,jj,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )227 ! 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. ) ) 228 222 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 229 223 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 230 224 !!gm end modif 231 END DO 232 END DO 233 END DO 225 END_3D 234 226 CALL lbc_lnk_multi( 'ldfslp', zwz, 'U', -1., zww, 'V', -1. ) ! lateral boundary conditions 235 227 ! 236 228 ! !* horizontal Shapiro filter 237 229 DO jk = 2, jpkm1 238 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 239 DO ji = 2, jpim1 240 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 241 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 242 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 243 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 244 & + 4.* zwz(ji ,jj ,jk) ) 245 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 246 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 247 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 248 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 249 & + 4.* zww(ji,jj ,jk) ) 250 END DO 251 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 252 242 DO jj = 3, jpj-2 ! other rows 253 DO ji = fs_2, fs_jpim1 ! vector opt.243 DO ji = 2, jpim1 ! vector opt. 254 244 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 255 245 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & … … 265 255 END DO 266 256 ! !* decrease along coastal boundaries 267 DO jj = 2, jpjm1 268 DO ji = fs_2, fs_jpim1 ! vector opt. 269 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 270 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 271 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 272 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 273 END DO 274 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 275 263 END DO 276 264 … … 279 267 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 280 268 ! 281 DO jk = 2, jpkm1 282 DO jj = 2, jpjm1 283 DO ji = fs_2, fs_jpim1 ! vector opt. 284 ! !* Local vertical density gradient evaluated from N^2 285 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 286 ! !* Slopes at w point 287 ! ! i- & j-gradient of density at w-points 288 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 289 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 290 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 291 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 292 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 293 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * wmask (ji,jj,jk) 294 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 295 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * wmask (ji,jj,jk) 296 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 297 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 298 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zai ) ) 299 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zaj ) ) 300 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 301 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 302 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 ) 303 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 304 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) 305 291 306 292 !!gm modif to suppress omlmask.... (as in Griffies operator) … … 311 297 ! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 312 298 !!gm end modif 313 END DO 314 END DO 315 END DO 299 END_3D 316 300 CALL lbc_lnk_multi( 'ldfslp', zwz, 'T', -1., zww, 'T', -1. ) ! lateral boundary conditions 317 301 ! 318 302 ! !* horizontal Shapiro filter 319 303 DO jk = 2, jpkm1 320 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 321 DO ji = 2, jpim1 322 zcofw = wmask(ji,jj,jk) * z1_16 323 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 324 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 325 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 326 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 327 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 328 329 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 330 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 331 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 332 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 333 & + 4.* zww(ji ,jj ,jk) ) * zcofw 334 END DO 335 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 336 318 DO jj = 3, jpj-2 ! other rows 337 DO ji = fs_2, fs_jpim1 ! vector opt.319 DO ji = 2, jpim1 ! vector opt. 338 320 zcofw = wmask(ji,jj,jk) * z1_16 339 321 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & … … 351 333 END DO 352 334 ! !* decrease in vicinity of topography 353 DO jj = 2, jpjm1 354 DO ji = fs_2, fs_jpim1 ! vector opt. 355 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 356 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 357 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 358 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 359 END DO 360 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 361 341 END DO 362 342 … … 365 345 CALL lbc_lnk_multi( 'ldfslp', uslp , 'U', -1. , vslp , 'V', -1. , wslpi, 'W', -1., wslpj, 'W', -1. ) 366 346 367 IF( ln_ctl) THEN347 IF(sn_cfctl%l_prtctl) THEN 368 348 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 369 349 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) … … 375 355 376 356 377 SUBROUTINE ldf_slp_triad ( kt )357 SUBROUTINE ldf_slp_triad ( kt, Kbb, Kmm ) 378 358 !!---------------------------------------------------------------------- 379 359 !! *** ROUTINE ldf_slp_triad *** … … 390 370 !!---------------------------------------------------------------------- 391 371 INTEGER, INTENT( in ) :: kt ! ocean time-step index 372 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 392 373 !! 393 374 INTEGER :: ji, jj, jk, jl, ip, jp, kp ! dummy loop indices … … 402 383 REAL(wp) :: zbeta0, ze3_e1, ze3_e2 403 384 REAL(wp), DIMENSION(jpi,jpj) :: z1_mlbw 404 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zalbet405 385 REAL(wp), DIMENSION(jpi,jpj,jpk,0:1) :: zdxrho , zdyrho, zdzrho ! Horizontal and vertical density gradients 406 386 REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) :: zti_mlb, ztj_mlb ! for Griffies operator only … … 416 396 ! 417 397 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln) 418 DO jk = 1, jpkm1 ! done each pair of triad 419 DO jj = 1, jpjm1 ! NB: not masked ==> a minimum value is set 420 DO ji = 1, fs_jpim1 ! vector opt. 421 zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! i-gradient of T & S at u-point 422 zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 423 zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! j-gradient of T & S at v-point 424 zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 425 zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj) 426 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj) 427 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 428 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 429 END DO 430 END DO 431 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 432 408 ! 433 409 IF( ln_zps .AND. l_grad_zps ) THEN ! partial steps: correction of i- & j-grad on bottom 434 DO jj = 1, jpjm1 435 DO ji = 1, jpim1 436 iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points) 437 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 438 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 439 zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj) 440 zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj) 441 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 442 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 443 END DO 444 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 445 419 ENDIF 446 420 ! … … 448 422 449 423 DO kp = 0, 1 !== unmasked before density i- j-, k-gradients ==! 450 DO jk = 1, jpkm1 ! done each pair of triad 451 DO jj = 1, jpj ! NB: not masked ==> a minimum value is set 452 DO ji = 1, jpi ! vector opt. 453 IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp 454 zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) ) 455 zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) ) 456 ELSE 457 zdkt = 0._wp ! 1st level gradient set to zero 458 zdks = 0._wp 459 ENDIF 460 zdzrho_raw = ( - rab_b(ji,jj,jk+kp,jp_tem) * zdkt & 461 & + rab_b(ji,jj,jk+kp,jp_sal) * zdks & 462 & ) / e3w_n(ji,jj,jk+kp) 463 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln 464 END DO 465 END DO 466 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 467 437 END DO 468 438 ! 469 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==! 470 DO ji = 1, jpi 471 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth 472 z1_mlbw(ji,jj) = 1._wp / gdepw_n(ji,jj,jk) 473 END DO 474 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 475 443 ! 476 444 ! !== intialisations to zero ==! … … 489 457 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 490 458 DO kp = 0, 1 ! with only the slope-max limit and MASKED 491 DO jj = 1, jpjm1 492 DO ji = 1, fs_jpim1 493 ip = jl ; jp = jl 494 ! 495 jk = nmln(ji+ip,jj) + 1 496 IF( jk > mbkt(ji+ip,jj) ) THEN ! ML reaches bottom 497 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 498 ELSE 499 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 500 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 501 & - ( gdept_n(ji+1,jj,jk-kp) - gdept_n(ji,jj,jk-kp) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk) 502 ze3_e1 = e3w_n(ji+ip,jj,jk-kp) * r1_e1u(ji,jj) 503 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1 , ABS( zti_g_raw ) ), zti_g_raw ) 504 ENDIF 505 ! 506 jk = nmln(ji,jj+jp) + 1 507 IF( jk > mbkt(ji,jj+jp) ) THEN !ML reaches bottom 508 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp 509 ELSE 510 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 511 & - ( gdept_n(ji,jj+1,jk-kp) - gdept_n(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 512 ze3_e2 = e3w_n(ji,jj+jp,jk-kp) / e2v(ji,jj) 513 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2 , ABS( ztj_g_raw ) ), ztj_g_raw ) 514 ENDIF 515 END DO 516 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 517 483 END DO 518 484 END DO … … 528 494 ! Must mask contribution to slope from dz/dx at constant s for triads jk=1,kp=0 that poke up though ocean surface 529 495 znot_thru_surface = REAL( 1-1/(jk+kp), wp ) !jk+kp=1,=0.; otherwise=1.0 530 DO jj = 1, jpjm1 531 DO ji = 1, fs_jpim1 ! vector opt. 532 ! 533 ! Calculate slope relative to geopotentials used for GM skew fluxes 534 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 535 ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 536 ! masked by umask taken at the level of dz(rho) 537 ! 538 ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti) 539 ! 540 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 541 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 542 ! 543 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 544 zti_coord = znot_thru_surface * ( gdept_n(ji+1,jj ,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) 545 ztj_coord = znot_thru_surface * ( gdept_n(ji ,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) ! unmasked 546 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 547 ztj_g_raw = ztj_raw - ztj_coord 548 ! additional limit required in bilaplacian case 549 ze3_e1 = e3w_n(ji+ip,jj ,jk+kp) * r1_e1u(ji,jj) 550 ze3_e2 = e3w_n(ji ,jj+jp,jk+kp) * r1_e2v(ji,jj) 551 ! NB: hard coded factor 5 (can be a namelist parameter...) 552 zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) 553 ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw ) 554 ! 555 ! Below ML use limited zti_g as is & mask 556 ! Inside ML replace by linearly reducing sx_mlb towards surface & mask 557 ! 558 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 559 zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1 560 ! ! otherwise zfact=0 561 zti_g_lim = ( zfacti * zti_g_lim & 562 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 563 & * gdepw_n(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 564 ztj_g_lim = ( zfactj * ztj_g_lim & 565 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 566 & * gdepw_n(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 567 ! 568 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim 569 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim 570 ! 571 ! Get coefficients of isoneutral diffusion tensor 572 ! 1. Utilise gradients *relative* to s-coordinate, so add t-point slopes (*subtract* depth gradients) 573 ! 2. We require that isoneutral diffusion gives no vertical buoyancy flux 574 ! i.e. 33 term = (real slope* 31, 13 terms) 575 ! To do this, retain limited sx**2 in vertical flux, but divide by real slope for 13/31 terms 576 ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 577 ! 578 zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces 579 ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 580 ! 581 IF( ln_triad_iso ) THEN 582 zti_raw = zti_lim*zti_lim / zti_raw 583 ztj_raw = ztj_lim*ztj_lim / ztj_raw 584 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 585 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 586 zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw 587 ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw 588 ENDIF 589 ! ! switching triad scheme 590 zisw = (1._wp - rn_sw_triad ) + rn_sw_triad & 591 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) ) ) 592 zjsw = (1._wp - rn_sw_triad ) + rn_sw_triad & 593 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) ) ) 594 ! 595 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim * zisw 596 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 597 ! 598 zbu = e1e2u(ji ,jj ) * e3u_n(ji ,jj ,jk ) 599 zbv = e1e2v(ji ,jj ) * e3v_n(ji ,jj ,jk ) 600 zbti = e1e2t(ji+ip,jj ) * e3w_n(ji+ip,jj ,jk+kp) 601 zbtj = e1e2t(ji ,jj+jp) * e3w_n(ji ,jj+jp,jk+kp) 602 ! 603 wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked 604 wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim 605 END DO 606 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 607 571 END DO 608 572 END DO … … 618 582 619 583 620 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 ) 621 585 !!---------------------------------------------------------------------- 622 586 !! *** ROUTINE ldf_slp_mxl *** … … 638 602 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_gru, p_grv ! i- & j-gradient of density (u- & v-pts) 639 603 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_dzr ! z-gradient of density (T-point) 604 INTEGER , INTENT(in) :: Kmm ! ocean time level indices 640 605 !! 641 606 INTEGER :: ji , jj , jk ! dummy loop indices … … 658 623 ! 659 624 ! !== surface mixed layer mask ! 660 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise 661 DO jj = 1, jpj 662 DO ji = 1, jpi 663 ik = nmln(ji,jj) - 1 664 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp 665 ELSE ; omlmask(ji,jj,jk) = 0._wp 666 ENDIF 667 END DO 668 END DO 669 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 670 631 671 632 … … 680 641 !----------------------------------------------------------------------- 681 642 ! 682 DO jj = 2, jpjm1 683 DO ji = 2, jpim1 684 ! !== Slope at u- & v-points just below the Mixed Layer ==! 685 ! 686 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 687 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 688 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 689 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 690 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 691 ! !- horizontal density gradient at u- & v-points 692 zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj) 693 zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj) 694 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 695 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 696 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,iku)* ABS( zau ) ) 697 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,ikv)* ABS( zav ) ) 698 ! !- Slope at u- & v-points (uslpml, vslpml) 699 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 700 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 701 ! 702 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! 703 ! 704 ik = MIN( nmln(ji,jj) + 1, jpk ) 705 ikm1 = MAX( 1, ik-1 ) 706 ! !- vertical density gradient for w-slope (from N^2) 707 zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 708 ! !- horizontal density i- & j-gradient at w-points 709 zci = MAX( umask(ji-1,jj,ik ) + umask(ji,jj,ik ) & 710 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 711 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 712 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) 713 zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) & 714 & + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1 ) ) / zci * tmask(ji,jj,ik) 715 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & 716 & + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik) 717 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 718 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 719 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zai ) ) 720 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zaj ) ) 721 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 722 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 723 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 724 END DO 725 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 726 685 !!gm this lbc_lnk should be useless.... 727 686 CALL lbc_lnk_multi( 'ldfslp', uslpml , 'U', -1. , vslpml , 'V', -1. , wslpiml, 'W', -1. , wslpjml, 'W', -1. ) … … 785 744 ! DO jk = 1, jpk 786 745 ! DO jj = 2, jpjm1 787 ! DO ji = fs_2, fs_jpim1 ! vector opt.788 ! uslp (ji,jj,jk) = - ( gdept _n(ji+1,jj,jk) - gdept_n(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)789 ! vslp (ji,jj,jk) = - ( gdept _n(ji,jj+1,jk) - gdept_n(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)790 ! 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.5791 ! 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 792 751 ! END DO 793 752 ! END DO -
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/LDF/ldftra.F90
r12178 r12928 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 ) ! Namelist namtra_ldf in reference namelist : Lateral physics on tracers155 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 REWIND( numnam_cfg ) ! Namelist namtra_ldf in configuration namelist : Lateral physics on tracers158 156 READ ( numnam_cfg, namtra_ldf, IOSTAT = ios, ERR = 902 ) 159 157 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_ldf in configuration namelist' ) … … 382 380 383 381 384 SUBROUTINE ldf_tra( kt )382 SUBROUTINE ldf_tra( kt, Kbb, Kmm ) 385 383 !!---------------------------------------------------------------------- 386 384 !! *** ROUTINE ldf_tra *** … … 403 401 !!---------------------------------------------------------------------- 404 402 INTEGER, INTENT(in) :: kt ! time step 403 INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices 405 404 ! 406 405 INTEGER :: ji, jj, jk ! dummy loop indices … … 411 410 ! ! =F(growth rate of baroclinic instability) 412 411 ! ! max value aeiv_0 ; decreased to 0 within 20N-20S 413 CALL ldf_eiv( kt, aei0, aeiu, aeiv )412 CALL ldf_eiv( kt, aei0, aeiu, aeiv, Kmm ) 414 413 ENDIF 415 414 ! … … 424 423 ahtv(:,:,1) = aeiv(:,:,1) 425 424 ELSE ! compute aht. 426 CALL ldf_eiv( kt, aht0, ahtu, ahtv )425 CALL ldf_eiv( kt, aht0, ahtu, ahtv, Kmm ) 427 426 ENDIF 428 427 ! … … 430 429 zaht_min = 0.2_wp * aht0 ! minimum value for aht 431 430 zDaht = aht0 - zaht_min 432 DO jj = 1, jpj 433 DO ji = 1, jpi 434 !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg) 435 !! ==>>> The Coriolis value is identical for t- & u_points, and for v- and f-points 436 zaht = ( 1._wp - MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * zDaht 437 zahf = ( 1._wp - MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * zDaht 438 ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) ! min value zaht_min 439 ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zahf ) ! increase within 20S-20N 440 END DO 441 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 442 439 DO jk = 1, jpkm1 ! deeper value = surface value + mask for all levels 443 440 ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk) … … 448 445 IF( ln_traldf_lap ) THEN ! laplacian operator |u| e /12 449 446 DO jk = 1, jpkm1 450 ahtu(:,:,jk) = ABS( u b(:,:,jk) ) * e1u(:,:) * r1_12 ! n.b. ub,vbare masked451 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 452 449 END DO 453 450 ELSEIF( ln_traldf_blp ) THEN ! bilaplacian operator sqrt( |u| e^3 /12 ) = sqrt( |u| e /12 ) * e 454 451 DO jk = 1, jpkm1 455 ahtu(:,:,jk) = SQRT( ABS( u b(:,:,jk) ) * e1u(:,:) * r1_12 ) * e1u(:,:)456 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(:,:) 457 454 END DO 458 455 ENDIF … … 510 507 ENDIF 511 508 ! 512 REWIND( numnam_ref ) ! Namelist namtra_eiv in reference namelist : eddy induced velocity param.513 509 READ ( numnam_ref, namtra_eiv, IOSTAT = ios, ERR = 901) 514 510 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_eiv in reference namelist' ) 515 511 ! 516 REWIND( numnam_cfg ) ! Namelist namtra_eiv in configuration namelist : eddy induced velocity param.517 512 READ ( numnam_cfg, namtra_eiv, IOSTAT = ios, ERR = 902 ) 518 513 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_eiv in configuration namelist' ) … … 625 620 626 621 627 SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv )622 SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv, Kmm ) 628 623 !!---------------------------------------------------------------------- 629 624 !! *** ROUTINE ldf_eiv *** … … 637 632 !!---------------------------------------------------------------------- 638 633 INTEGER , INTENT(in ) :: kt ! ocean time-step index 634 INTEGER , INTENT(in ) :: Kmm ! ocean time level indices 639 635 REAL(wp) , INTENT(inout) :: paei0 ! max value [m2/s] 640 636 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: paeiu, paeiv ! eiv coefficient [m2/s] … … 651 647 ! ! Compute lateral diffusive coefficient at T-point 652 648 IF( ln_traldf_triad ) THEN 653 DO jk = 1, jpk 654 DO jj = 2, jpjm1 655 DO ji = 2, jpim1 656 ! Take the max of N^2 and zero then take the vertical sum 657 ! of the square root of the resulting N^2 ( required to compute 658 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 659 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 660 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk) 661 ! Compute elements required for the inverse time scale of baroclinic 662 ! eddies using the isopycnal slopes calculated in ldfslp.F : 663 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 664 ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk) 665 zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w 666 zhw(ji,jj) = zhw(ji,jj) + ze3w 667 END DO 668 END DO 669 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 670 662 ELSE 671 DO jk = 1, jpk 672 DO jj = 2, jpjm1 673 DO ji = 2, jpim1 674 ! Take the max of N^2 and zero then take the vertical sum 675 ! of the square root of the resulting N^2 ( required to compute 676 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 677 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 678 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk) 679 ! Compute elements required for the inverse time scale of baroclinic 680 ! eddies using the isopycnal slopes calculated in ldfslp.F : 681 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 682 ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk) 683 zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 684 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w 685 zhw(ji,jj) = zhw(ji,jj) + ze3w 686 END DO 687 END DO 688 END DO 689 ENDIF 690 691 DO jj = 2, jpjm1 692 DO ji = fs_2, fs_jpim1 ! vector opt. 693 zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 694 ! Rossby radius at w-point taken betwenn 2 km and 40km 695 zRo(ji,jj) = MAX( 2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 ) ) 696 ! Compute aeiw by multiplying Ro^2 and T^-1 697 zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 698 END DO 699 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 700 686 701 687 ! !== Bound on eiv coeff. ==! 702 688 z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) ) 703 DO jj = 2, jpjm1 704 DO ji = fs_2, fs_jpim1 ! vector opt. 705 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease 706 zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0 707 END DO 708 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 709 693 CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1. ) ! lateral boundary condition 710 694 ! 711 DO jj = 2, jpjm1 !== aei at u- and v-points ==! 712 DO ji = fs_2, fs_jpim1 ! vector opt. 713 paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj ) ) * umask(ji,jj,1) 714 paeiv(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji ,jj+1) ) * vmask(ji,jj,1) 715 END DO 716 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 717 699 CALL lbc_lnk_multi( 'ldftra', paeiu(:,:,1), 'U', 1. , paeiv(:,:,1), 'V', 1. ) ! lateral boundary condition 718 700 … … 725 707 726 708 727 SUBROUTINE ldf_eiv_trp( kt, kit000, pu n, pvn, pwn, cdtype)709 SUBROUTINE ldf_eiv_trp( kt, kit000, pu, pv, pw, cdtype, Kmm, Krhs ) 728 710 !!---------------------------------------------------------------------- 729 711 !! *** ROUTINE ldf_eiv_trp *** … … 741 723 !! velocity and heat transport (call ldf_eiv_dia) 742 724 !! 743 !! ** Action : pun, pvn increased by the eiv transport 744 !!---------------------------------------------------------------------- 745 INTEGER , INTENT(in ) :: kt ! ocean time-step index 746 INTEGER , INTENT(in ) :: kit000 ! first time step index 747 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 748 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pun ! in : 3 ocean transport components [m3/s] 749 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pvn ! out: 3 ocean transport components [m3/s] 750 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] 751 734 !! 752 735 INTEGER :: ji, jj, jk ! dummy loop indices … … 766 749 zpsi_uw(:,:,jpk) = 0._wp ; zpsi_vw(:,:,jpk) = 0._wp 767 750 ! 768 DO jk = 2, jpkm1 769 DO jj = 1, jpjm1 770 DO ji = 1, fs_jpim1 ! vector opt. 771 zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk ) + wslpi(ji+1,jj,jk) ) & 772 & * ( aeiu (ji,jj,jk-1) + aeiu (ji ,jj,jk) ) * umask(ji,jj,jk) 773 zpsi_vw(ji,jj,jk) = - r1_4 * e1v(ji,jj) * ( wslpj(ji,jj,jk ) + wslpj(ji,jj+1,jk) ) & 774 & * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj ,jk) ) * vmask(ji,jj,jk) 775 END DO 776 END DO 777 END DO 778 ! 779 DO jk = 1, jpkm1 780 DO jj = 1, jpjm1 781 DO ji = 1, fs_jpim1 ! vector opt. 782 pun(ji,jj,jk) = pun(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 783 pvn(ji,jj,jk) = pvn(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 784 END DO 785 END DO 786 END DO 787 DO jk = 1, jpkm1 788 DO jj = 2, jpjm1 789 DO ji = fs_2, fs_jpim1 ! vector opt. 790 pwn(ji,jj,jk) = pwn(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj ,jk) & 791 & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji ,jj-1,jk) ) 792 END DO 793 END DO 794 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 795 766 ! 796 767 ! ! diagnose the eddy induced velocity and associated heat transport 797 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 ) 798 769 ! 799 770 END SUBROUTINE ldf_eiv_trp 800 771 801 772 802 SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw )773 SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw, Kmm ) 803 774 !!---------------------------------------------------------------------- 804 775 !! *** ROUTINE ldf_eiv_dia *** … … 811 782 !!---------------------------------------------------------------------- 812 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 813 785 ! 814 786 INTEGER :: ji, jj, jk ! dummy loop indices … … 831 803 ! 832 804 DO jk = 1, jpkm1 ! e2u e3u u_eiv = -dk[psi_uw] 833 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) ) 834 806 END DO 835 807 CALL iom_put( "uoce_eiv", zw3d ) 836 808 ! 837 809 DO jk = 1, jpkm1 ! e1v e3v v_eiv = -dk[psi_vw] 838 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) ) 839 811 END DO 840 812 CALL iom_put( "voce_eiv", zw3d ) 841 813 ! 842 DO jk = 1, jpkm1 ! e1 e2 w_eiv = dk[psix] + dk[psix] 843 DO jj = 2, jpjm1 844 DO ji = fs_2, fs_jpim1 ! vector opt. 845 zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk) - psi_vw(ji ,jj-1,jk) & 846 & + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj ,jk) ) / e1e2t(ji,jj) 847 END DO 848 END DO 849 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 850 818 CALL lbc_lnk( 'ldftra', zw3d, 'T', 1. ) ! lateral boundary condition 851 819 CALL iom_put( "woce_eiv", zw3d ) 852 820 ! 853 ! 854 zztmp = 0.5_wp * rau0 * rcp 821 IF( iom_use('weiv_masstr') ) THEN ! vertical mass transport & its square value 822 zw2d(:,:) = rho0 * e1e2t(:,:) 823 DO jk = 1, jpk 824 zw3d(:,:,jk) = zw3d(:,:,jk) * zw2d(:,:) 825 END DO 826 CALL iom_put( "weiv_masstr" , zw3d ) 827 ENDIF 828 ! 829 IF( iom_use('ueiv_masstr') ) THEN 830 zw3d(:,:,:) = 0.e0 831 DO jk = 1, jpkm1 832 zw3d(:,:,jk) = rho0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) 833 END DO 834 CALL iom_put( "ueiv_masstr", zw3d ) ! mass transport in i-direction 835 ENDIF 836 ! 837 zztmp = 0.5_wp * rho0 * rcp 855 838 IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN 856 839 zw2d(:,:) = 0._wp 857 840 zw3d(:,:,:) = 0._wp 858 DO jk = 1, jpkm1 859 DO jj = 2, jpjm1 860 DO ji = fs_2, fs_jpim1 ! vector opt. 861 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) & 862 & * ( tsn (ji,jj,jk,jp_tem) + tsn (ji+1,jj,jk,jp_tem) ) 863 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 864 END DO 865 END DO 866 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 867 846 CALL lbc_lnk( 'ldftra', zw2d, 'U', -1. ) 868 847 CALL lbc_lnk( 'ldftra', zw3d, 'U', -1. ) … … 870 849 CALL iom_put( "ueiv_heattr3d", zztmp * zw3d ) ! heat transport in i-direction 871 850 ENDIF 851 ! 852 IF( iom_use('veiv_masstr') ) THEN 853 zw3d(:,:,:) = 0.e0 854 DO jk = 1, jpkm1 855 zw3d(:,:,jk) = rho0 * ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) 856 END DO 857 CALL iom_put( "veiv_masstr", zw3d ) ! mass transport in i-direction 858 ENDIF 859 ! 872 860 zw2d(:,:) = 0._wp 873 861 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_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) & 878 & * ( tsn (ji,jj,jk,jp_tem) + tsn (ji,jj+1,jk,jp_tem) ) 879 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 880 END DO 881 END DO 882 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 883 867 CALL lbc_lnk( 'ldftra', zw2d, 'V', -1. ) 884 868 CALL iom_put( "veiv_heattr", zztmp * zw2d ) ! heat transport in j-direction 885 869 CALL iom_put( "veiv_heattr", zztmp * zw3d ) ! heat transport in j-direction 886 870 ! 887 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 ) 888 872 ! 889 873 zztmp = 0.5_wp * 0.5 … … 891 875 zw2d(:,:) = 0._wp 892 876 zw3d(:,:,:) = 0._wp 893 DO jk = 1, jpkm1 894 DO jj = 2, jpjm1 895 DO ji = fs_2, fs_jpim1 ! vector opt. 896 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) & 897 & * ( tsn (ji,jj,jk,jp_sal) + tsn (ji+1,jj,jk,jp_sal) ) 898 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 899 END DO 900 END DO 901 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 902 882 CALL lbc_lnk( 'ldftra', zw2d, 'U', -1. ) 903 883 CALL lbc_lnk( 'ldftra', zw3d, 'U', -1. ) 904 884 CALL iom_put( "ueiv_salttr", zztmp * zw2d ) ! salt transport in i-direction 905 CALL iom_put( "ueiv_salttr3d", zztmp * zw3d ) 885 CALL iom_put( "ueiv_salttr3d", zztmp * zw3d ) ! salt transport in i-direction 906 886 ENDIF 907 887 zw2d(:,:) = 0._wp 908 888 zw3d(:,:,:) = 0._wp 909 DO jk = 1, jpkm1 910 DO jj = 2, jpjm1 911 DO ji = fs_2, fs_jpim1 ! vector opt. 912 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) & 913 & * ( tsn (ji,jj,jk,jp_sal) + tsn (ji,jj+1,jk,jp_sal) ) 914 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 915 END DO 916 END DO 917 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 918 894 CALL lbc_lnk( 'ldftra', zw2d, 'V', -1. ) 919 895 CALL iom_put( "veiv_salttr", zztmp * zw2d ) ! salt transport in j-direction 920 896 CALL iom_put( "veiv_salttr", zztmp * zw3d ) ! salt transport in j-direction 921 897 ! 922 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 ) 923 899 ! 924 900 !
Note: See TracChangeset
for help on using the changeset viewer.
