- Timestamp:
- 2020-12-02T16:32:24+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/LDF
- Files:
-
- 2 edited
-
ldfc1d_c2d.F90 (modified) (8 diffs)
-
ldftra.F90 (modified) (37 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/LDF/ldfc1d_c2d.F90
r13497 r14017 2 2 !!====================================================================== 3 3 !! *** MODULE ldfc1d_c2d *** 4 !! Ocean physics: profile and horizontal shape of lateral eddy coefficients 4 !! Ocean physics: profile and horizontal shape of lateral eddy coefficients 5 5 !!===================================================================== 6 6 !! History : 3.7 ! 2013-12 (G. Madec) restructuration/simplification of aht/aeiv specification, … … 9 9 10 10 !!---------------------------------------------------------------------- 11 !! ldf_c1d : ah reduced by 1/4 on the vertical (tanh profile, inflection at 300m) 11 !! ldf_c1d : ah reduced by 1/4 on the vertical (tanh profile, inflection at 300m) 12 12 !! ldf_c2d : ah = F(e1,e2) (laplacian or = F(e1^3,e2^3) (bilaplacian) 13 13 !!---------------------------------------------------------------------- … … 29 29 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 30 30 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 31 31 32 32 !! * Substitutions 33 33 # include "do_loop_substitute.h90" … … 42 42 !!---------------------------------------------------------------------- 43 43 !! *** ROUTINE ldf_c1d *** 44 !! 44 !! 45 45 !! ** Purpose : 1D eddy diffusivity/viscosity coefficients 46 46 !! 47 47 !! ** Method : 1D eddy diffusivity coefficients F( depth ) 48 !! Reduction by zratio from surface to bottom 49 !! hyperbolic tangent profile with inflection point 48 !! Reduction by zratio from surface to bottom 49 !! hyperbolic tangent profile with inflection point 50 50 !! at zh=500m and a width of zw=200m 51 51 !! … … 95 95 END_3D 96 96 ! Lateral boundary conditions 97 CALL lbc_lnk_multi( 'ldfc1d_c2d', pah1, 'U', 1.0_wp , pah2, 'V', 1.0_wp ) 97 CALL lbc_lnk_multi( 'ldfc1d_c2d', pah1, 'U', 1.0_wp , pah2, 'V', 1.0_wp ) 98 98 ! 99 99 CASE DEFAULT ! error … … 107 107 !!---------------------------------------------------------------------- 108 108 !! *** ROUTINE ldf_c2d *** 109 !! 109 !! 110 110 !! ** Purpose : 2D eddy diffusivity/viscosity coefficients 111 111 !! … … 113 113 !! laplacian operator : ah proportional to the scale factor [m2/s] 114 114 !! bilaplacian operator : ah proportional to the (scale factor)^3 [m4/s] 115 !! In both cases, pah0 is the maximum value reached by the coefficient 115 !! In both cases, pah0 is the maximum value reached by the coefficient 116 116 !! at the Equator in case of e1=ra*rad= ~111km, not over the whole domain. 117 117 !! … … 140 140 END_2D 141 141 CASE( 'TRA' ) ! U- and V-points 142 DO_2D( 1, 1, 1, 1 ) 142 ! NOTE: [tiling-comms-merge] Change needed to preserve results with respect to the trunk 143 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 143 144 pah1(ji,jj,1) = pUfac * MAX( e1u(ji,jj), e2u(ji,jj) )**knn 144 145 pah2(ji,jj,1) = pUfac * MAX( e1v(ji,jj), e2v(ji,jj) )**knn -
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/LDF/ldftra.F90
r13558 r14017 2 2 !!====================================================================== 3 3 !! *** MODULE ldftra *** 4 !! Ocean physics: lateral diffusivity coefficients 4 !! Ocean physics: lateral diffusivity coefficients 5 5 !!===================================================================== 6 6 !! History : ! 1997-07 (G. Madec) from inimix.F split in 2 routines 7 7 !! NEMO 1.0 ! 2002-09 (G. Madec) F90: Free form and module 8 !! 2.0 ! 2005-11 (G. Madec) 8 !! 2.0 ! 2005-11 (G. Madec) 9 9 !! 3.7 ! 2013-12 (F. Lemarie, G. Madec) restructuration/simplification of aht/aeiv specification, 10 10 !! ! add velocity dependent coefficient and optional read in file … … 13 13 !!---------------------------------------------------------------------- 14 14 !! ldf_tra_init : initialization, namelist read, and parameters control 15 !! ldf_tra : update lateral eddy diffusivity coefficients at each time step 16 !! ldf_eiv_init : initialization of the eiv coeff. from namelist choices 15 !! ldf_tra : update lateral eddy diffusivity coefficients at each time step 16 !! ldf_eiv_init : initialization of the eiv coeff. from namelist choices 17 17 !! ldf_eiv : time evolution of the eiv coefficients (function of the growth rate of baroclinic instability) 18 18 !! ldf_eiv_trp : add to the input ocean transport the contribution of the EIV parametrization … … 23 23 USE phycst ! physical constants 24 24 USE ldfslp ! lateral diffusion: slope of iso-neutral surfaces 25 USE ldfc1d_c2d ! lateral diffusion: 1D & 2D cases 25 USE ldfc1d_c2d ! lateral diffusion: 1D & 2D cases 26 26 USE diaptr 27 27 ! … … 40 40 PUBLIC ldf_eiv_trp ! called by traadv.F90 41 41 PUBLIC ldf_eiv_dia ! called by traldf_iso and traldf_iso_triad.F90 42 43 ! !!* Namelist namtra_ldf : lateral mixing on tracers * 42 43 ! !!* Namelist namtra_ldf : lateral mixing on tracers * 44 44 ! != Operator type =! 45 45 LOGICAL , PUBLIC :: ln_traldf_OFF !: no operator: No explicit diffusion … … 52 52 ! != iso-neutral options =! 53 53 ! LOGICAL , PUBLIC :: ln_traldf_triad !: griffies triad scheme (see ldfslp) 54 LOGICAL , PUBLIC :: ln_traldf_msc !: Method of Stabilizing Correction 54 LOGICAL , PUBLIC :: ln_traldf_msc !: Method of Stabilizing Correction 55 55 ! LOGICAL , PUBLIC :: ln_triad_iso !: pure horizontal mixing in ML (see ldfslp) 56 56 ! LOGICAL , PUBLIC :: ln_botmix_triad !: mixing on bottom (see ldfslp) … … 59 59 ! != Coefficients =! 60 60 INTEGER , PUBLIC :: nn_aht_ijk_t !: choice of time & space variations of the lateral eddy diffusivity coef. 61 ! ! time invariant coefficients: aht_0 = 1/2 Ud*Ld (lap case) 61 ! ! time invariant coefficients: aht_0 = 1/2 Ud*Ld (lap case) 62 62 ! ! bht_0 = 1/12 Ud*Ld^3 (blp case) 63 63 REAL(wp), PUBLIC :: rn_Ud !: lateral diffusive velocity [m/s] … … 72 72 REAL(wp), PUBLIC :: rn_Ue !: lateral diffusive velocity [m/s] 73 73 REAL(wp), PUBLIC :: rn_Le !: lateral diffusive length [m] 74 74 75 75 ! ! Flag to control the type of lateral diffusive operator 76 76 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in specification of lateral diffusion … … 106 106 !!---------------------------------------------------------------------- 107 107 !! *** ROUTINE ldf_tra_init *** 108 !! 108 !! 109 109 !! ** Purpose : initializations of the tracer lateral mixing coeff. 110 110 !! … … 116 116 !! nn_aht_ijk_t = 0 => = constant 117 117 !! ! 118 !! = 10 => = F(z) : constant with a reduction of 1/4 with depth 118 !! = 10 => = F(z) : constant with a reduction of 1/4 with depth 119 119 !! ! 120 120 !! =-20 => = F(i,j) = shape read in 'eddy_diffusivity.nc' file … … 126 126 !! = 31 = F(i,j,k,t) = F(local velocity) ( 1/2 |u|e laplacian operator 127 127 !! or 1/12 |u|e^3 bilaplacian operator ) 128 !! * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init 129 !! 128 !! * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init 129 !! 130 130 !! ** action : ahtu, ahtv initialized one for all or l_ldftra_time set to true 131 131 !! aeiu, aeiv initialized one for all or l_ldfeiv_time set to true … … 148 148 WRITE(numout,*) '~~~~~~~~~~~~ ' 149 149 ENDIF 150 150 151 151 ! 152 152 ! Choice of lateral tracer physics … … 182 182 ! 183 183 ! 184 ! Operator and its acting direction (set nldf_tra) 184 ! Operator and its acting direction (set nldf_tra) 185 185 ! ================================= 186 186 ! … … 210 210 ENDIF 211 211 IF ( ln_zps ) THEN ! z-coordinate with partial step 212 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 212 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 213 213 IF ( ln_traldf_hor ) nldf_tra = np_lap ! horizontal (no rotation) 214 214 IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard (rotation) … … 231 231 ENDIF 232 232 IF ( ln_zps ) THEN ! z-coordinate with partial step 233 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 233 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 234 234 IF ( ln_traldf_hor ) nldf_tra = np_blp ! horizontal (no rotation) 235 235 IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation) … … 249 249 ! 250 250 IF( nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. & 251 & nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it ) l_ldfslp = .TRUE. ! slope of neutral surfaces required 251 & nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it ) l_ldfslp = .TRUE. ! slope of neutral surfaces required 252 252 ! 253 253 IF( ln_traldf_blp .AND. ( ln_traldf_iso .OR. ln_traldf_triad) ) THEN ! iso-neutral bilaplacian need MSC … … 270 270 271 271 ! 272 ! Space/time variation of eddy coefficients 272 ! Space/time variation of eddy coefficients 273 273 ! =========================================== 274 274 ! … … 286 286 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays') 287 287 ! 288 ahtu(:,:,jpk) = 0._wp ! last level always 0 288 ahtu(:,:,jpk) = 0._wp ! last level always 0 289 289 ahtv(:,:,jpk) = 0._wp 290 290 !. … … 363 363 END SELECT 364 364 ! 365 IF( .NOT.l_ldftra_time ) THEN !* No time variation 365 IF( .NOT.l_ldftra_time ) THEN !* No time variation 366 366 IF( ln_traldf_lap ) THEN ! laplacian operator (mask only) 367 367 ahtu(:,:,1:jpkm1) = ahtu(:,:,1:jpkm1) * umask(:,:,1:jpkm1) … … 381 381 !!---------------------------------------------------------------------- 382 382 !! *** ROUTINE ldf_tra *** 383 !! 383 !! 384 384 !! ** Purpose : update at kt the tracer lateral mixing coeff. (aht and aeiv) 385 385 !! … … 395 395 !! * time varying EIV coefficients: call to ldf_eiv routine 396 396 !! 397 !! ** action : ahtu, ahtv update at each time step 398 !! aeiu, aeiv - - - - (if ln_ldfeiv=T) 397 !! ** action : ahtu, ahtv update at each time step 398 !! aeiu, aeiv - - - - (if ln_ldfeiv=T) 399 399 !!---------------------------------------------------------------------- 400 400 INTEGER, INTENT(in) :: kt ! time step … … 420 420 ahtu(:,:,1) = aeiu(:,:,1) 421 421 ahtv(:,:,1) = aeiv(:,:,1) 422 ELSE ! compute aht. 422 ELSE ! compute aht. 423 423 CALL ldf_eiv( kt, aht0, ahtu, ahtv, Kmm ) 424 424 ENDIF 425 425 ! 426 z1_f20 = 1._wp / ( 2._wp * omega * SIN( rad * 20._wp ) ) ! 1 / ff(20 degrees) 426 z1_f20 = 1._wp / ( 2._wp * omega * SIN( rad * 20._wp ) ) ! 1 / ff(20 degrees) 427 427 zaht_min = 0.2_wp * aht0 ! minimum value for aht 428 zDaht = aht0 - zaht_min 429 DO_2D( 1, 1, 1, 1 ) 428 zDaht = aht0 - zaht_min 429 ! NOTE: [tiling-comms-merge] Change needed to preserve results with respect to the trunk 430 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 430 431 !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg) 431 432 !! ==>>> The Coriolis value is identical for t- & u_points, and for v- and f-points … … 479 480 !! nn_aei_ijk_t = 0 => = constant 480 481 !! ! 481 !! = 10 => = F(z) : constant with a reduction of 1/4 with depth 482 !! = 10 => = F(z) : constant with a reduction of 1/4 with depth 482 483 !! ! 483 484 !! =-20 => = F(i,j) = shape read in 'eddy_diffusivity.nc' file … … 546 547 ! != Specification of space-time variations of eaiu, aeiv 547 548 ! 548 aeiu(:,:,jpk) = 0._wp ! last level always 0 549 aeiu(:,:,jpk) = 0._wp ! last level always 0 549 550 aeiv(:,:,jpk) = 0._wp 550 551 ! ! value of EIV coef. (laplacian operator) … … 608 609 END SELECT 609 610 ! 610 IF( .NOT.l_ldfeiv_time ) THEN !* mask if No time variation 611 IF( .NOT.l_ldfeiv_time ) THEN !* mask if No time variation 611 612 DO jk = 1, jpkm1 612 613 aeiu(:,:,jk) = aeiu(:,:,jk) * umask(:,:,jk) … … 616 617 ! 617 618 ENDIF 618 ! 619 ! 619 620 END SUBROUTINE ldf_eiv_init 620 621 … … 648 649 IF( ln_traldf_triad ) THEN 649 650 DO_3D( 0, 0, 0, 0, 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 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 653 654 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 654 655 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm) 655 656 ! Compute elements required for the inverse time scale of baroclinic 656 ! eddies using the isopycnal slopes calculated in ldfslp.F : 657 ! eddies using the isopycnal slopes calculated in ldfslp.F : 657 658 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 658 659 ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) … … 662 663 ELSE 663 664 DO_3D( 0, 0, 0, 0, 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 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 667 668 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 668 669 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm) 669 670 ! Compute elements required for the inverse time scale of baroclinic 670 ! eddies using the isopycnal slopes calculated in ldfslp.F : 671 ! eddies using the isopycnal slopes calculated in ldfslp.F : 671 672 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 672 673 ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) … … 692 693 END_2D 693 694 CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1.0_wp ) ! lateral boundary condition 694 ! 695 ! 695 696 DO_2D( 0, 0, 0, 0 ) !== aei at u- and v-points ==! 696 697 paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj ) ) * umask(ji,jj,1) … … 703 704 paeiv(:,:,jk) = paeiv(:,:,1) * vmask(:,:,jk) 704 705 END DO 705 ! 706 ! 706 707 END SUBROUTINE ldf_eiv 707 708 … … 710 711 !!---------------------------------------------------------------------- 711 712 !! *** ROUTINE ldf_eiv_trp *** 712 !! 713 !! ** Purpose : add to the input ocean transport the contribution of 713 !! 714 !! ** Purpose : add to the input ocean transport the contribution of 714 715 !! the eddy induced velocity parametrization. 715 716 !! 716 717 !! ** Method : The eddy induced transport is computed from a flux stream- 717 718 !! function which depends on the slope of iso-neutral surfaces 718 !! (see ldf_slp). For example, in the i-k plan : 719 !! (see ldf_slp). For example, in the i-k plan : 719 720 !! psi_uw = mk(aeiu) e2u mi(wslpi) [in m3/s] 720 721 !! Utr_eiv = - dk[psi_uw] … … 725 726 !! ** Action : pu, pv increased by the eiv transport 726 727 !!---------------------------------------------------------------------- 727 INTEGER 728 INTEGER 729 INTEGER 730 CHARACTER(len=3) 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]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(A2D(nn_hls),jpk), INTENT(inout) :: pu ! in : 3 ocean transport components [m3/s] 733 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: pv ! out: 3 ocean transport components [m3/s] 734 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: pw ! increased by the eiv [m3/s] 734 735 !! 735 736 INTEGER :: ji, jj, jk ! dummy loop indices 736 737 REAL(wp) :: zuwk, zuwk1, zuwi, zuwi1 ! local scalars 737 738 REAL(wp) :: zvwk, zvwk1, zvwj, zvwj1 ! - - 738 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 739 !!---------------------------------------------------------------------- 740 ! 741 IF( kt == kit000 ) THEN 742 IF(lwp) WRITE(numout,*) 743 IF(lwp) WRITE(numout,*) 'ldf_eiv_trp : eddy induced advection on ', cdtype,' :' 744 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ add to velocity fields the eiv component' 745 ENDIF 746 747 739 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpsi_uw, zpsi_vw 740 !!---------------------------------------------------------------------- 741 ! 742 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 743 IF( kt == kit000 ) THEN 744 IF(lwp) WRITE(numout,*) 745 IF(lwp) WRITE(numout,*) 'ldf_eiv_trp : eddy induced advection on ', cdtype,' :' 746 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ add to velocity fields the eiv component' 747 ENDIF 748 ENDIF 749 750 748 751 zpsi_uw(:,:, 1 ) = 0._wp ; zpsi_vw(:,:, 1 ) = 0._wp 749 752 zpsi_uw(:,:,jpk) = 0._wp ; zpsi_vw(:,:,jpk) = 0._wp … … 781 784 !! 782 785 !!---------------------------------------------------------------------- 783 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT(inout) :: psi_uw, psi_vw ! streamfunction [m3/s]784 INTEGER 786 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: psi_uw, psi_vw ! streamfunction [m3/s] 787 INTEGER , INTENT(in ) :: Kmm ! ocean time level indices 785 788 ! 786 789 INTEGER :: ji, jj, jk ! dummy loop indices 787 790 REAL(wp) :: zztmp ! local scalar 788 REAL(wp), DIMENSION( jpi,jpj) :: zw2d ! 2D workspace789 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zw3d ! 3D workspace791 REAL(wp), DIMENSION(A2D(nn_hls)) :: zw2d ! 2D workspace 792 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zw3d ! 3D workspace 790 793 !!---------------------------------------------------------------------- 791 794 ! 792 795 !!gm I don't like this routine.... Crazy way of doing things, not optimal at all... 793 !!gm to be redesigned.... 796 !!gm to be redesigned.... 794 797 ! !== eiv stream function: output ==! 795 CALL lbc_lnk_multi( 'ldftra', psi_uw, 'U', -1.0_wp , psi_vw, 'V', -1.0_wp )796 !797 798 !!gm CALL iom_put( "psi_eiv_uw", psi_uw ) ! output 798 799 !!gm CALL iom_put( "psi_eiv_vw", psi_vw ) … … 802 803 zw3d(:,:,jpk) = 0._wp ! bottom value always 0 803 804 ! 804 DO jk = 1, jpkm1! e2u e3u u_eiv = -dk[psi_uw]805 zw3d( :,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u(:,:,jk,Kmm) )806 END DO805 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! e2u e3u u_eiv = -dk[psi_uw] 806 zw3d(ji,jj,jk) = ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm) ) 807 END_3D 807 808 CALL iom_put( "uoce_eiv", zw3d ) 808 809 ! 809 DO jk = 1, jpkm1! e1v e3v v_eiv = -dk[psi_vw]810 zw3d( :,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v(:,:,jk,Kmm) )811 END DO810 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! e1v e3v v_eiv = -dk[psi_vw] 811 zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm) ) 812 END_3D 812 813 CALL iom_put( "voce_eiv", zw3d ) 813 814 ! … … 816 817 & + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj ,jk) ) / e1e2t(ji,jj) 817 818 END_3D 818 CALL lbc_lnk( 'ldftra', zw3d, 'T', 1.0_wp ) ! lateral boundary condition819 819 CALL iom_put( "woce_eiv", zw3d ) 820 820 ! 821 821 IF( iom_use('weiv_masstr') ) THEN ! vertical mass transport & its square value 822 zw2d(:,:) = rho0 * e1e2t(:,:) 822 DO_2D( 0, 0, 0, 0 ) 823 zw2d(ji,jj) = rho0 * e1e2t(ji,jj) 824 END_2D 823 825 DO jk = 1, jpk 824 826 zw3d(:,:,jk) = zw3d(:,:,jk) * zw2d(:,:) 825 827 END DO 826 CALL iom_put( "weiv_masstr" , zw3d ) 828 CALL iom_put( "weiv_masstr" , zw3d ) 827 829 ENDIF 828 830 ! … … 830 832 zw3d(:,:,:) = 0.e0 831 833 DO jk = 1, jpkm1 832 zw3d(:,:,jk) = rho0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) 834 zw3d(:,:,jk) = rho0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) 833 835 END DO 834 836 CALL iom_put( "ueiv_masstr", zw3d ) ! mass transport in i-direction 835 837 ENDIF 836 838 ! 837 zztmp = 0.5_wp * rho0 * rcp 839 zztmp = 0.5_wp * rho0 * rcp 838 840 IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN 839 zw2d(:,:) = 0._wp 840 zw3d(:,:,:) = 0._wp 841 zw2d(:,:) = 0._wp 842 zw3d(:,:,:) = 0._wp 841 843 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 842 844 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) ) 845 & * ( ts (ji,jj,jk,jp_tem,Kmm) + ts (ji+1,jj,jk,jp_tem,Kmm) ) 844 846 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 845 847 END_3D 846 CALL lbc_lnk( 'ldftra', zw2d, 'U', -1.0_wp )847 CALL lbc_lnk( 'ldftra', zw3d, 'U', -1.0_wp )848 848 CALL iom_put( "ueiv_heattr" , zztmp * zw2d ) ! heat transport in i-direction 849 849 CALL iom_put( "ueiv_heattr3d", zztmp * zw3d ) ! heat transport in i-direction … … 853 853 zw3d(:,:,:) = 0.e0 854 854 DO jk = 1, jpkm1 855 zw3d(:,:,jk) = rho0 * ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) 855 zw3d(:,:,jk) = rho0 * ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) 856 856 END DO 857 857 CALL iom_put( "veiv_masstr", zw3d ) ! mass transport in i-direction 858 858 ENDIF 859 859 ! 860 zw2d(:,:) = 0._wp 861 zw3d(:,:,:) = 0._wp 860 zw2d(:,:) = 0._wp 861 zw3d(:,:,:) = 0._wp 862 862 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 863 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) ) 864 & * ( ts (ji,jj,jk,jp_tem,Kmm) + ts (ji,jj+1,jk,jp_tem,Kmm) ) 865 865 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 866 866 END_3D 867 CALL lbc_lnk( 'ldftra', zw2d, 'V', -1.0_wp ) 868 CALL iom_put( "veiv_heattr", zztmp * zw2d ) ! heat transport in j-direction 869 CALL iom_put( "veiv_heattr", zztmp * zw3d ) ! heat transport in j-direction 867 CALL iom_put( "veiv_heattr" , zztmp * zw2d ) ! heat transport in j-direction 868 CALL iom_put( "veiv_heattr3d", zztmp * zw3d ) ! heat transport in j-direction 870 869 ! 871 870 IF( iom_use( 'sophteiv' ) ) CALL dia_ptr_hst( jp_tem, 'eiv', 0.5 * zw3d ) … … 873 872 zztmp = 0.5_wp * 0.5 874 873 IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN 875 zw2d(:,:) = 0._wp 876 zw3d(:,:,:) = 0._wp 874 zw2d(:,:) = 0._wp 875 zw3d(:,:,:) = 0._wp 877 876 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 878 877 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) ) 878 & * ( ts (ji,jj,jk,jp_sal,Kmm) + ts (ji+1,jj,jk,jp_sal,Kmm) ) 880 879 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 881 880 END_3D 882 CALL lbc_lnk( 'ldftra', zw2d, 'U', -1.0_wp )883 CALL lbc_lnk( 'ldftra', zw3d, 'U', -1.0_wp )884 881 CALL iom_put( "ueiv_salttr", zztmp * zw2d ) ! salt transport in i-direction 885 882 CALL iom_put( "ueiv_salttr3d", zztmp * zw3d ) ! salt transport in i-direction 886 883 ENDIF 887 zw2d(:,:) = 0._wp 888 zw3d(:,:,:) = 0._wp 884 zw2d(:,:) = 0._wp 885 zw3d(:,:,:) = 0._wp 889 886 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 890 887 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) ) 888 & * ( ts (ji,jj,jk,jp_sal,Kmm) + ts (ji,jj+1,jk,jp_sal,Kmm) ) 892 889 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 893 890 END_3D 894 CALL lbc_lnk( 'ldftra', zw2d, 'V', -1.0_wp ) 895 CALL iom_put( "veiv_salttr", zztmp * zw2d ) ! salt transport in j-direction 896 CALL iom_put( "veiv_salttr", zztmp * zw3d ) ! salt transport in j-direction 891 CALL iom_put( "veiv_salttr" , zztmp * zw2d ) ! salt transport in j-direction 892 CALL iom_put( "veiv_salttr3d", zztmp * zw3d ) ! salt transport in j-direction 897 893 ! 898 894 IF( iom_use( 'sopsteiv' ) ) CALL dia_ptr_hst( jp_sal, 'eiv', 0.5 * zw3d )
Note: See TracChangeset
for help on using the changeset viewer.
