- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/LDF/ldftra.F90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/LDF/ldftra.F90
r10425 r13463 94 94 95 95 !! * Substitutions 96 # include "vectopt_loop_substitute.h90" 96 # include "do_loop_substitute.h90" 97 # include "domzgr_substitute.h90" 97 98 !!---------------------------------------------------------------------- 98 99 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 152 153 ! ================================= 153 154 ! 154 REWIND( numnam_ref ) ! Namelist namtra_ldf in reference namelist : Lateral physics on tracers155 155 READ ( numnam_ref, namtra_ldf, IOSTAT = ios, ERR = 901) 156 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldf in reference namelist', lwp ) 157 REWIND( numnam_cfg ) ! Namelist namtra_ldf in configuration namelist : Lateral physics on tracers 156 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldf in reference namelist' ) 158 157 READ ( numnam_cfg, namtra_ldf, IOSTAT = ios, ERR = 902 ) 159 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_ldf in configuration namelist' , lwp)158 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_ldf in configuration namelist' ) 160 159 IF(lwm) WRITE( numond, namtra_ldf ) 161 160 ! … … 318 317 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F(i,j) read in eddy_diffusivity.nc file' 319 318 CALL iom_open( 'eddy_diffusivity_2D.nc', inum ) 320 CALL iom_get ( inum, jpdom_ data, 'ahtu_2D', ahtu(:,:,1))321 CALL iom_get ( inum, jpdom_ data, 'ahtv_2D', ahtv(:,:,1))319 CALL iom_get ( inum, jpdom_global, 'ahtu_2D', ahtu(:,:,1), cd_type = 'U', psgn = 1._wp ) 320 CALL iom_get ( inum, jpdom_global, 'ahtv_2D', ahtv(:,:,1), cd_type = 'V', psgn = 1._wp ) 322 321 CALL iom_close( inum ) 323 322 DO jk = 2, jpkm1 … … 346 345 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F(i,j,k) read in eddy_diffusivity.nc file' 347 346 CALL iom_open( 'eddy_diffusivity_3D.nc', inum ) 348 CALL iom_get ( inum, jpdom_ data, 'ahtu_3D', ahtu)349 CALL iom_get ( inum, jpdom_ data, 'ahtv_3D', ahtv)347 CALL iom_get ( inum, jpdom_global, 'ahtu_3D', ahtu, cd_type = 'U', psgn = 1._wp ) 348 CALL iom_get ( inum, jpdom_global, 'ahtv_3D', ahtv, cd_type = 'V', psgn = 1._wp ) 350 349 CALL iom_close( inum ) 351 350 ! … … 382 381 383 382 384 SUBROUTINE ldf_tra( kt )383 SUBROUTINE ldf_tra( kt, Kbb, Kmm ) 385 384 !!---------------------------------------------------------------------- 386 385 !! *** ROUTINE ldf_tra *** … … 403 402 !!---------------------------------------------------------------------- 404 403 INTEGER, INTENT(in) :: kt ! time step 404 INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices 405 405 ! 406 406 INTEGER :: ji, jj, jk ! dummy loop indices … … 411 411 ! ! =F(growth rate of baroclinic instability) 412 412 ! ! max value aeiv_0 ; decreased to 0 within 20N-20S 413 CALL ldf_eiv( kt, aei0, aeiu, aeiv )413 CALL ldf_eiv( kt, aei0, aeiu, aeiv, Kmm ) 414 414 ENDIF 415 415 ! … … 424 424 ahtv(:,:,1) = aeiv(:,:,1) 425 425 ELSE ! compute aht. 426 CALL ldf_eiv( kt, aht0, ahtu, ahtv )426 CALL ldf_eiv( kt, aht0, ahtu, ahtv, Kmm ) 427 427 ENDIF 428 428 ! … … 430 430 zaht_min = 0.2_wp * aht0 ! minimum value for aht 431 431 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 432 DO_2D( 1, 1, 1, 1 ) 433 !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg) 434 !! ==>>> The Coriolis value is identical for t- & u_points, and for v- and f-points 435 zaht = ( 1._wp - MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * zDaht 436 zahf = ( 1._wp - MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * zDaht 437 ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) ! min value zaht_min 438 ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zahf ) ! increase within 20S-20N 439 END_2D 442 440 DO jk = 1, jpkm1 ! deeper value = surface value + mask for all levels 443 441 ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk) … … 448 446 IF( ln_traldf_lap ) THEN ! laplacian operator |u| e /12 449 447 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_12448 ahtu(:,:,jk) = ABS( uu(:,:,jk,Kbb) ) * e1u(:,:) * r1_12 ! n.b. uu,vv are masked 449 ahtv(:,:,jk) = ABS( vv(:,:,jk,Kbb) ) * e2v(:,:) * r1_12 452 450 END DO 453 451 ELSEIF( ln_traldf_blp ) THEN ! bilaplacian operator sqrt( |u| e^3 /12 ) = sqrt( |u| e /12 ) * e 454 452 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(:,:)453 ahtu(:,:,jk) = SQRT( ABS( uu(:,:,jk,Kbb) ) * e1u(:,:) * r1_12 ) * e1u(:,:) 454 ahtv(:,:,jk) = SQRT( ABS( vv(:,:,jk,Kbb) ) * e2v(:,:) * r1_12 ) * e2v(:,:) 457 455 END DO 458 456 ENDIF … … 510 508 ENDIF 511 509 ! 512 REWIND( numnam_ref ) ! Namelist namtra_eiv in reference namelist : eddy induced velocity param.513 510 READ ( numnam_ref, namtra_eiv, IOSTAT = ios, ERR = 901) 514 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_eiv in reference namelist', lwp ) 515 ! 516 REWIND( numnam_cfg ) ! Namelist namtra_eiv in configuration namelist : eddy induced velocity param. 511 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_eiv in reference namelist' ) 512 ! 517 513 READ ( numnam_cfg, namtra_eiv, IOSTAT = ios, ERR = 902 ) 518 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_eiv in configuration namelist' , lwp)514 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_eiv in configuration namelist' ) 519 515 IF(lwm) WRITE ( numond, namtra_eiv ) 520 516 … … 576 572 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F(i,j) read in eddy_diffusivity_2D.nc file' 577 573 CALL iom_open ( 'eddy_induced_velocity_2D.nc', inum ) 578 CALL iom_get ( inum, jpdom_ data, 'aeiu', aeiu(:,:,1))579 CALL iom_get ( inum, jpdom_ data, 'aeiv', aeiv(:,:,1))574 CALL iom_get ( inum, jpdom_global, 'aeiu', aeiu(:,:,1), cd_type = 'U', psgn = 1._wp ) 575 CALL iom_get ( inum, jpdom_global, 'aeiv', aeiv(:,:,1), cd_type = 'V', psgn = 1._wp ) 580 576 CALL iom_close( inum ) 581 577 DO jk = 2, jpkm1 … … 600 596 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file' 601 597 CALL iom_open ( 'eddy_induced_velocity_3D.nc', inum ) 602 CALL iom_get ( inum, jpdom_ data, 'aeiu', aeiu)603 CALL iom_get ( inum, jpdom_ data, 'aeiv', aeiv)598 CALL iom_get ( inum, jpdom_global, 'aeiu', aeiu, cd_type = 'U', psgn = 1._wp ) 599 CALL iom_get ( inum, jpdom_global, 'aeiv', aeiv, cd_type = 'V', psgn = 1._wp ) 604 600 CALL iom_close( inum ) 605 601 ! … … 625 621 626 622 627 SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv )623 SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv, Kmm ) 628 624 !!---------------------------------------------------------------------- 629 625 !! *** ROUTINE ldf_eiv *** … … 637 633 !!---------------------------------------------------------------------- 638 634 INTEGER , INTENT(in ) :: kt ! ocean time-step index 635 INTEGER , INTENT(in ) :: Kmm ! ocean time level indices 639 636 REAL(wp) , INTENT(inout) :: paei0 ! max value [m2/s] 640 637 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: paeiu, paeiv ! eiv coefficient [m2/s] … … 651 648 ! ! Compute lateral diffusive coefficient at T-point 652 649 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 650 DO_3D( 0, 0, 0, 0, 1, jpk ) 651 ! Take the max of N^2 and zero then take the vertical sum 652 ! of the square root of the resulting N^2 ( required to compute 653 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 654 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 655 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm) 656 ! Compute elements required for the inverse time scale of baroclinic 657 ! eddies using the isopycnal slopes calculated in ldfslp.F : 658 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 659 ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 660 zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w 661 zhw(ji,jj) = zhw(ji,jj) + ze3w 662 END_3D 670 663 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 664 DO_3D( 0, 0, 0, 0, 1, jpk ) 665 ! Take the max of N^2 and zero then take the vertical sum 666 ! of the square root of the resulting N^2 ( required to compute 667 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 668 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 669 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm) 670 ! Compute elements required for the inverse time scale of baroclinic 671 ! eddies using the isopycnal slopes calculated in ldfslp.F : 672 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 673 ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) 674 zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 675 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w 676 zhw(ji,jj) = zhw(ji,jj) + ze3w 677 END_3D 678 ENDIF 679 680 DO_2D( 0, 0, 0, 0 ) 681 zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 682 ! Rossby radius at w-point taken betwenn 2 km and 40km 683 zRo(ji,jj) = MAX( 2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 ) ) 684 ! Compute aeiw by multiplying Ro^2 and T^-1 685 zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 686 END_2D 700 687 701 688 ! !== Bound on eiv coeff. ==! 702 689 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 709 CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1. ) ! lateral boundary condition 690 DO_2D( 0, 0, 0, 0 ) 691 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease 692 zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0 693 END_2D 694 CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1.0_wp ) ! lateral boundary condition 710 695 ! 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 717 CALL lbc_lnk_multi( 'ldftra', paeiu(:,:,1), 'U', 1. , paeiv(:,:,1), 'V', 1. ) ! lateral boundary condition 696 DO_2D( 0, 0, 0, 0 ) 697 paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj ) ) * umask(ji,jj,1) 698 paeiv(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji ,jj+1) ) * vmask(ji,jj,1) 699 END_2D 700 CALL lbc_lnk_multi( 'ldftra', paeiu(:,:,1), 'U', 1.0_wp , paeiv(:,:,1), 'V', 1.0_wp ) ! lateral boundary condition 718 701 719 702 DO jk = 2, jpkm1 !== deeper values equal the surface one ==! … … 725 708 726 709 727 SUBROUTINE ldf_eiv_trp( kt, kit000, pu n, pvn, pwn, cdtype)710 SUBROUTINE ldf_eiv_trp( kt, kit000, pu, pv, pw, cdtype, Kmm, Krhs ) 728 711 !!---------------------------------------------------------------------- 729 712 !! *** ROUTINE ldf_eiv_trp *** … … 741 724 !! velocity and heat transport (call ldf_eiv_dia) 742 725 !! 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] 726 !! ** Action : pu, pv increased by the eiv transport 727 !!---------------------------------------------------------------------- 728 INTEGER , INTENT(in ) :: kt ! ocean time-step index 729 INTEGER , INTENT(in ) :: kit000 ! first time step index 730 INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices 731 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 732 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components [m3/s] 733 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pv ! out: 3 ocean transport components [m3/s] 734 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the eiv [m3/s] 751 735 !! 752 736 INTEGER :: ji, jj, jk ! dummy loop indices … … 766 750 zpsi_uw(:,:,jpk) = 0._wp ; zpsi_vw(:,:,jpk) = 0._wp 767 751 ! 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 752 DO_3D( 1, 0, 1, 0, 2, jpkm1 ) 753 zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk ) + wslpi(ji+1,jj,jk) ) & 754 & * ( aeiu (ji,jj,jk-1) + aeiu (ji ,jj,jk) ) * wumask(ji,jj,jk) 755 zpsi_vw(ji,jj,jk) = - r1_4 * e1v(ji,jj) * ( wslpj(ji,jj,jk ) + wslpj(ji,jj+1,jk) ) & 756 & * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj ,jk) ) * wvmask(ji,jj,jk) 757 END_3D 758 ! 759 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 760 pu(ji,jj,jk) = pu(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 761 pv(ji,jj,jk) = pv(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 762 END_3D 763 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 764 pw(ji,jj,jk) = pw(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj ,jk) & 765 & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji ,jj-1,jk) ) 766 END_3D 795 767 ! 796 768 ! ! 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 )769 IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 798 770 ! 799 771 END SUBROUTINE ldf_eiv_trp 800 772 801 773 802 SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw )774 SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw, Kmm ) 803 775 !!---------------------------------------------------------------------- 804 776 !! *** ROUTINE ldf_eiv_dia *** … … 811 783 !!---------------------------------------------------------------------- 812 784 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: psi_uw, psi_vw ! streamfunction [m3/s] 785 INTEGER , INTENT(in ) :: Kmm ! ocean time level indices 813 786 ! 814 787 INTEGER :: ji, jj, jk ! dummy loop indices … … 821 794 !!gm to be redesigned.... 822 795 ! !== eiv stream function: output ==! 823 CALL lbc_lnk_multi( 'ldftra', psi_uw, 'U', -1. , psi_vw, 'V', -1.)796 CALL lbc_lnk_multi( 'ldftra', psi_uw, 'U', -1.0_wp , psi_vw, 'V', -1.0_wp ) 824 797 ! 825 798 !!gm CALL iom_put( "psi_eiv_uw", psi_uw ) ! output … … 831 804 ! 832 805 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) )806 zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u(:,:,jk,Kmm) ) 834 807 END DO 835 808 CALL iom_put( "uoce_eiv", zw3d ) 836 809 ! 837 810 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) )811 zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v(:,:,jk,Kmm) ) 839 812 END DO 840 813 CALL iom_put( "voce_eiv", zw3d ) 841 814 ! 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 815 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 816 zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk) - psi_vw(ji ,jj-1,jk) & 817 & + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj ,jk) ) / e1e2t(ji,jj) 818 END_3D 819 CALL lbc_lnk( 'ldftra', zw3d, 'T', 1.0_wp ) ! lateral boundary condition 820 CALL iom_put( "woce_eiv", zw3d ) 821 ! 822 IF( iom_use('weiv_masstr') ) THEN ! vertical mass transport & its square value 823 zw2d(:,:) = rho0 * e1e2t(:,:) 824 DO jk = 1, jpk 825 zw3d(:,:,jk) = zw3d(:,:,jk) * zw2d(:,:) 848 826 END DO 849 END DO 850 CALL lbc_lnk( 'ldftra', zw3d, 'T', 1. ) ! lateral boundary condition 851 CALL iom_put( "woce_eiv", zw3d ) 852 ! 853 ! 854 zztmp = 0.5_wp * rau0 * rcp 827 CALL iom_put( "weiv_masstr" , zw3d ) 828 ENDIF 829 ! 830 IF( iom_use('ueiv_masstr') ) THEN 831 zw3d(:,:,:) = 0.e0 832 DO jk = 1, jpkm1 833 zw3d(:,:,jk) = rho0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) 834 END DO 835 CALL iom_put( "ueiv_masstr", zw3d ) ! mass transport in i-direction 836 ENDIF 837 ! 838 zztmp = 0.5_wp * rho0 * rcp 855 839 IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN 856 840 zw2d(:,:) = 0._wp 857 841 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 867 CALL lbc_lnk( 'ldftra', zw2d, 'U', -1. ) 868 CALL lbc_lnk( 'ldftra', zw3d, 'U', -1. ) 842 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 843 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1) - psi_uw(ji ,jj,jk) ) & 844 & * ( ts (ji,jj,jk,jp_tem,Kmm) + ts (ji+1,jj,jk,jp_tem,Kmm) ) 845 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 846 END_3D 847 CALL lbc_lnk( 'ldftra', zw2d, 'U', -1.0_wp ) 848 CALL lbc_lnk( 'ldftra', zw3d, 'U', -1.0_wp ) 869 849 CALL iom_put( "ueiv_heattr" , zztmp * zw2d ) ! heat transport in i-direction 870 850 CALL iom_put( "ueiv_heattr3d", zztmp * zw3d ) ! heat transport in i-direction 871 851 ENDIF 852 ! 853 IF( iom_use('veiv_masstr') ) THEN 854 zw3d(:,:,:) = 0.e0 855 DO jk = 1, jpkm1 856 zw3d(:,:,jk) = rho0 * ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) 857 END DO 858 CALL iom_put( "veiv_masstr", zw3d ) ! mass transport in i-direction 859 ENDIF 860 ! 872 861 zw2d(:,:) = 0._wp 873 862 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 883 CALL lbc_lnk( 'ldftra', zw2d, 'V', -1. ) 863 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 864 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj ,jk) ) & 865 & * ( ts (ji,jj,jk,jp_tem,Kmm) + ts (ji,jj+1,jk,jp_tem,Kmm) ) 866 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 867 END_3D 868 CALL lbc_lnk( 'ldftra', zw2d, 'V', -1.0_wp ) 884 869 CALL iom_put( "veiv_heattr", zztmp * zw2d ) ! heat transport in j-direction 885 870 CALL iom_put( "veiv_heattr", zztmp * zw3d ) ! heat transport in j-direction 886 871 ! 887 IF( ln_diaptr )CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d )872 IF( iom_use( 'sophteiv' ) ) CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d ) 888 873 ! 889 874 zztmp = 0.5_wp * 0.5 … … 891 876 zw2d(:,:) = 0._wp 892 877 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 902 CALL lbc_lnk( 'ldftra', zw2d, 'U', -1. ) 903 CALL lbc_lnk( 'ldftra', zw3d, 'U', -1. ) 878 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 879 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1) - psi_uw(ji ,jj,jk) ) & 880 & * ( ts (ji,jj,jk,jp_sal,Kmm) + ts (ji+1,jj,jk,jp_sal,Kmm) ) 881 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 882 END_3D 883 CALL lbc_lnk( 'ldftra', zw2d, 'U', -1.0_wp ) 884 CALL lbc_lnk( 'ldftra', zw3d, 'U', -1.0_wp ) 904 885 CALL iom_put( "ueiv_salttr", zztmp * zw2d ) ! salt transport in i-direction 905 CALL iom_put( "ueiv_salttr3d", zztmp * zw3d ) 886 CALL iom_put( "ueiv_salttr3d", zztmp * zw3d ) ! salt transport in i-direction 906 887 ENDIF 907 888 zw2d(:,:) = 0._wp 908 889 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 918 CALL lbc_lnk( 'ldftra', zw2d, 'V', -1. ) 890 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 891 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj ,jk) ) & 892 & * ( ts (ji,jj,jk,jp_sal,Kmm) + ts (ji,jj+1,jk,jp_sal,Kmm) ) 893 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 894 END_3D 895 CALL lbc_lnk( 'ldftra', zw2d, 'V', -1.0_wp ) 919 896 CALL iom_put( "veiv_salttr", zztmp * zw2d ) ! salt transport in j-direction 920 897 CALL iom_put( "veiv_salttr", zztmp * zw3d ) ! salt transport in j-direction 921 898 ! 922 IF( ln_diaptr) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )899 IF( iom_use( 'sopsteiv' ) ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d ) 923 900 ! 924 901 !
Note: See TracChangeset
for help on using the changeset viewer.
