- Timestamp:
- 2011-11-10T18:23:33+01:00 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r2980 r3077 410 410 INTEGER :: iku, ikv ! local integer 411 411 REAL(wp) :: zfacti, zfactj ! local scalars 412 REAL(wp) :: znot_thru_surface ! local scalars 412 413 REAL(wp) :: zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj 413 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_ lim2, zti_g_raw, zti_g_lim414 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_ lim2, ztj_g_raw, ztj_g_lim414 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_g_raw, zti_g_lim 415 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim 415 416 REAL(wp) :: zdzrho_raw 416 417 REAL(wp) :: zbeta0 … … 531 532 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes) 532 533 DO jk = 1, jpkm1 534 ! Must mask contribution to slope from dz/dx at constant s for triads jk=1,kp=0 that poke up though ocean surface 535 znot_thru_surface = REAL( 1-1/(jk+kp), wp ) !jk+kp=1,=0.; otherwise=1.0 533 536 DO jj = 1, jpjm1 534 537 DO ji = 1, fs_jpim1 ! vector opt. … … 543 546 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 544 547 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 545 zti_coord = ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 546 ztj_coord = ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) ! unmasked 548 549 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 550 zti_coord = znot_thru_surface * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 551 ztj_coord = znot_thru_surface * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) ! unmasked 547 552 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 548 553 ztj_g_raw = ztj_raw - ztj_coord … … 550 555 ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 551 556 ! 552 ! Below ML use limited zti_g as is 553 ! Inside ML replace by linearly reducing sx_mlb towards surface 557 ! Below ML use limited zti_g as is & mask 558 ! Inside ML replace by linearly reducing sx_mlb towards surface & mask 554 559 ! 555 560 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 556 561 zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1 557 562 ! ! otherwise zfact=0 558 zti_g_lim = zfacti * zti_g_lim &563 zti_g_lim = ( zfacti * zti_g_lim & 559 564 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 560 & * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) 561 ztj_g_lim = zfactj * ztj_g_lim &565 & * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 566 ztj_g_lim = ( zfactj * ztj_g_lim & 562 567 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 563 & * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) 564 ! 565 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim * umask(ji,jj,jk+kp) ! masked566 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim * vmask(ji,jj,jk+kp)568 & * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 569 ! 570 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim 571 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim 567 572 ! 568 573 ! Get coefficients of isoneutral diffusion tensor … … 575 580 zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces 576 581 ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 577 zti_lim2 = zti_lim * zti_lim ! square of limited slopes ! masked <<==578 ztj_lim2 = ztj_lim * ztj_lim579 582 ! 580 583 IF( ln_triad_iso ) THEN 581 zti_raw = zti_lim2/ zti_raw582 ztj_raw = ztj_lim2/ ztj_raw584 zti_raw = ( zti_lim*zti_lim ) / zti_raw 585 ztj_raw = ( ztj_lim*ztj_lim ) / ztj_raw 583 586 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 584 587 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) … … 597 600 ! 598 601 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked 599 ! agn may need to change to using ztj_g_lim**2, as wslp2 just used for eddy growth rate, needs *real* slope 600 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2 ! masked 601 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2 602 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * ( zti_g_lim * zti_g_lim ) ! masked 603 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ( ztj_g_lim * ztj_g_lim ) 602 604 END DO 603 605 END DO
Note: See TracChangeset
for help on using the changeset viewer.