- Timestamp:
- 2011-10-20T18:00:23+02:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r2840 r2969 402 402 !!---------------------------------------------------------------------- 403 403 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 404 USE oce , ONLY: zalbet => ua ! use ua as workspace404 USE oce , ONLY: zalbet => ua ! use ua as workspace 405 405 USE wrk_nemo, ONLY: z1_mlbw => wrk_2d_1 406 406 !! 407 INTEGER, INTENT( in ) :: kt ! ocean time-step index407 INTEGER, INTENT( in ) :: kt ! ocean time-step index 408 408 !! 409 409 INTEGER :: ji, jj, jk, jl, ip, jp, kp ! dummy loop indices … … 425 425 !--------------------------------! 426 426 ! 427 CALL eos_alpbet( tsb, zalbet, zbeta0 ) 428 ! 429 DO jl = 0, 1 !== unmasked before density i- j-, k-gradients ==!427 CALL eos_alpbet( tsb, zalbet, zbeta0 ) !== before local thermal/haline expension ratio at T-points ==! 428 ! 429 DO jl = 0, 1 !== unmasked before density i- j-, k-gradients ==! 430 430 ! 431 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln)432 DO jk = 1, jpkm1 433 DO jj = 1, jpjm1 434 DO ji = 1, fs_jpim1 ! vector opt.431 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln) 432 DO jk = 1, jpkm1 ! done each pair of triad 433 DO jj = 1, jpjm1 ! NB: not masked ==> a minimum value is set 434 DO ji = 1, fs_jpim1 ! vector opt. 435 435 zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! i-gradient of T & S at u-point 436 436 zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) … … 445 445 END DO 446 446 ! 447 IF( ln_zps.and.l n_grad_zps ) THEN! partial steps: correction of i- & j-grad on bottom447 IF( ln_zps.and.l_grad_zps ) THEN ! partial steps: correction of i- & j-grad on bottom 448 448 # if defined key_vectopt_loop 449 449 DO jj = 1, 1 450 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)450 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 451 451 # else 452 452 DO jj = 1, jpjm1 … … 466 466 END DO 467 467 468 DO kp = 0, 1 !== unmasked before density i- j-, k-gradients ==!469 DO jk = 1, jpkm1 470 DO jj = 1, jpj 471 DO ji = 1, jpi ! vector opt.472 IF( jk+kp > 1 ) THEN 468 DO kp = 0, 1 !== unmasked before density i- j-, k-gradients ==! 469 DO jk = 1, jpkm1 ! done each pair of triad 470 DO jj = 1, jpj ! NB: not masked ==> a minimum value is set 471 DO ji = 1, jpi ! vector opt. 472 IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp 473 473 zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) ) 474 474 zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) ) … … 484 484 END DO 485 485 ! 486 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==!486 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==! 487 487 DO ji = 1, jpi 488 488 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth … … 491 491 END DO 492 492 ! 493 ! !== intialisations to zero ==!494 ! 495 wslp2 (:,:,:) = 0._wp 493 ! !== intialisations to zero ==! 494 ! 495 wslp2 (:,:,:) = 0._wp ! wslp2 will be cumulated 3D field set to zero 496 496 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 497 497 triadj_g(:,:,1,:,:) = 0._wp ; triadj_g(:,:,jpk,:,:) = 0._wp 498 !!gm _iso set to zero missing498 !!gm _iso set to zero missing 499 499 triadi (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 500 500 triadj (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp … … 504 504 !-------------------------------------! 505 505 ! 506 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base507 DO kp = 0, 1 ! with only the slope-max limit and MASKED506 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 507 DO kp = 0, 1 ! with only the slope-max limit and MASKED 508 508 DO jj = 1, jpjm1 509 509 DO ji = 1, fs_jpim1 … … 527 527 !-------------------------------------! 528 528 ! 529 DO kp = 0, 1 ! k-index of triads529 DO kp = 0, 1 ! k-index of triads 530 530 DO jl = 0, 1 531 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes)531 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes) 532 532 DO jk = 1, jpkm1 533 533 DO jj = 1, jpjm1 534 DO ji = 1, fs_jpim1 ! vector opt.534 DO ji = 1, fs_jpim1 ! vector opt. 535 535 ! 536 536 ! Calculate slope relative to geopotentials used for GM skew fluxes … … 575 575 zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces 576 576 ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 577 zti_lim2 = zti_lim * zti_lim ! square of limited slopes ! masked <<==577 zti_lim2 = zti_lim * zti_lim ! square of limited slopes ! masked <<== 578 578 ztj_lim2 = ztj_lim * ztj_lim 579 579 ! … … 596 596 zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 597 597 ! 598 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked599 ! agn may need to change to using ztj_g_lim**2, as wslp2 just used for eddy growth rate, needs *real* slope598 !!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 600 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2 ! masked 601 601 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2 … … 637 637 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_dzr ! z-gradient of density (T-point) 638 638 !! 639 INTEGER :: ji , jj , jk ! dummy loop indices640 INTEGER :: iku, ikv, ik, ikm1 ! local integers639 INTEGER :: ji , jj , jk ! dummy loop indices 640 INTEGER :: iku, ikv, ik, ikm1 ! local integers 641 641 REAL(wp) :: zeps, zm1_g, zm1_2g ! local scalars 642 642 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - … … 654 654 wslpjml(1,:) = 0._wp ; wslpjml(jpi,:) = 0._wp 655 655 ! 656 ! !== surface mixed layer mask !657 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise656 ! !== surface mixed layer mask ! 657 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise 658 658 # if defined key_vectopt_loop 659 659 DO jj = 1, 1 660 DO ji = 1, jpij ! vector opt. (forced unrolling)660 DO ji = 1, jpij ! vector opt. (forced unrolling) 661 661 # else 662 662 DO jj = 1, jpj … … 689 689 DO ji = 2, jpim1 690 690 # endif 691 ! !== Slope at u- & v-points just below the Mixed Layer ==!691 ! !== Slope at u- & v-points just below the Mixed Layer ==! 692 692 ! 693 ! 693 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 694 694 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 695 695 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 696 696 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 697 697 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 698 ! 698 ! !- horizontal density gradient at u- & v-points 699 699 zau = p_gru(ji,jj,iku) / e1u(ji,jj) 700 700 zav = p_grv(ji,jj,ikv) / e2v(ji,jj) 701 ! 702 ! 701 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 702 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 703 703 zbu = MIN( zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) ) 704 704 zbv = MIN( zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) ) 705 ! 705 ! !- Slope at u- & v-points (uslpml, vslpml) 706 706 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 707 707 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 708 708 ! 709 ! !== i- & j-slopes at w-points just below the Mixed Layer ==!709 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! 710 710 ! 711 711 ik = MIN( nmln(ji,jj) + 1, jpk ) 712 712 ikm1 = MAX( 1, ik-1 ) 713 ! 713 ! !- vertical density gradient for w-slope (from N^2) 714 714 zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 715 ! 715 ! !- horizontal density i- & j-gradient at w-points 716 716 zci = MAX( umask(ji-1,jj,ik ) + umask(ji,jj,ik ) & 717 717 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) … … 722 722 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & 723 723 & + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik) 724 ! 725 ! 724 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 725 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 726 726 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zai ) ) 727 727 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj ) ) 728 ! 728 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 729 729 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 730 730 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 731 731 END DO 732 732 END DO 733 !!gm this lbc_lnk should be useless....733 !!gm this lbc_lnk should be useless.... 734 734 CALL lbc_lnk( uslpml , 'U', -1. ) ; CALL lbc_lnk( vslpml , 'V', -1. ) ! lateral boundary cond. (sign change) 735 735 CALL lbc_lnk( wslpiml, 'W', -1. ) ; CALL lbc_lnk( wslpjml, 'W', -1. ) ! lateral boundary conditions … … 765 765 IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 766 766 ! 767 ! IF( ( ln_traldf_hor .OR. ln_dynldf_hor ) .AND. ln_sco ) &768 ! CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator in s-coordinate not supported' )769 !770 767 ELSE ! Madec operator : slopes at u-, v-, and w-points 771 768 ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) , & … … 780 777 wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp 781 778 782 !!gm I no longer understand this.....779 !!gm I no longer understand this..... 783 780 IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 784 781 IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' … … 813 810 LOGICAL, PUBLIC, PARAMETER :: lk_ldfslp = .FALSE. !: slopes flag 814 811 CONTAINS 815 SUBROUTINE ldf_slp( kt, prd, pn2 ) 812 SUBROUTINE ldf_slp( kt, prd, pn2 ) ! Dummy routine 816 813 INTEGER, INTENT(in) :: kt 817 814 REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2 … … 822 819 WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt 823 820 END SUBROUTINE ldf_slp_grif 824 SUBROUTINE ldf_slp_init ! Dummy routine821 SUBROUTINE ldf_slp_init ! Dummy routine 825 822 END SUBROUTINE ldf_slp_init 826 823 #endif
Note: See TracChangeset
for help on using the changeset viewer.
