- Timestamp:
- 2010-11-10T16:38:27+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r2344 r2371 9 9 !! NEMO 0.5 ! 2002-10 (G. Madec) Free form, F90 10 10 !! 1.0 ! 2005-10 (A. Beckmann) correction for s-coordinates 11 !! 3.3 ! 20 06-10 (C. Harris, G. Nurser) add ldf_slp_grif (Griffies operator)11 !! 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) add Griffies operator 12 12 !!---------------------------------------------------------------------- 13 13 #if defined key_ldfslp || defined key_esopa … … 15 15 !! 'key_ldfslp' Rotation of lateral mixing tensor 16 16 !!---------------------------------------------------------------------- 17 !! ldf_slp : compute the slopes of neutral surface 18 !! ldf_slp_mxl : compute the slopes of iso-neutral surface 17 !! ldf_slp_grif : calculates the triads of isoneutral slopes (Griffies operator) 18 !! ldf_slp : calculates the slopes of neutral surface (Madec operator) 19 !! ldf_slp_mxl : calculates the slopes at the base of the mixed layer (Madec operator) 19 20 !! ldf_slp_init : initialization of the slopes computation 20 21 !!---------------------------------------------------------------------- 21 22 USE oce ! ocean dynamics and tracers 22 23 USE dom_oce ! ocean space and time domain 23 USE ldftra_oce 24 USE ldfdyn_oce 24 USE ldftra_oce ! lateral diffusion: traceur 25 USE ldfdyn_oce ! lateral diffusion: dynamics 25 26 USE phycst ! physical constants 26 27 USE zdfmxl ! mixed layer depth 27 USE eosbn2 28 USE eosbn2 ! equation of states 28 29 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 30 USE in_out_manager ! I/O manager … … 33 34 PRIVATE 34 35 35 PUBLIC ldf_slp ! routine called by step.F90 36 PUBLIC ldf_slp_init ! routine called by opa.F90 37 PUBLIC ldf_slp_grif ! routine called by step.F90 38 39 LOGICAL , PUBLIC, PARAMETER :: lk_ldfslp = .TRUE. !: slopes flag 40 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: uslp, wslpi !: i_slope at U- and W-points 41 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: vslp, wslpj !: j-slope at V- and W-points 42 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: wslp2 !: wslp**2 from Griffies quarter cells 43 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: alpha, beta !: alpha,beta at T points 44 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: tfw,sfw,ftu,fsu,ftv,fsv,ftud,fsud,ftvd,fsvd 45 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: psix_eiv 46 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: psiy_eiv 47 48 REAL(wp), DIMENSION(jpi,jpj,jpk) :: omlmask ! mask of the surface mixed layer at T-pt 49 REAL(wp), DIMENSION(jpi,jpj) :: uslpml, wslpiml ! i_slope at U- and W-points just below the mixed layer 50 REAL(wp), DIMENSION(jpi,jpj) :: vslpml, wslpjml ! j_slope at V- and W-points just below the mixed layer 36 PUBLIC ldf_slp ! routine called by step.F90 37 PUBLIC ldf_slp_grif ! routine called by step.F90 38 PUBLIC ldf_slp_init ! routine called by opa.F90 39 40 LOGICAL , PUBLIC, PARAMETER :: lk_ldfslp = .TRUE. !: slopes flag 41 ! !! Madec operator 42 REAL(wp), PUBLIC, DIMENSION(:,:,:) , ALLOCATABLE :: uslp, wslpi !: i_slope at U- and W-points 43 REAL(wp), PUBLIC, DIMENSION(:,:,:) , ALLOCATABLE :: vslp, wslpj !: j-slope at V- and W-points 44 ! !! Griffies operator 45 REAL(wp), PUBLIC, DIMENSION(:,:,:) , ALLOCATABLE :: wslp2 !: wslp**2 from Griffies quarter cells 46 REAL(wp), PUBLIC, DIMENSION(:,:,:,:,:), ALLOCATABLE :: triadi_g, triadj_g !: skew flux slopes relative to geopotentials 47 REAL(wp), PUBLIC, DIMENSION(:,:,:,:,:), ALLOCATABLE :: triadi , triadj !: isoneutral slopes relative to model-coordinate 48 49 ! !! Madec operator 50 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: omlmask ! mask of the surface mixed layer at T-pt 51 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: uslpml, wslpiml ! i_slope at U- and W-points just below the mixed layer 52 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: vslpml, wslpjml ! j_slope at V- and W-points just below the mixed layer 53 54 REAL(wp) :: repsln = 1.e-25_wp ! tiny value used as minium of di(rho), dj(rho) and dk(rho) 51 55 52 56 !! * Substitutions … … 67 71 !! 68 72 !! ** Purpose : Compute the slopes of neutral surface (slope of isopycnal 69 !! surfaces referenced locally) ( 'key_traldfiso').73 !! surfaces referenced locally) (ln_traldf_iso=T). 70 74 !! 71 75 !! ** Method : The slope in the i-direction is computed at U- and … … 347 351 ! 348 352 END SUBROUTINE ldf_slp 349 353 350 354 351 355 SUBROUTINE ldf_slp_grif ( kt ) 352 !!---------------------------------------------------------------------- 353 !! *** ROUTINE ldf_slp_grif *** 354 !! 355 !! ** Purpose : Compute the squared slopes of neutral surfaces (slope 356 !! of iso-pycnal surfaces referenced locally) ('key_traldfiso') 357 !! at W-points using the Griffies quarter-cells. Also calculates 358 !! alpha and beta at T-points for use in the Griffies isopycnal 359 !! scheme. 360 !! 361 !! ** Method : The slope in the i-direction is computed at U- and 362 !! W-points (uslp, wslpi) and the slope in the j-direction is 363 !! computed at V- and W-points (vslp, wslpj). 364 !! 365 !! ** Action : - alpha, beta 366 !! wslp2 squared slope of neutral surfaces at w-points. 367 !!---------------------------------------------------------------------- 368 USE oce, zdit => ua ! use ua as workspace 369 USE oce, zdis => va ! use va as workspace 370 USE oce, zdjt => ta ! use ta as workspace 371 USE oce, zdjs => sa ! use sa as workspace 372 !! 373 INTEGER, INTENT( in ) :: kt ! ocean time-step index 374 !! 375 INTEGER :: ji, jj, jk, ip, jp, kp ! dummy loop indices 376 INTEGER :: iku, ikv ! local integer 377 REAL(wp) :: zt, zs, zh, zt2, zsp5, zp1t1 ! local scalars 378 REAL(wp) :: zdenr, zrhotmp, zdndt, zdddt ! - - 379 REAL(wp) :: zdnds, zddds, znum, zden ! - - 380 REAL(wp) :: zslope, za_sxe, zslopec, zdsloper ! - - 381 REAL(wp) :: zfact, zepsln, zatempw,zatempu,zatempv ! - - 382 REAL(wp) :: ze1ur, ze2vr, ze3wr, zdxt, zdxs, zdyt, zdys, zdzt, zdzs, zvolf 383 REAL(wp) :: zr_slpmax, zdxrho, zdyrho, zabs_dzrho 384 REAL(wp), DIMENSION(jpi,jpj,jpk,0:1,0:1) :: zsx,zsy 385 REAL(wp), DIMENSION(jpi,jpj ,0:1,0:1) :: zsx_ml_base, zsy_ml_base 386 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdkt, zdks 387 REAL(wp), DIMENSION(jpi,jpj) :: zr_ml_basew 388 !!---------------------------------------------------------------------- 389 390 ! 0. Local constant initialization 391 ! -------------------------------- 392 zr_slpmax = 1.0_wp/slpmax 393 394 ! zslopec=0.004 395 ! zdsloper=1000.0 396 zepsln=1e-25 397 398 SELECT CASE ( nn_eos ) 399 400 CASE ( 0 ) ! Jackett and McDougall (1994) formulation 401 402 ! ! =============== 403 DO jk = 1, jpk ! Horizontal slab 404 ! ! =============== 405 DO jj = 1, jpjm1 406 DO ji = 1, fs_jpim1 407 zt = tb(ji,jj,jk) ! potential temperature 408 zs = sb(ji,jj,jk) - 35.0 ! salinity anomaly (s-35) 409 zh = fsdept(ji,jj,jk) ! depth in meters 410 411 beta(ji,jj,jk) = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt & 412 & - 0.301985e-05 ) * zt & 413 & + 0.785567e-03 & 414 & + ( 0.515032e-08 * zs & 415 & + 0.788212e-08 * zt - 0.356603e-06 ) * zs & 416 & +( ( 0.121551e-17 * zh & 417 & - 0.602281e-15 * zs & 418 & - 0.175379e-14 * zt + 0.176621e-12 ) * zh & 419 & + 0.408195e-10 * zs & 420 & + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt & 421 & - 0.121555e-07 ) * zh 422 423 alpha(ji,jj,jk) = - beta(ji,jj,jk) * & 424 & (((( - 0.255019e-07 * zt + 0.298357e-05 ) * zt & 425 & - 0.203814e-03 ) * zt & 426 & + 0.170907e-01 ) * zt & 427 & + 0.665157e-01 & 428 & + ( - 0.678662e-05 * zs & 429 & - 0.846960e-04 * zt + 0.378110e-02 ) * zs & 430 & + ( ( - 0.302285e-13 * zh & 431 & - 0.251520e-11 * zs & 432 & + 0.512857e-12 * zt * zt ) * zh & 433 & - 0.164759e-06 * zs & 434 & +( 0.791325e-08 * zt - 0.933746e-06 ) * zt & 435 & + 0.380374e-04 ) * zh) 436 437 END DO 438 END DO 439 END DO 440 441 CASE ( 1 ) 442 443 alpha(:,:,:)=-rn_alpha 444 beta(:,:,:)=0.0 445 446 CASE ( 2 ) 447 448 alpha(:,:,:)=-rn_alpha 449 beta (:,:,:)=rn_beta 450 451 CASE ( 3 ) 452 453 DO jk =1, jpk 454 DO jj = 1, jpjm1 455 DO ji = 1, fs_jpim1 456 zt = tb(ji,jj,jk) 457 zs = sb(ji,jj,jk) 458 zh = fsdept(ji,jj,jk) 459 zt2 = zt**2 460 zsp5 = sqrt(ABS(zs)) 461 zp1t1=zh*zt 462 znum=9.99843699e+02+zt*(7.35212840e+00+zt*(-5.45928211e-02+3.98476704e-04*zt)) & 463 +zs*(2.96938239e+00-7.23268813e-03*zt+2.12382341e-03*zs) & 464 +zh*(1.04004591e-02+1.03970529e-07*zt2+5.18761880e-06*zs+ & 465 zh*(-3.24041825e-08-1.23869360e-11*zt2)) 466 zden=1.00000000e+00+zt*(7.28606739e-03+zt*(-4.60835542e-05+zt*(3.68390573e-07+zt*1.80809186e-10))) & 467 +zs*(2.14691708e-03+zt*(-9.27062484e-06-1.78343643e-10*zt2)+zsp5*(4.76534122e-06+1.63410736e-09*zt2)) & 468 + zh*(5.30848875e-06+zh*zt*(-3.03175128e-16*zt2-1.27934137e-17*zh)) 469 zdenr=1.0/zden 470 zrhotmp=znum*zdenr 471 zdndt=7.35212840e+00+zt*(-1.091856422e-01+1.195430112e-03*zt)-7.23268813e-03*zs & 472 +zp1t1*(2.07941058e-07-2.4773872e-11*zh) 473 zdddt=7.28606739e-03+zt*(-9.21671084e-05+zt*(1.105171719e-06+7.23236744e-10*zt)) & 474 +zs*(-9.27062484e-06+zt*(-5.35030929e-10*zt+3.26821472e-09*zsp5)) & 475 +zh*zh*(-9.09525384e-16*zt2-1.27934137e-17*zh) 476 zdnds=2.96938239e+00-7.23268813e-03*zt+2*2.12382341e-03*zs+5.18761880e-06*zh 477 zddds=2.14691708e-03+zt*(-9.27062484e-06-1.78343643e-10*zt2)+zsp5*(7.14801183e-06+2.45116104e-09*zt2) 478 alpha(ji,jj,jk)=(zdndt-zrhotmp*zdddt)*zdenr 479 beta(ji,jj,jk)=zdenr*(zdnds-zrhotmp*zddds) 480 481 END DO 482 END DO 483 END DO 484 485 CASE DEFAULT 486 487 IF(lwp) WRITE(numout,cform_err) 488 IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos 489 nstop = nstop + 1 490 491 END SELECT 492 493 CALL lbc_lnk( alpha, 'T', 1. ) 494 CALL lbc_lnk( beta, 'T', 1. ) 495 496 497 ! --------------------------------------- 498 ! 1. Horizontal tracer gradients at T-level jk 499 ! --------------------------------------- 500 DO jk = 1, jpkm1 501 DO jj = 1, jpjm1 502 DO ji = 1, fs_jpim1 ! vector opt. 503 ! i-gradient of T and S at jj 504 zdit (ji,jj,jk) = ( tb(ji+1,jj,jk)-tb(ji,jj,jk) ) * umask(ji,jj,jk) 505 zdis (ji,jj,jk) = ( sb(ji+1,jj,jk)-sb(ji,jj,jk) ) * umask(ji,jj,jk) 506 ! j-gradient of T and S at jj 507 zdjt (ji,jj,jk) = ( tb(ji,jj+1,jk)-tb(ji,jj,jk) ) * vmask(ji,jj,jk) 508 zdjs (ji,jj,jk) = ( sb(ji,jj+1,jk)-sb(ji,jj,jk) ) * vmask(ji,jj,jk) 509 END DO 510 END DO 511 END DO 512 513 IF( ln_zps ) THEN ! partial steps correction at the last level 356 !!---------------------------------------------------------------------- 357 !! *** ROUTINE ldf_slp_grif *** 358 !! 359 !! ** Purpose : Compute the squared slopes of neutral surfaces (slope 360 !! of iso-pycnal surfaces referenced locally) (ln_traldf_grif=T) 361 !! at W-points using the Griffies quarter-cells. 362 !! 363 !! ** Method : calculates alpha and beta at T-points 364 !! 365 !! ** Action : - triadi_g, triadj_g T-pts i- and j-slope triads relative to geopot. (used for eiv) 366 !! - triadi , triadj T-pts i- and j-slope triads relative to model-coordinate 367 !! - wslp2 squared slope of neutral surfaces at w-points. 368 !!---------------------------------------------------------------------- 369 USE oce, zdit => ua ! use ua as workspace 370 USE oce, zdis => va ! use va as workspace 371 USE oce, zdjt => ta ! use ta as workspace 372 USE oce, zdjs => sa ! use sa as workspace 373 !! 374 INTEGER, INTENT( in ) :: kt ! ocean time-step index 375 !! 376 INTEGER :: ji, jj, jk, jl, ip, jp, kp ! dummy loop indices 377 INTEGER :: iku, ikv ! temporary integer 378 REAL(wp) :: zfacti, zfactj, zatempw,zatempu,zatempv ! local scalars 379 REAL(wp) :: zbu, zbv, zbti, zbtj 380 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_lim2, zti_g_raw, zti_g_lim 381 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_lim2, ztj_g_raw, ztj_g_lim 382 REAL(wp) :: zdzrho_raw 383 REAL(wp), DIMENSION(jpi,jpj,jpk,0:1) :: zdzrho, zdyrho, zdxrho ! Horizontal and vertical density gradients 384 REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) :: zti_mlb, ztj_mlb 385 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdkt, zdks 386 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zalpha, zbeta ! alpha, beta at T points, at depth fsgdept 387 REAL(wp), DIMENSION(jpi,jpj) :: z1_mlbw 388 !!---------------------------------------------------------------------- 389 390 !--------------------------------! 391 ! Some preliminary calculation ! 392 !--------------------------------! 393 ! 394 CALL eos_alpbet( tsb, zalpha, zbeta ) !== before thermal and haline expension coeff. at T-points ==! 395 ! 396 DO jk = 1, jpkm1 !== before lateral T & S gradients at T-level jk ==! 397 DO jj = 1, jpjm1 398 DO ji = 1, fs_jpim1 ! vector opt. 399 zdit(ji,jj,jk) = ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) * umask(ji,jj,jk) ! i-gradient of T and S at jj 400 zdis(ji,jj,jk) = ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) * umask(ji,jj,jk) 401 zdjt(ji,jj,jk) = ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) * vmask(ji,jj,jk) ! j-gradient of T and S at jj 402 zdjs(ji,jj,jk) = ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) * vmask(ji,jj,jk) 403 END DO 404 END DO 405 END DO 406 IF( ln_zps ) THEN ! partial steps: correction at the last level 514 407 # if defined key_vectopt_loop 515 DO jj = 1, 1516 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)408 DO jj = 1, 1 409 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 517 410 # else 518 519 411 DO jj = 1, jpjm1 412 DO ji = 1, jpim1 520 413 # endif 521 ! last ocean level 522 iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj ) ) - 1 523 ikv = MIN( mbathy(ji,jj), mbathy(ji ,jj+1) ) - 1 524 ! i-gradient of T and S 525 zdit (ji,jj,iku) = gtsu(ji,jj,jp_tem) 526 zdis (ji,jj,iku) = gtsu(ji,jj,jp_sal) 527 ! j-gradient of T and S 528 zdjt (ji,jj,ikv) = gtsv(ji,jj,jp_tem) 529 zdjs (ji,jj,ikv) = gtsv(ji,jj,jp_sal) 530 END DO 531 END DO 532 ENDIF 533 534 ! --------------------------------------- 535 ! 1. Vertical tracer gradient at w-level jk 536 ! --------------------------------------- 537 DO jk = 2, jpk 538 zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk) 539 zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk) 540 END DO 541 542 zdkt(:,:,1) = 0.0 543 zdks(:,:,1) = 0.0 544 545 ! --------------------------------------- 546 ! Depth of the w-point below ML base 547 ! --------------------------------------- 548 DO jj = 1, jpj 549 DO ji = 1, jpi 550 jk = nmln(ji,jj) 551 zr_ml_basew(ji,jj)=1.0/fsdepw(ji,jj,jk+1) 552 END DO 553 END DO 554 555 556 wslp2(:,:,:)=0.0 557 tfw(:,:,:) = 0.0 558 sfw(:,:,:) = 0.0 559 ftu(:,:,:) = 0.0 560 fsu(:,:,:) = 0.0 561 ftv(:,:,:) = 0.0 562 fsv(:,:,:) = 0.0 563 ftud(:,:,:) = 0.0 564 fsud(:,:,:) = 0.0 565 ftvd(:,:,:) = 0.0 566 fsvd(:,:,:) = 0.0 567 psix_eiv(:,:,:) = 0.0 568 psiy_eiv(:,:,:) = 0.0 569 570 ! ---------------------------------------------------------------------- 571 ! x-z plane 572 ! ---------------------------------------------------------------------- 573 574 ! calculate limited triad x-slopes zsx in interior (1=<jk=<jpk-1) 575 DO ip=0,1 576 DO kp=0,1 577 578 DO jk = 1, jpkm1 579 DO jj = 1, jpjm1 580 DO ji = 1, fs_jpim1 ! vector opt. 581 582 ze1ur=1.0/e1u(ji,jj) 583 zdxt=zdit(ji,jj,jk)*ze1ur 584 zdxs=zdis(ji,jj,jk)*ze1ur 585 586 ze3wr=1.0/fse3w(ji+ip,jj,jk+kp) 587 zdzt=zdkt(ji+ip,jj,jk+kp)*ze3wr 588 zdzs=zdks(ji+ip,jj,jk+kp)*ze3wr 589 ! Calculate the horizontal and vertical density gradient 590 zdxrho = alpha(ji+ip,jj,jk)*zdxt+beta(ji+ip,jj,jk)*zdxs 591 zabs_dzrho = ABS(alpha(ji+ip,jj,jk)*zdzt+beta(ji+ip,jj,jk)*zdzs)+zepsln 592 ! Limit by slpmax, and mask by psi-point 593 zsx(ji+ip,jj,jk,1-ip,kp) = umask(ji,jj,jk+kp) & 594 & *zdxrho/MAX(zabs_dzrho,zr_slpmax*ABS(zdxrho)) 595 END DO 596 END DO 597 END DO 598 599 END DO 600 END DO 601 602 ! calculate slope of triad immediately below mixed-layer base 603 DO ip=0,1 604 DO kp=0,1 605 DO jj = 1, jpjm1 606 DO ji = 1, fs_jpim1 607 jk = nmln(ji+ip,jj) 608 zsx_ml_base(ji+ip,jj,1-ip,kp)=zsx(ji+ip,jj,jk+1-kp,1-ip,kp) 609 END DO 610 END DO 611 END DO 612 END DO 613 614 ! Below ML use limited zsx as is 615 ! Inside ML replace by linearly reducing zsx_ml_base towards surface 616 DO ip=0,1 617 DO kp=0,1 618 619 DO jk = 1, jpkm1 620 621 DO jj = 1, jpjm1 622 623 DO ji = 1, fs_jpim1 ! vector opt. 624 ! k index of uppermost point(s) of triad is jk+kp-1 625 ! must be .ge. nmln(ji,jj) for zfact=1. 626 ! otherwise zfact=0. 627 zfact = 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)) 628 zsx(ji+ip,jj,jk,1-ip,kp) = zfact*zsx(ji+ip,jj,jk,1-ip,kp) + & 629 & (1.0_wp-zfact)*(fsdepw(ji+ip,jj,jk+kp)*zr_ml_basew(ji+ip,jj))*zsx_ml_base(ji+ip,jj,1-ip,kp) 630 END DO 631 632 END DO 633 634 END DO 635 END DO 636 END DO 637 638 ! Use zsx to calculate fluxes and save averaged slopes psix_eiv at psi-points 639 DO ip=0,1 640 DO kp=0,1 641 642 DO jk = 1, jpkm1 643 644 DO jj = 1, jpjm1 645 646 DO ji = 1, fs_jpim1 ! vector opt. 647 648 ze1ur=1.0/e1u(ji,jj) 649 zdxt=zdit(ji,jj,jk)*ze1ur 650 zdxs=zdis(ji,jj,jk)*ze1ur 651 652 ze3wr=1.0/fse3w(ji+ip,jj,jk+kp) 653 zdzt=zdkt(ji+ip,jj,jk+kp)*ze3wr 654 zdzs=zdks(ji+ip,jj,jk+kp)*ze3wr 655 zslope=zsx(ji+ip,jj,jk,1-ip,kp) 656 657 zvolf = 0.25_wp*e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) 658 659 ftu(ji,jj,jk)= ftu(ji,jj,jk)+zslope*zdzt*zvolf*ze1ur 660 fsu(ji,jj,jk)= fsu(ji,jj,jk)+zslope*zdzs*zvolf*ze1ur 661 ftud(ji,jj,jk)=ftud(ji,jj,jk)+fsahtu(ji,jj,jk)*zdxt*zvolf*ze1ur 662 fsud(ji,jj,jk)=fsud(ji,jj,jk)+fsahtu(ji,jj,jk)*zdxs*zvolf*ze1ur 663 tfw(ji+ip,jj,jk+kp)=tfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zslope*zdxt 664 sfw(ji+ip,jj,jk+kp)=sfw(ji+ip,jj,jk+kp)+(zvolf*ze3wr)*zslope*zdxs 665 wslp2(ji+ip,jj,jk+kp)=wslp2(ji+ip,jj,jk+kp)+ & 666 & ((zvolf*ze3wr)/(e1t(ji+ip,jj)*e2t(ji+ip,jj)))*(zslope)**2 667 psix_eiv(ji,jj,jk+kp) = psix_eiv(ji,jj,jk+kp) + 0.25_wp*zslope 668 669 END DO 670 END DO 671 672 END DO 673 END DO 674 END DO 675 676 ! ---------------------------------------------------------------------- 677 ! y-z plane 678 ! ---------------------------------------------------------------------- 679 680 ! calculate limited triad y-slopes zsy in interior (1=<jk=<jpk-1) 681 DO jp=0,1 682 DO kp=0,1 683 684 DO jk = 1, jpkm1 685 DO jj = 1, jpjm1 686 DO ji = 1, fs_jpim1 ! vector opt. 687 688 ze2vr=1.0/e2v(ji,jj) 689 zdyt=zdjt(ji,jj,jk)*ze2vr 690 zdys=zdjs(ji,jj,jk)*ze2vr 691 692 ze3wr=1.0/fse3w(ji,jj+jp,jk+kp) 693 zdzt=zdkt(ji,jj+jp,jk+kp)*ze3wr 694 zdzs=zdks(ji,jj+jp,jk+kp)*ze3wr 695 ! Calculate the horizontal and vertical density gradient 696 zdyrho = alpha(ji,jj+jp,jk)*zdyt+beta(ji,jj+jp,jk)*zdys 697 zabs_dzrho = ABS(alpha(ji,jj+jp,jk)*zdzt+beta(ji,jj+jp,jk)*zdzs)+zepsln 698 ! Limit by slpmax, and mask by psi-point 699 zsy(ji,jj+jp,jk,1-jp,kp) = vmask(ji,jj,jk+kp) & 700 & *zdyrho/MAX(zabs_dzrho,zr_slpmax*ABS(zdyrho)) 701 END DO 702 END DO 703 END DO 704 705 END DO 706 END DO 707 708 ! calculate slope of triad immediately below mixed-layer base 709 DO jp=0,1 710 DO kp=0,1 711 DO jj = 1, jpjm1 712 DO ji = 1, fs_jpim1 713 jk = nmln(ji,jj+jp) 714 zsy_ml_base(ji,jj+jp,1-jp,kp)=zsy(ji,jj+jp,jk+1-kp,1-jp,kp) 715 END DO 716 END DO 717 END DO 718 END DO 719 720 ! Below ML use limited zsy as is 721 ! Inside ML replace by linearly reducing zsy_ml_base towards surface 722 DO jp=0,1 723 DO kp=0,1 724 DO jk = 1, jpkm1 725 DO jj = 1, jpjm1 726 DO ji = 1, fs_jpim1 ! vector opt. 727 ! k index of uppermost point(s) of triad is jk+kp-1 728 ! must be .ge. nmln(ji,jj) for zfact=1. 729 ! otherwise zfact=0. 730 zfact = 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)) 731 zsy(ji,jj+jp,jk,1-jp,kp) = zfact*zsy(ji,jj+jp,jk,1-jp,kp) + & 732 & (1.0_wp-zfact)*(fsdepw(ji,jj+jp,jk+kp)*zr_ml_basew(ji,jj+jp))*zsy_ml_base(ji,jj+jp,1-jp,kp) 733 END DO 734 735 END DO 736 737 END DO 738 END DO 739 END DO 740 741 ! Use zsy to calculate fluxes and save averaged slopes psiy_eiv at psi-points 742 DO jp=0,1 743 DO kp=0,1 744 DO jk = 1, jpkm1 745 DO jj = 1, jpjm1 746 DO ji = 1, fs_jpim1 ! vector opt. 747 748 ze2vr=1.0/e2v(ji,jj) 749 zdyt=zdjt(ji,jj,jk)*ze2vr 750 zdys=zdjs(ji,jj,jk)*ze2vr 751 752 ze3wr=1.0/fse3w(ji,jj+jp,jk+kp) 753 zdzt=zdkt(ji,jj+jp,jk+kp)*ze3wr 754 zdzs=zdks(ji,jj+jp,jk+kp)*ze3wr 755 zslope=zsy(ji,jj+jp,jk,1-jp,kp) 756 757 zvolf = 0.25_wp*e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) 758 759 ftv(ji,jj,jk)= ftv(ji,jj,jk)+zslope*zdzt*zvolf*ze2vr 760 fsv(ji,jj,jk)= fsv(ji,jj,jk)+zslope*zdzs*zvolf*ze2vr 761 ftvd(ji,jj,jk)=ftvd(ji,jj,jk)+fsahtv(ji,jj,jk)*zdyt*zvolf*ze2vr 762 fsvd(ji,jj,jk)=fsvd(ji,jj,jk)+fsahtv(ji,jj,jk)*zdys*zvolf*ze2vr 763 tfw(ji,jj+jp,jk+kp)=tfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zslope*zdyt 764 sfw(ji,jj+jp,jk+kp)=sfw(ji,jj+jp,jk+kp)+(zvolf*ze3wr)*zslope*zdys 765 wslp2(ji,jj+jp,jk+kp)=wslp2(ji,jj+jp,jk+kp)+ & 766 & ((zvolf*ze3wr)/(e1t(ji,jj+jp)*e2t(ji,jj+jp)))*(zslope)**2 767 psiy_eiv(ji,jj,jk+kp) = psiy_eiv(ji,jj,jk+kp) + 0.25_wp*zslope 768 769 END DO 770 END DO 771 END DO 772 END DO 773 END DO 774 775 tfw(:,:,1)=0.e0 776 sfw(:,:,1)=0.e0 777 wslp2(:,:,1)=0.e0 778 779 CALL lbc_lnk( wslp2, 'W', 1. ) 780 CALL lbc_lnk( tfw, 'W', 1. ) 781 CALL lbc_lnk( sfw, 'W', 1. ) 782 CALL lbc_lnk( ftu, 'U', -1. ) 783 CALL lbc_lnk( fsu, 'U', -1. ) 784 CALL lbc_lnk( ftv, 'V', -1. ) 785 CALL lbc_lnk( fsv, 'V', -1. ) 786 CALL lbc_lnk( ftud, 'U', -1. ) 787 CALL lbc_lnk( fsud, 'U', -1. ) 788 CALL lbc_lnk( ftvd, 'V', -1. ) 789 CALL lbc_lnk( fsvd, 'V', -1. ) 790 CALL lbc_lnk( psix_eiv, 'U', -1. ) 791 CALL lbc_lnk( psiy_eiv, 'V', -1. ) 414 iku = MAX( MIN( mbathy(ji,jj), mbathy(ji+1,jj ) ) - 1, 2 ) ! last ocean level 415 ikv = MAX( MIN( mbathy(ji,jj), mbathy(ji,jj+1 ) ) - 1, 2 ) 416 zdit (ji,jj,iku) = gtsu(ji,jj,jp_tem) ! i-gradient of T and S 417 zdis (ji,jj,iku) = gtsu(ji,jj,jp_sal) 418 zdjt (ji,jj,ikv) = gtsv(ji,jj,jp_tem) ! j-gradient of T and S 419 zdjs (ji,jj,ikv) = gtsv(ji,jj,jp_sal) 420 END DO 421 END DO 422 ENDIF 423 ! 424 zdkt(:,:,1) = 0._wp !== before vertical T & S gradient at w-level ==! 425 zdks(:,:,1) = 0._wp 426 DO jk = 2, jpk 427 zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk) 428 zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk) 429 END DO 430 ! 431 ! 432 DO jl = 0, 1 !== density i-, j-, and k-gradients ==! 433 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln) 434 DO jk = 1, jpkm1 ! done each pair of triad 435 DO jj = 1, jpjm1 ! NB: not masked due to the minimum value set 436 DO ji = 1, fs_jpim1 ! vector opt. 437 zdxrho_raw = ( zalpha(ji+ip,jj ,jk) * zdit(ji,jj,jk) + zbeta(ji+ip,jj ,jk) * zdis(ji,jj,jk) ) / e1u(ji,jj) 438 zdyrho_raw = ( zalpha(ji ,jj+jp,jk) * zdjt(ji,jj,jk) + zbeta(ji ,jj+jp,jk) * zdjs(ji,jj,jk) ) / e2v(ji,jj) 439 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 440 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 441 END DO 442 END DO 443 END DO 444 END DO 445 DO kp = 0, 1 !== density i-, j-, and k-gradients ==! 446 DO jk = 1, jpkm1 ! done each pair of triad 447 DO jj = 1, jpj ! NB: not masked due to the minimum value set 448 DO ji = 1, jpi ! vector opt. 449 zdzrho_raw = ( zalpha(ji,jj,jk) * zdkt(ji,jj,jk+kp) + zbeta(ji,jj,jk) * zdks(ji,jj,jk+kp) ) & 450 & / fse3w(ji,jj,jk+kp) 451 zdzrho(ji ,jj ,jk, kp) = - MIN( - repsln, zdzrho_raw ) ! force zdzrho >= repsln 452 END DO 453 END DO 454 END DO 455 END DO 456 ! 457 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==! 458 DO ji = 1, jpi ! i.e. 1 / (hmld+e3t(nmln)) where hmld=depw(nmln) 459 jk = MIN( nmln(ji,jj), mbathy(ji,jj) - 1 ) + 1 ! MIN in case ML depth is the ocean depth 460 z1_mlbw(ji,jj) = 1._wp / fsdepw(ji,jj,jk) 461 END DO 462 END DO 463 ! 464 ! !== intialisations to zero ==! 465 ! 466 wslp2 (:,:,:) = 0._wp ! wslp2 will be cumulated 3D field set to zero 467 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 468 triadj_g(:,:,1,:,:) = 0._wp ; triadj_g(:,:,jpk,:,:) = 0._wp 469 !!gm _iso set to zero missing 470 triadi (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 471 triadj (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp 472 473 !-------------------------------------! 474 ! Triads just below the Mixed Layer ! 475 !-------------------------------------! 476 ! 477 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 478 DO kp = 0, 1 ! with only the slope-max limit and MASKED 479 DO jj = 1, jpjm1 480 DO ji = 1, fs_jpim1 481 ip = jl ; jp = jl 482 jk = MIN( nmln(ji+ip,jj), mbathy(ji+ip,jj) - 1 ) + 1 ! ML level+1 (MIN in case ML depth is the ocean depth) 483 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 484 & + ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj) ) * umask(ji,jj,jk) 485 jk = MIN( nmln(ji,jj+jp), mbathy(ji,jj+jp) - 1 ) + 1 486 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 487 & + ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 488 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 489 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 490 END DO 491 END DO 492 END DO 493 END DO 494 495 !-------------------------------------! 496 ! Triads with surface limits ! 497 !-------------------------------------! 498 ! 499 DO kp = 0, 1 ! k-index of triads 500 DO jl = 0, 1 501 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes) 502 DO jk = 1, jpkm1 503 DO jj = 1, jpjm1 504 DO ji = 1, fs_jpim1 ! vector opt. 505 ! 506 ! Calculate slope relative to geopotentials used for GM skew fluxes 507 ! For s-coordinate, subtract slope at t-points (equivalent to *adding* gradient of depth) 508 ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 509 ! masked by umask taken at the level of dz(rho) 510 ! 511 ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti) 512 ! 513 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 514 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 515 zti_coord = ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 516 ztj_coord = ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 517 ! unmasked 518 zti_g_raw = zti_raw + zti_coord ! ref to geopot surfaces 519 ztj_g_raw = ztj_raw + ztj_coord 520 zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 521 ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 522 ! 523 ! Below ML use limited zti_g as is 524 ! Inside ML replace by linearly reducing sx_mlb towards surface 525 ! 526 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 527 zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1 528 ! ! otherwise zfact=0 529 zti_g_lim = zfacti * zti_g_lim & 530 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 531 & * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) 532 ztj_g_lim = zfactj * ztj_g_lim & 533 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 534 & * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ! masked 535 ! 536 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim * umask(ji,jj,jk+kp) 537 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim * vmask(ji,jj,jk+kp) 538 ! 539 ! Get coefficients of isoneutral diffusion tensor 540 ! 1. Utilise gradients *relative* to s-coordinate, so add t-point slopes (*subtract* depth gradients) 541 ! 2. We require that isoneutral diffusion gives no vertical buoyancy flux 542 ! i.e. 33 term = (real slope* 31, 13 terms) 543 ! To do this, retain limited sx**2 in vertical flux, but divide by real slope for 13/31 terms 544 ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 545 ! 546 zti_lim = zti_g_lim - zti_coord ! remove the coordinate slope ==> relative to coordinate surfaces 547 ztj_lim = ztj_g_lim - ztj_coord 548 zti_lim2 = zti_lim * zti_lim * umask(ji,jj,jk+kp) ! square of limited slopes ! masked <<== 549 ztj_lim2 = ztj_lim * ztj_lim * vmask(ji,jj,jk+kp) 550 ! 551 zbu = e1u(ji ,jj) * e2u(ji ,jj) * fse3u(ji ,jj,jk ) 552 zbv = e1v(ji ,jj) * e2v(ji ,jj) * fse3v(ji ,jj,jk ) 553 zbti = e1t(ji+ip,jj) * e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp) 554 zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 555 ! 556 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim2 / zti_raw ! masked 557 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim2 / ztj_raw 558 ! 559 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked 560 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2 ! masked 561 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2 562 END DO 563 END DO 564 END DO 565 END DO 566 END DO 567 ! 568 wslp2(:,:,1) = 0._wp ! force the surface wslp to zero 569 570 CALL lbc_lnk( wslp2, 'W', 1. ) ! lateral boundary confition on wslp2 only ==>>> gm : necessary ? to be checked 792 571 ! 793 572 END SUBROUTINE ldf_slp_grif … … 1008 787 !!---------------------------------------------------------------------- 1009 788 INTEGER :: ji, jj, jk ! dummy loop indices 789 INTEGER :: ierr ! local integer 1010 790 !!---------------------------------------------------------------------- 1011 791 1012 792 IF(lwp) THEN 1013 793 WRITE(numout,*) 1014 WRITE(numout,*) 'ldf_slp : direction of lateral mixing'1015 WRITE(numout,*) '~~~~~~~ '794 WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing' 795 WRITE(numout,*) '~~~~~~~~~~~~' 1016 796 ENDIF 1017 1018 ! Direction of lateral diffusion (tracers and/or momentum) 1019 ! ------------------------------ 1020 ! set the slope to zero (even in s-coordinates) 1021 1022 uslp (:,:,:) = 0.e0 1023 vslp (:,:,:) = 0.e0 1024 wslpi(:,:,:) = 0.e0 1025 wslpj(:,:,:) = 0.e0 1026 1027 uslpml (:,:) = 0.e0 1028 vslpml (:,:) = 0.e0 1029 wslpiml(:,:) = 0.e0 1030 wslpjml(:,:) = 0.e0 1031 1032 IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 1033 IF(lwp) THEN 1034 WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 797 798 IF( ln_traldf_grif ) THEN ! Griffies operator : triad of slopes 799 ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , wslp2(jpi,jpj,jpk) , STAT=ierr ) 800 ALLOCATE( triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , STAT=ierr ) 801 IF( ierr > 0 ) THEN 802 CALL ctl_stop( 'ldf_slp_init : unable to allocate Griffies operator slope ' ) ; RETURN 1035 803 ENDIF 1036 1037 ! geopotential diffusion in s-coordinates on tracers and/or momentum 1038 ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 1039 ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 1040 1041 ! set the slope of diffusion to the slope of s-surfaces 1042 ! ( c a u t i o n : minus sign as fsdep has positive value ) 1043 DO jk = 1, jpk 1044 DO jj = 2, jpjm1 1045 DO ji = fs_2, fs_jpim1 ! vector opt. 1046 uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk) 1047 vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk) 1048 wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 1049 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 804 ! 805 IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 806 ! 807 IF( ( ln_traldf_hor .AND. ln_dynldf_hor ) .AND. ln_sco ) & 808 & CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator ', & 809 & 'in s-coordinate not supported' ) 810 ! 811 ELSE ! Madec operator : slopes at u-, v-, and w-points 812 ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) , & 813 & omlmask(jpi,jpj,jpk) , uslpml(jpi,jpj) , vslpml(jpi,jpj) , wslpiml(jpi,jpj) , wslpjml(jpi,jpj) , STAT=ierr ) 814 IF( ierr > 0 ) THEN 815 CALL ctl_stop( 'ldf_slp_init : unable to allocate Madec operator slope ' ) ; RETURN 816 ENDIF 817 818 ! Direction of lateral diffusion (tracers and/or momentum) 819 ! ------------------------------ 820 uslp (:,:,:) = 0._wp ; uslpml (:,:) = 0._wp ! set the slope to zero (even in s-coordinates) 821 vslp (:,:,:) = 0._wp ; vslpml (:,:) = 0._wp 822 wslpi(:,:,:) = 0._wp ; wslpiml(:,:) = 0._wp 823 wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp 824 825 !!gm I no longer understand this..... 826 IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 827 IF(lwp) THEN 828 WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 829 ENDIF 830 831 ! geopotential diffusion in s-coordinates on tracers and/or momentum 832 ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 833 ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 834 835 ! set the slope of diffusion to the slope of s-surfaces 836 ! ( c a u t i o n : minus sign as fsdep has positive value ) 837 DO jk = 1, jpk 838 DO jj = 2, jpjm1 839 DO ji = fs_2, fs_jpim1 ! vector opt. 840 uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk) 841 vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk) 842 wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 843 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 844 END DO 1050 845 END DO 1051 846 END DO 1052 END DO 1053 ! Lateral boundary conditions on the slopes 1054 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 1055 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 1056 ENDIF 1057 ! 847 ! Lateral boundary conditions on the slopes 848 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 849 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 850 ENDIF 851 ENDIF ! 1058 852 END SUBROUTINE ldf_slp_init 1059 853
Note: See TracChangeset
for help on using the changeset viewer.
