Changeset 14072 for NEMO/trunk/src/OCE/LDF/ldftra.F90
- Timestamp:
- 2020-12-04T08:48:38+01:00 (3 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/OCE/LDF/ldftra.F90 (modified) (37 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/LDF/ldftra.F90
r13982 r14072 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 428 zDaht = aht0 - zaht_min 429 429 ! NOTE: [tiling-comms-merge] Change needed to preserve results with respect to the trunk 430 430 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) … … 480 480 !! nn_aei_ijk_t = 0 => = constant 481 481 !! ! 482 !! = 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 483 483 !! ! 484 484 !! =-20 => = F(i,j) = shape read in 'eddy_diffusivity.nc' file … … 547 547 ! != Specification of space-time variations of eaiu, aeiv 548 548 ! 549 aeiu(:,:,jpk) = 0._wp ! last level always 0 549 aeiu(:,:,jpk) = 0._wp ! last level always 0 550 550 aeiv(:,:,jpk) = 0._wp 551 551 ! ! value of EIV coef. (laplacian operator) … … 609 609 END SELECT 610 610 ! 611 IF( .NOT.l_ldfeiv_time ) THEN !* mask if No time variation 611 IF( .NOT.l_ldfeiv_time ) THEN !* mask if No time variation 612 612 DO jk = 1, jpkm1 613 613 aeiu(:,:,jk) = aeiu(:,:,jk) * umask(:,:,jk) … … 617 617 ! 618 618 ENDIF 619 ! 619 ! 620 620 END SUBROUTINE ldf_eiv_init 621 621 … … 649 649 IF( ln_traldf_triad ) THEN 650 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 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 654 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 655 655 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm) 656 656 ! Compute elements required for the inverse time scale of baroclinic 657 ! eddies using the isopycnal slopes calculated in ldfslp.F : 657 ! eddies using the isopycnal slopes calculated in ldfslp.F : 658 658 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 659 659 ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) … … 663 663 ELSE 664 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 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 668 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 669 669 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w(ji,jj,jk,Kmm) 670 670 ! Compute elements required for the inverse time scale of baroclinic 671 ! eddies using the isopycnal slopes calculated in ldfslp.F : 671 ! eddies using the isopycnal slopes calculated in ldfslp.F : 672 672 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 673 673 ze3w = e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) … … 693 693 END_2D 694 694 CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1.0_wp ) ! lateral boundary condition 695 ! 695 ! 696 696 DO_2D( 0, 0, 0, 0 ) !== aei at u- and v-points ==! 697 697 paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj ) ) * umask(ji,jj,1) … … 704 704 paeiv(:,:,jk) = paeiv(:,:,1) * vmask(:,:,jk) 705 705 END DO 706 ! 706 ! 707 707 END SUBROUTINE ldf_eiv 708 708 … … 711 711 !!---------------------------------------------------------------------- 712 712 !! *** ROUTINE ldf_eiv_trp *** 713 !! 714 !! ** Purpose : add to the input ocean transport the contribution of 713 !! 714 !! ** Purpose : add to the input ocean transport the contribution of 715 715 !! the eddy induced velocity parametrization. 716 716 !! 717 717 !! ** Method : The eddy induced transport is computed from a flux stream- 718 718 !! function which depends on the slope of iso-neutral surfaces 719 !! (see ldf_slp). For example, in the i-k plan : 719 !! (see ldf_slp). For example, in the i-k plan : 720 720 !! psi_uw = mk(aeiu) e2u mi(wslpi) [in m3/s] 721 721 !! Utr_eiv = - dk[psi_uw] … … 748 748 ENDIF 749 749 750 750 751 751 zpsi_uw(:,:, 1 ) = 0._wp ; zpsi_vw(:,:, 1 ) = 0._wp 752 752 zpsi_uw(:,:,jpk) = 0._wp ; zpsi_vw(:,:,jpk) = 0._wp … … 794 794 ! 795 795 !!gm I don't like this routine.... Crazy way of doing things, not optimal at all... 796 !!gm to be redesigned.... 796 !!gm to be redesigned.... 797 797 ! !== eiv stream function: output ==! 798 798 !!gm CALL iom_put( "psi_eiv_uw", psi_uw ) ! output … … 826 826 zw3d(:,:,jk) = zw3d(:,:,jk) * zw2d(:,:) 827 827 END DO 828 CALL iom_put( "weiv_masstr" , zw3d ) 828 CALL iom_put( "weiv_masstr" , zw3d ) 829 829 ENDIF 830 830 ! … … 832 832 zw3d(:,:,:) = 0.e0 833 833 DO jk = 1, jpkm1 834 zw3d(:,:,jk) = rho0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) 834 zw3d(:,:,jk) = rho0 * ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) 835 835 END DO 836 836 CALL iom_put( "ueiv_masstr", zw3d ) ! mass transport in i-direction 837 837 ENDIF 838 838 ! 839 zztmp = 0.5_wp * rho0 * rcp 839 zztmp = 0.5_wp * rho0 * rcp 840 840 IF( iom_use('ueiv_heattr') .OR. iom_use('ueiv_heattr3d') ) THEN 841 zw2d(:,:) = 0._wp 842 zw3d(:,:,:) = 0._wp 841 zw2d(:,:) = 0._wp 842 zw3d(:,:,:) = 0._wp 843 843 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 844 844 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_uw(ji,jj,jk+1) - psi_uw(ji ,jj,jk) ) & 845 & * ( 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) ) 846 846 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 847 847 END_3D … … 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 … … 872 872 zztmp = 0.5_wp * 0.5 873 873 IF( iom_use('ueiv_salttr') .OR. iom_use('ueiv_salttr3d')) THEN 874 zw2d(:,:) = 0._wp 875 zw3d(:,:,:) = 0._wp 874 zw2d(:,:) = 0._wp 875 zw3d(:,:,:) = 0._wp 876 876 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 877 877 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) * ( psi_uw(ji,jj,jk+1) - psi_uw(ji ,jj,jk) ) & 878 & * ( 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) ) 879 879 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 880 880 END_3D … … 882 882 CALL iom_put( "ueiv_salttr3d", zztmp * zw3d ) ! salt transport in i-direction 883 883 ENDIF 884 zw2d(:,:) = 0._wp 885 zw3d(:,:,:) = 0._wp 884 zw2d(:,:) = 0._wp 885 zw3d(:,:,:) = 0._wp 886 886 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 887 887 zw3d(ji,jj,jk) = zw3d(ji,jj,jk) + ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj ,jk) ) & 888 & * ( 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) ) 889 889 zw2d(ji,jj) = zw2d(ji,jj) + zw3d(ji,jj,jk) 890 890 END_3D
Note: See TracChangeset
for help on using the changeset viewer.
