Changeset 9490 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF
- Timestamp:
- 2018-04-23T10:44:07+02:00 (6 years ago)
- Location:
- branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF
- Files:
-
- 4 edited
-
ldfc1d_c2d.F90 (modified) (10 diffs)
-
ldfdyn.F90 (modified) (9 diffs)
-
ldfslp.F90 (modified) (4 diffs)
-
ldftra.F90 (modified) (24 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldfc1d_c2d.F90
r9094 r9490 26 26 PUBLIC ldf_c2d ! called by ldftra and ldfdyn modules 27 27 28 REAL(wp) :: r1_2 = 0.5_wp ! =1/2 29 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 30 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 28 31 29 !! * Substitutions30 # include "vectopt_loop_substitute.h90"31 32 !!---------------------------------------------------------------------- 32 33 !! NEMO/OPA 3.7 , NEMO Consortium (2015) … … 36 37 CONTAINS 37 38 38 SUBROUTINE ldf_c1d( cd_type, p rat, pahs1, pahs2, pah1, pah2 )39 SUBROUTINE ldf_c1d( cd_type, pahs1, pahs2, pah1, pah2 ) 39 40 !!---------------------------------------------------------------------- 40 41 !! *** ROUTINE ldf_c1d *** … … 43 44 !! 44 45 !! ** Method : 1D eddy diffusivity coefficients F( depth ) 45 !! Reduction by pratfrom surface to bottom46 !! Reduction by zratio from surface to bottom 46 47 !! hyperbolic tangent profile with inflection point 47 48 !! at zh=500m and a width of zw=200m … … 50 51 !! DYN pah1, pah2 defined at T- and F-points 51 52 !!---------------------------------------------------------------------- 52 CHARACTER(len=2) , INTENT(in ) :: cd_type ! DYNamique or TRAcers 53 REAL(wp) , INTENT(in ) :: prat ! ratio surface/deep values [-] 53 CHARACTER(len=3) , INTENT(in ) :: cd_type ! DYNamique or TRAcers 54 54 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pahs1, pahs2 ! surface value of eddy coefficient [m2/s] 55 55 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pah1 , pah2 ! eddy coefficient [m2/s] … … 58 58 REAL(wp) :: zh, zc, zdep1 ! local scalars 59 59 REAL(wp) :: zw , zdep2 ! - - 60 REAL(wp) :: zratio ! - - 60 61 !!---------------------------------------------------------------------- 61 62 IF(lwp) THEN 63 WRITE(numout,*) 64 WRITE(numout,*) ' ldf_c1d : set a given profile to eddy diffusivity/viscosity coefficients' 65 WRITE(numout,*) ' ~~~~~~~' 66 ENDIF 67 62 ! 63 IF(lwp) WRITE(numout,*) 64 IF(lwp) WRITE(numout,*) ' ldf_c1d : set a given profile to eddy mixing coefficients' 65 ! 68 66 ! initialization of the profile 67 zratio = 0.25_wp ! surface/bottom ratio 69 68 zh = 500._wp ! depth of the inflection point [m] 70 69 zw = 1._wp / 200._wp ! width^-1 - - - [1/m] 71 70 ! ! associated coefficient [-] 72 zc = ( 1._wp - prat) / ( 1._wp + TANH( zh * zw) )71 zc = ( 1._wp - zratio ) / ( 1._wp + TANH( zh * zw) ) 73 72 ! 74 73 ! … … 76 75 ! 77 76 CASE( 'DYN' ) ! T- and F-points 78 DO jk = 1, jpk! pah1 at T-point79 pah1(:,:,jk) = pahs1(:,:) * ( prat + zc * ( 1._wp + TANH( - ( gdept_n(:,:,jk) - zh ) * zw) ) ) * tmask(:,:,jk)77 DO jk = jpkm1, 1, -1 ! pah1 at T-point 78 pah1(:,:,jk) = pahs1(:,:) * ( zratio + zc * ( 1._wp + TANH( - ( gdept_0(:,:,jk) - zh ) * zw) ) ) 80 79 END DO 81 DO jk = 1, jpk! pah2 at F-point (zdep2 is an approximation in zps-coord.)80 DO jk = jpkm1, 1, -1 ! pah2 at F-point (zdep2 is an approximation in zps-coord.) 82 81 DO jj = 1, jpjm1 83 DO ji = 1, fs_jpim184 zdep2 = ( gdept_ n(ji,jj+1,jk) + gdept_n(ji+1,jj+1,jk) &85 & + gdept_ n(ji,jj ,jk) + gdept_n(ji+1,jj ,jk) ) * 0.25_wp86 pah2(ji,jj,jk) = pahs2(ji,jj) * ( prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) * fmask(ji,jj,jk)82 DO ji = 1, jpim1 83 zdep2 = ( gdept_0(ji,jj+1,jk) + gdept_0(ji+1,jj+1,jk) & 84 & + gdept_0(ji,jj ,jk) + gdept_0(ji+1,jj ,jk) ) * r1_4 85 pah2(ji,jj,jk) = pahs2(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) 87 86 END DO 88 87 END DO … … 91 90 ! 92 91 CASE( 'TRA' ) ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.) 93 DO jk = 1, jpk92 DO jk = jpkm1, 1, -1 94 93 DO jj = 1, jpjm1 95 DO ji = 1, fs_jpim196 zdep1 = ( gdept_ n(ji,jj,jk) + gdept_n(ji+1,jj,jk) ) * 0.5_wp97 zdep2 = ( gdept_ n(ji,jj,jk) + gdept_n(ji,jj+1,jk) ) * 0.5_wp98 pah1(ji,jj,jk) = pahs1(ji,jj) * ( prat + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) ) ) * umask(ji,jj,jk)99 pah2(ji,jj,jk) = pahs2(ji,jj) * ( prat + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) * vmask(ji,jj,jk)94 DO ji = 1, jpim1 95 zdep1 = ( gdept_0(ji,jj,jk) + gdept_0(ji+1,jj,jk) ) * 0.5_wp 96 zdep2 = ( gdept_0(ji,jj,jk) + gdept_0(ji,jj+1,jk) ) * 0.5_wp 97 pah1(ji,jj,jk) = pahs1(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) ) ) 98 pah2(ji,jj,jk) = pahs2(ji,jj) * ( zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) ) ) 100 99 END DO 101 100 END DO … … 104 103 CALL lbc_lnk_multi( pah1, 'U', 1. , pah2, 'V', 1. ) 105 104 ! 105 CASE DEFAULT ! error 106 CALL ctl_stop( 'ldf_c1d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' ) 106 107 END SELECT 107 108 ! … … 109 110 110 111 111 SUBROUTINE ldf_c2d( cd_type, cd_op, pah0, pah1, pah2 )112 SUBROUTINE ldf_c2d( cd_type, pUfac, knn, pah1, pah2 ) 112 113 !!---------------------------------------------------------------------- 113 114 !! *** ROUTINE ldf_c2d *** … … 124 125 !! DYN pah1, pah2 defined at T- and F-points 125 126 !!---------------------------------------------------------------------- 126 CHARACTER(len=3) 127 CHARACTER(len=3) , INTENT(in ) :: cd_op ! operator:LAPlacian BiLaPlacian128 REAL(wp) , INTENT(in ) :: pah0 ! eddy coefficient [m2/s or m4/s]129 REAL(wp), DIMENSION( jpi,jpj,jpk), INTENT( out) :: pah1, pah2 ! eddy coefficient[m2/s or m4/s]127 CHARACTER(len=3) , INTENT(in ) :: cd_type ! DYNamique or TRAcers 128 REAL(wp) , INTENT(in ) :: pUfac ! =1/2*Uc LAPlacian BiLaPlacian 129 INTEGER , INTENT(in ) :: knn ! characteristic velocity [m/s] 130 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: pah1, pah2 ! eddy coefficients [m2/s or m4/s] 130 131 ! 131 132 INTEGER :: ji, jj, jk ! dummy loop indices 132 REAL(wp) :: za00, zd_max, zemax1, zemax2 ! local scalar133 INTEGER :: inn ! local integer 133 134 !!---------------------------------------------------------------------- 134 135 ! 135 zd_max = ra * rad ! = 1 degree at the equator in meters 136 IF(lwp) WRITE(numout,*) 137 IF(lwp) WRITE(numout,*) ' ldf_c2d : aht = Ufac * max(e1,e2) with Ufac = ', pUfac, ' m/s' 136 138 ! 137 IF(lwp) THEN138 WRITE(numout,*)139 WRITE(numout,*) ' ldf_c2d : aht = rn_aht0 * max(e1,e2)/e_equator ( laplacian) '140 WRITE(numout,*) ' ~~~~~~~ or = rn_bht0 * [max(e1,e2)/e_equator]**3 (bilaplacian)'141 WRITE(numout,*)142 ENDIF143 139 ! 144 SELECT CASE( cd_type ) !== surface values ==! ( depending on DYN/TRA)140 SELECT CASE( cd_type ) !== surface values ==! (chosen grid point function of DYN or TRA) 145 141 ! 146 CASE( 'DYN' ) ! T- and F-points 147 IF( cd_op == 'LAP' ) THEN ! laplacian operator 148 IF(lwp) WRITE(numout,*) ' momentum laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)' 149 za00 = pah0 / zd_max 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1) 153 zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1) 154 pah1(ji,jj,1) = za00 * zemax1 155 pah2(ji,jj,1) = za00 * zemax2 156 END DO 142 CASE( 'DYN' ) ! T- and F-points 143 DO jj = 1, jpj 144 DO ji = 1, jpi 145 pah1(ji,jj,1) = pUfac * MAX( e1t(ji,jj) , e2t(ji,jj) )**knn 146 pah2(ji,jj,1) = pUfac * MAX( e1f(ji,jj) , e2f(ji,jj) )**knn 157 147 END DO 158 ELSEIF( cd_op == 'BLP' ) THEN ! bilaplacian operator 159 IF(lwp) WRITE(numout,*) ' momentum bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3' 160 za00 = pah0 / ( zd_max * zd_max * zd_max ) 161 DO jj = 1, jpj 162 DO ji = 1, jpi 163 zemax1 = MAX( e1t(ji,jj), e2t(ji,jj) ) * tmask(ji,jj,1) 164 zemax2 = MAX( e1f(ji,jj), e2f(ji,jj) ) * fmask(ji,jj,1) 165 pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 166 pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 167 END DO 148 END DO 149 CASE( 'TRA' ) ! U- and V-points 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 pah1(ji,jj,1) = pUfac * MAX( e1u(ji,jj), e2u(ji,jj) )**knn 153 pah2(ji,jj,1) = pUfac * MAX( e1v(ji,jj), e2v(ji,jj) )**knn 168 154 END DO 169 ELSE ! NO diffusion/viscosity170 CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' )171 ENDIF172 ! ! deeper values (LAP and BLP cases)173 DO jk = 2, jpk174 pah1(:,:,jk) = pah1(:,:,1) * tmask(:,:,jk)175 pah2(:,:,jk) = pah2(:,:,1) * fmask(:,:,jk)176 155 END DO 177 ! 178 CASE( 'TRA' ) ! U- and V-points (approximation in zps-coord.) 179 IF( cd_op == 'LAP' ) THEN ! laplacian operator 180 IF(lwp) WRITE(numout,*) ' tracer laplacian coeffcients = rn_aht0/e_equ * max(e1,e2)' 181 za00 = pah0 / zd_max 182 DO jj = 1, jpj 183 DO ji = 1, jpi 184 zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) 185 zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) 186 pah1(ji,jj,1) = za00 * zemax1 187 pah2(ji,jj,1) = za00 * zemax2 188 END DO 189 END DO 190 ELSEIF( cd_op == 'BLP' ) THEN ! bilaplacian operator (NB: square root of the coeff) 191 IF(lwp) WRITE(numout,*) ' tracer bilaplacian coeffcients = rn_bht0/e_equ * max(e1,e2)**3' 192 za00 = pah0 / ( zd_max * zd_max * zd_max ) 193 DO jj = 1, jpj 194 DO ji = 1, jpi 195 zemax1 = MAX( e1u(ji,jj), e2u(ji,jj) ) * umask(ji,jj,1) 196 zemax2 = MAX( e1v(ji,jj), e2v(ji,jj) ) * vmask(ji,jj,1) 197 pah1(ji,jj,1) = za00 * zemax1 * zemax1 * zemax1 198 pah2(ji,jj,1) = za00 * zemax2 * zemax2 * zemax2 199 END DO 200 END DO 201 ELSE ! NO diffusion/viscosity 202 CALL ctl_stop( 'ldf_c2d: ', cd_op, ' case. Unknown lateral operator ' ) 203 ENDIF 204 ! ! deeper values (LAP and BLP cases) 205 DO jk = 2, jpk 206 pah1(:,:,jk) = pah1(:,:,1) * umask(:,:,jk) 207 pah2(:,:,jk) = pah2(:,:,1) * vmask(:,:,jk) 208 END DO 209 ! 156 CASE DEFAULT ! error 157 CALL ctl_stop( 'ldf_c2d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' ) 210 158 END SELECT 159 ! !== deeper values = surface one ==! (except jpk) 160 DO jk = 2, jpkm1 161 pah1(:,:,jk) = pah1(:,:,1) 162 pah2(:,:,jk) = pah2(:,:,1) 163 END DO 211 164 ! 212 165 END SUBROUTINE ldf_c2d -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
r9169 r9490 17 17 USE dom_oce ! ocean space and time domain 18 18 USE phycst ! physical constants 19 USE ldfslp ! lateral diffusion: slopes of mixing orientation 19 20 USE ldfc1d_c2d ! lateral diffusion: 1D and 2D cases 20 21 ! … … 31 32 PUBLIC ldf_dyn ! called by step.F90 32 33 33 ! 34 ! !!* Namelist namdyn_ldf : lateral mixing on momentum * 34 35 LOGICAL , PUBLIC :: ln_dynldf_NONE !: No operator (i.e. no explicit diffusion) 35 36 LOGICAL , PUBLIC :: ln_dynldf_lap !: laplacian operator … … 37 38 LOGICAL , PUBLIC :: ln_dynldf_lev !: iso-level direction 38 39 LOGICAL , PUBLIC :: ln_dynldf_hor !: horizontal (geopotential) direction 39 LOGICAL , PUBLIC :: ln_dynldf_iso !: iso-neutral direction 40 ! LOGICAL , PUBLIC :: ln_dynldf_iso !: iso-neutral direction (see ldfslp) 40 41 INTEGER , PUBLIC :: nn_ahm_ijk_t !: choice of time & space variations of the lateral eddy viscosity coef. 41 REAL(wp), PUBLIC :: rn_ahm_0 !: lateral laplacian eddy viscosity [m2/s] 42 REAL(wp), PUBLIC :: rn_ahm_b !: lateral laplacian background eddy viscosity [m2/s] 43 REAL(wp), PUBLIC :: rn_bhm_0 !: lateral bilaplacian eddy viscosity [m4/s] 44 !! If nn_ahm_ijk_t = 32 a time and space varying Smagorinsky viscosity 45 !! will be computed. 46 REAL(wp), PUBLIC :: rn_csmc !: Smagorinsky constant of proportionality 47 REAL(wp), PUBLIC :: rn_minfac !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 48 REAL(wp), PUBLIC :: rn_maxfac !: Multiplicative factor of theorectical maximum Smagorinsky viscosity 49 50 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 51 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy diffusivity coef. at U- and V-points [m2/s or m4/s] 42 ! ! time invariant coefficients: aht = 1/2 Ud*Ld (lap case) 43 ! ! bht = 1/12 Ud*Ld^3 (blp case) 44 REAL(wp), PUBLIC :: rn_Uv !: lateral viscous velocity [m/s] 45 REAL(wp), PUBLIC :: rn_Lv !: lateral viscous length [m] 46 ! ! Smagorinsky viscosity (nn_ahm_ijk_t = 32) 47 REAL(wp), PUBLIC :: rn_csmc !: Smagorinsky constant of proportionality 48 REAL(wp), PUBLIC :: rn_minfac !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 49 REAL(wp), PUBLIC :: rn_maxfac !: Multiplicative factor of theorectical maximum Smagorinsky viscosity 50 ! ! iso-neutral laplacian (ln_dynldf_lap=ln_dynldf_iso=T) 51 REAL(wp), PUBLIC :: rn_ahm_b !: lateral laplacian background eddy viscosity [m2/s] 52 53 ! !!* Parameter to control the type of lateral viscous operator 54 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 !: error in setting the operator 55 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 !: without operator (i.e. no lateral viscous trend) 56 ! !! laplacian ! bilaplacian ! 57 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 !: iso-level operator 58 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 !: iso-neutral or geopotential operator 59 ! 60 INTEGER , PUBLIC :: nldf_dyn !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 61 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 62 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy viscosity coef. at T- and F-points [m2/s or m4/s] 53 64 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dtensq !: horizontal tension squared (Smagorinsky only) 54 65 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dshesq !: horizontal shearing strain squared (Smagorinsky only) 55 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: esqt, esqf !: Square of the local gridscale (e1e2/(e1+e2))**2 56 67 57 REAL(wp) :: r1_ 12 = 1._wp / 12._wp ! =1/1268 REAL(wp) :: r1_2 = 0.5_wp ! =1/2 58 69 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 59 70 REAL(wp) :: r1_8 = 0.125_wp ! =1/8 71 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 60 72 REAL(wp) :: r1_288 = 1._wp / 288._wp ! =1/( 12^2 * 2 ) 61 73 … … 92 104 !! or L^4|D|/8 bilaplacian operator ) 93 105 !!---------------------------------------------------------------------- 94 INTEGER :: ji, jj, jk ! dummy loop indices 95 INTEGER :: ierr, inum, ios ! local integer 96 REAL(wp) :: zah0 ! local scalar 97 ! 106 INTEGER :: ji, jj, jk ! dummy loop indices 107 INTEGER :: ioptio, ierr, inum, ios, inn ! local integer 108 REAL(wp) :: zah0, zah_max, zUfac ! local scalar 109 CHARACTER(len=5) :: cl_Units ! units (m2/s or m4/s) 110 !! 98 111 NAMELIST/namdyn_ldf/ ln_dynldf_NONE, ln_dynldf_lap, ln_dynldf_blp, & ! type of operator 99 & ln_dynldf_lev , ln_dynldf_hor, ln_dynldf_iso, & ! acting direction of the operator100 & nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0, & ! lateral eddy coefficient101 & rn_csmc , rn_minfac, rn_maxfac! Smagorinsky settings112 & ln_dynldf_lev , ln_dynldf_hor, ln_dynldf_iso, & ! acting direction of the operator 113 & nn_ahm_ijk_t , rn_Uv , rn_Lv, rn_ahm_b, & ! lateral eddy coefficient 114 & rn_csmc , rn_minfac , rn_maxfac ! Smagorinsky settings 102 115 !!---------------------------------------------------------------------- 103 116 ! … … 129 142 WRITE(numout,*) ' coefficients :' 130 143 WRITE(numout,*) ' type of time-space variation nn_ahm_ijk_t = ', nn_ahm_ijk_t 131 WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0 = ', rn_ahm_0, ' m2/s' 132 WRITE(numout,*) ' background viscosity (iso case) rn_ahm_b = ', rn_ahm_b, ' m2/s' 133 WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_bhm_0 = ', rn_bhm_0, ' m4/s' 144 WRITE(numout,*) ' lateral viscous velocity (if cst) rn_Uv = ', rn_Uv, ' m/s' 145 WRITE(numout,*) ' lateral viscous length (if cst) rn_Lv = ', rn_Lv, ' m' 146 WRITE(numout,*) ' background viscosity (iso-lap case) rn_ahm_b = ', rn_ahm_b, ' m2/s' 147 ! 134 148 WRITE(numout,*) ' Smagorinsky settings (nn_ahm_ijk_t = 32) :' 135 149 WRITE(numout,*) ' Smagorinsky coefficient rn_csmc = ', rn_csmc 136 WRITE(numout,*) ' factor multiplier for theorectical lower limit for ' 137 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_minfac = ', rn_minfac 138 WRITE(numout,*) ' factor multiplier for theorectical lower upper for ' 139 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_maxfac = ', rn_maxfac 150 WRITE(numout,*) ' factor multiplier for eddy visc.' 151 WRITE(numout,*) ' lower limit (default 1.0) rn_minfac = ', rn_minfac 152 WRITE(numout,*) ' upper limit (default 1.0) rn_maxfac = ', rn_maxfac 140 153 ENDIF 141 154 142 ! ! Parameter control 155 ! 156 ! !== type of lateral operator used ==! (set nldf_dyn) 157 ! !=====================================! 158 ! 159 nldf_dyn = np_ERROR 160 ioptio = 0 161 IF( ln_dynldf_NONE ) THEN ; nldf_dyn = np_no_ldf ; ioptio = ioptio + 1 ; ENDIF 162 IF( ln_dynldf_lap ) THEN ; ioptio = ioptio + 1 ; ENDIF 163 IF( ln_dynldf_blp ) THEN ; ioptio = ioptio + 1 ; ENDIF 164 IF( ioptio /= 1 ) CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 165 ! 166 IF(.NOT.ln_dynldf_NONE ) THEN !== direction ==>> type of operator ==! 167 ioptio = 0 168 IF( ln_dynldf_lev ) ioptio = ioptio + 1 169 IF( ln_dynldf_hor ) ioptio = ioptio + 1 170 IF( ln_dynldf_iso ) ioptio = ioptio + 1 171 IF( ioptio /= 1 ) CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 direction options (level/hor/iso)' ) 172 ! 173 ! ! Set nldf_dyn, the type of lateral diffusion, from ln_dynldf_... logicals 174 ierr = 0 175 IF( ln_dynldf_lap ) THEN ! laplacian operator 176 IF( ln_zco ) THEN ! z-coordinate 177 IF ( ln_dynldf_lev ) nldf_dyn = np_lap ! iso-level = horizontal (no rotation) 178 IF ( ln_dynldf_hor ) nldf_dyn = np_lap ! iso-level = horizontal (no rotation) 179 IF ( ln_dynldf_iso ) nldf_dyn = np_lap_i ! iso-neutral ( rotation) 180 ENDIF 181 IF( ln_zps ) THEN ! z-coordinate with partial step 182 IF ( ln_dynldf_lev ) nldf_dyn = np_lap ! iso-level (no rotation) 183 IF ( ln_dynldf_hor ) nldf_dyn = np_lap ! iso-level (no rotation) 184 IF ( ln_dynldf_iso ) nldf_dyn = np_lap_i ! iso-neutral ( rotation) 185 ENDIF 186 IF( ln_sco ) THEN ! s-coordinate 187 IF ( ln_dynldf_lev ) nldf_dyn = np_lap ! iso-level = horizontal (no rotation) 188 IF ( ln_dynldf_hor ) nldf_dyn = np_lap_i ! horizontal ( rotation) 189 IF ( ln_dynldf_iso ) nldf_dyn = np_lap_i ! iso-neutral ( rotation) 190 ENDIF 191 ENDIF 192 ! 193 IF( ln_dynldf_blp ) THEN ! bilaplacian operator 194 IF( ln_zco ) THEN ! z-coordinate 195 IF( ln_dynldf_lev ) nldf_dyn = np_blp ! iso-level = horizontal (no rotation) 196 IF( ln_dynldf_hor ) nldf_dyn = np_blp ! iso-level = horizontal (no rotation) 197 IF( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 198 ENDIF 199 IF( ln_zps ) THEN ! z-coordinate with partial step 200 IF( ln_dynldf_lev ) nldf_dyn = np_blp ! iso-level (no rotation) 201 IF( ln_dynldf_hor ) nldf_dyn = np_blp ! iso-level (no rotation) 202 IF( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 203 ENDIF 204 IF( ln_sco ) THEN ! s-coordinate 205 IF( ln_dynldf_lev ) nldf_dyn = np_blp ! iso-level (no rotation) 206 IF( ln_dynldf_hor ) ierr = 2 ! horizontal ( rotation) 207 IF( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 208 ENDIF 209 ENDIF 210 ! 211 IF( ierr == 2 ) CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 212 ! 213 IF( nldf_dyn == np_lap_i ) l_ldfslp = .TRUE. ! rotation require the computation of the slopes 214 ! 215 ENDIF 216 ! 217 IF(lwp) THEN 218 WRITE(numout,*) 219 SELECT CASE( nldf_dyn ) 220 CASE( np_no_ldf ) ; WRITE(numout,*) ' ==>>> NO lateral viscosity' 221 CASE( np_lap ) ; WRITE(numout,*) ' ==>>> iso-level laplacian operator' 222 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> rotated laplacian operator with iso-level background' 223 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> iso-level bi-laplacian operator' 224 END SELECT 225 WRITE(numout,*) 226 ENDIF 227 228 ! 229 ! !== Space/time variation of eddy coefficients ==! 230 ! !=================================================! 231 ! 232 l_ldfdyn_time = .FALSE. ! no time variation except in case defined below 233 ! 143 234 IF( ln_dynldf_NONE ) THEN 144 235 IF(lwp) WRITE(numout,*) ' ==>>> No viscous operator selected. ahmt and ahmf are not allocated' 145 l_ldfdyn_time = .FALSE.146 236 RETURN 147 ENDIF 148 ! 149 IF( ln_dynldf_blp .AND. ln_dynldf_iso ) THEN ! iso-neutral bilaplacian not implemented 150 CALL ctl_stop( 'dyn_ldf_init: iso-neutral bilaplacian not coded yet' ) 151 ENDIF 152 153 ! ... Space/Time variation of eddy coefficients 154 ! ! allocate the ahm arrays 155 ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 156 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 157 ! 158 ahmt(:,:,jpk) = 0._wp ! last level always 0 159 ahmf(:,:,jpk) = 0._wp 160 ! 161 ! ! value of eddy mixing coef. 162 IF ( ln_dynldf_lap ) THEN ; zah0 = rn_ahm_0 ! laplacian operator 163 ELSEIF( ln_dynldf_blp ) THEN ; zah0 = ABS( rn_bhm_0 ) ! bilaplacian operator 164 ELSE ! NO viscous operator 165 CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' ) 166 ENDIF 167 ! 168 l_ldfdyn_time = .FALSE. ! no time variation except in case defined below 169 ! 170 IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! only if a lateral diffusion operator is used 171 ! 172 SELECT CASE( nn_ahm_ijk_t ) ! Specification of space time variations of ahmt, ahmf 237 ! 238 ELSE !== a lateral diffusion operator is used ==! 239 ! 240 ! ! allocate the ahm arrays 241 ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 242 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 243 ! 244 ahmt(:,:,jpk) = 0._wp ! last level always 0 245 ahmf(:,:,jpk) = 0._wp 246 ! 247 ! ! value of lap/blp eddy mixing coef. 248 IF( ln_dynldf_lap ) THEN ; zUfac = r1_2 *rn_Uv ; inn = 1 ; cl_Units = ' m2/s' ! laplacian 249 ELSEIF( ln_dynldf_blp ) THEN ; zUfac = r1_12*rn_Uv ; inn = 3 ; cl_Units = ' m4/s' ! bilaplacian 250 ENDIF 251 zah0 = zUfac * rn_Lv**inn ! mixing coefficient 252 zah_max = zUfac * (ra*rad)**inn ! maximum reachable coefficient (value at the Equator) 253 ! 254 SELECT CASE( nn_ahm_ijk_t ) !* Specification of space-time variations of ahmt, ahmf 173 255 ! 174 256 CASE( 0 ) !== constant ==! 175 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = constant '176 ahmt(:,:, :) = zah0 * tmask(:,:,:)177 ahmf(:,:, :) = zah0 * fmask(:,:,:)257 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity. = constant = ', zah0, cl_Units 258 ahmt(:,:,1:jpkm1) = zah0 259 ahmf(:,:,1:jpkm1) = zah0 178 260 ! 179 261 CASE( 10 ) !== fixed profile ==! 180 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = F( depth )' 181 ahmt(:,:,1) = zah0 * tmask(:,:,1) ! constant surface value 182 ahmf(:,:,1) = zah0 * fmask(:,:,1) 183 CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 262 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( depth )' 263 IF(lwp) WRITE(numout,*) ' surface viscous coef. = constant = ', zah0, cl_Units 264 ahmt(:,:,1) = zah0 ! constant surface value 265 ahmf(:,:,1) = zah0 266 CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 184 267 ! 185 268 CASE ( -20 ) !== fixed horizontal shape read in file ==! 186 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef.= F(i,j) read in eddy_viscosity.nc file'269 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F(i,j) read in eddy_viscosity.nc file' 187 270 CALL iom_open( 'eddy_viscosity_2D.nc', inum ) 188 271 CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) ) 189 272 CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) ) 190 273 CALL iom_close( inum ) 191 !!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ???192 !! do we introduce a scaling by the max value of the array, and then multiply by zah0 ????193 !! better: check that the max is <=1 i.e. it is a shape from 0 to 1, not a coef that has physical dimension194 274 DO jk = 2, jpkm1 195 ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk)196 ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk)275 ahmt(:,:,jk) = ahmt(:,:,1) 276 ahmf(:,:,jk) = ahmf(:,:,1) 197 277 END DO 198 278 ! 199 279 CASE( 20 ) !== fixed horizontal shape ==! 200 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)' 201 IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor 202 IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor^3 280 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)' 281 IF(lwp) WRITE(numout,*) ' using a fixed viscous velocity = ', rn_Uv ,' m/s and Lv = Max(e1,e2)' 282 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)' 283 CALL ldf_c2d( 'DYN', zUfac , inn , ahmt, ahmf ) ! surface value proportional to scale factor^inn 203 284 ! 204 285 CASE( -30 ) !== fixed 3D shape read in file ==! 205 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'286 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F(i,j,k) read in eddy_viscosity_3D.nc file' 206 287 CALL iom_open( 'eddy_viscosity_3D.nc', inum ) 207 288 CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt ) 208 289 CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf ) 209 290 CALL iom_close( inum ) 210 !!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ????211 !! do we introduce a scaling by the max value of the array, and then multiply by zah0 ????212 DO jk = 1, jpkm1213 ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk)214 ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk)215 END DO216 291 ! 217 292 CASE( 30 ) !== fixed 3D shape ==! 218 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef.= F( latitude, longitude, depth )'219 IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor220 IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor221 ! ! reduction with depth222 CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )293 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( latitude, longitude, depth )' 294 IF(lwp) WRITE(numout,*) ' using a fixed viscous velocity = ', rn_Uv ,' m/s and Ld = Max(e1,e2)' 295 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)' 296 CALL ldf_c2d( 'DYN', zUfac , inn , ahmt, ahmf ) ! surface value proportional to scale factor^inn 297 CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) ! reduction with depth 223 298 ! 224 299 CASE( 31 ) !== time varying 3D field ==! 225 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef.= F( latitude, longitude, depth , time )'226 IF(lwp) WRITE(numout,*) ' proportional to the velocity : |u|e/12 or |u|e^3/12'300 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( latitude, longitude, depth , time )' 301 IF(lwp) WRITE(numout,*) ' proportional to the local velocity : 1/2 |u|e (lap) or 1/12 |u|e^3 (blp)' 227 302 ! 228 303 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 229 304 ! 230 305 CASE( 32 ) !== time varying 3D field ==! 231 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = F( latitude, longitude, depth , time )' 232 IF(lwp) WRITE(numout,*) ' proportional to the local deformation rate and gridscale (Smagorinsky)' 233 IF(lwp) WRITE(numout,*) ' : L^2|D| or L^4|D|/8' 306 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( latitude, longitude, depth , time )' 307 IF(lwp) WRITE(numout,*) ' proportional to the local deformation rate and gridscale (Smagorinsky)' 234 308 ! 235 309 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 236 310 ! 237 ! allocate arrays used in ldf_dyn.311 ! ! allocate arrays used in ldf_dyn. 238 312 ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr ) 239 313 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays') 240 314 ! 241 ! Set local gridscale values 242 DO jj = 2, jpjm1 315 DO jj = 2, jpjm1 ! Set local gridscale values 243 316 DO ji = fs_2, fs_jpim1 244 317 esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2 … … 251 324 END SELECT 252 325 ! 253 IF( ln_dynldf_blp .AND. .NOT. l_ldfdyn_time ) THEN ! bilapcian and no time variation: 254 ahmt(:,:,:) = SQRT( ahmt(:,:,:) ) ! take the square root of the coefficient 255 ahmf(:,:,:) = SQRT( ahmf(:,:,:) ) 326 IF( .NOT.l_ldfdyn_time ) THEN !* No time variation 327 IF( ln_dynldf_lap ) THEN ! laplacian operator (mask only) 328 ahmt(:,:,1:jpkm1) = ahmt(:,:,1:jpkm1) * tmask(:,:,1:jpkm1) 329 ahmf(:,:,1:jpkm1) = ahmf(:,:,1:jpkm1) * fmask(:,:,1:jpkm1) 330 ELSEIF( ln_dynldf_blp ) THEN ! bilaplacian operator (square root + mask) 331 ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) 332 ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1) 333 ENDIF 256 334 ENDIF 257 335 ! … … 341 419 & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) & 342 420 & * r1_e2t(ji,jj) * e1t(ji,jj) ) * tmask(ji,jj,jk) 343 dtensq(ji,jj) = zdb *zdb421 dtensq(ji,jj) = zdb * zdb 344 422 END DO 345 423 END DO … … 351 429 & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) & 352 430 & * r1_e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,jk) 353 dshesq(ji,jj) = zdb *zdb431 dshesq(ji,jj) = zdb * zdb 354 432 END DO 355 433 END DO … … 385 463 ! 386 464 IF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8) 387 ! = sqrt( A_lap_smag L^2/8 ) 388 ! stability limits already applied to laplacian values 389 ! effective default limits are |U|L^3/12 < B_hm < L^4/(32*2dt) 390 ! 465 ! ! = sqrt( A_lap_smag L^2/8 ) 466 ! ! stability limits already applied to laplacian values 467 ! ! effective default limits are 1/12 |U|L^3 < B_hm < 1//(32*2dt) L^4 391 468 DO jk = 1, jpkm1 392 !393 469 DO jj = 2, jpjm1 394 470 DO ji = fs_2, fs_jpim1 395 ! 396 ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 397 ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 398 ! 471 ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 472 ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 399 473 END DO 400 474 END DO -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r9124 r9490 22 22 USE oce ! ocean dynamics and tracers 23 23 USE dom_oce ! ocean space and time domain 24 USE ldfdyn ! lateral diffusion: eddy viscosity coef.24 ! USE ldfdyn ! lateral diffusion: eddy viscosity coef. 25 25 USE phycst ! physical constants 26 26 USE zdfmxl ! mixed layer depth … … 43 43 LOGICAL , PUBLIC :: l_ldfslp = .FALSE. !: slopes flag 44 44 45 LOGICAL , PUBLIC :: ln_traldf_iso = .TRUE. !: iso-neutral direction 46 LOGICAL , PUBLIC :: ln_traldf_triad = .FALSE. !: griffies triad scheme 47 48 LOGICAL , PUBLIC :: ln_triad_iso = .FALSE. !: pure horizontal mixing in ML 49 LOGICAL , PUBLIC :: ln_botmix_triad = .FALSE. !: mixing on bottom 50 REAL(wp), PUBLIC :: rn_sw_triad = 1._wp !: =1 switching triads ; =0 all four triads used 51 REAL(wp), PUBLIC :: rn_slpmax = 0.01_wp !: slope limit 45 LOGICAL , PUBLIC :: ln_traldf_iso = .TRUE. !: iso-neutral direction (nam_traldf namelist) 46 LOGICAL , PUBLIC :: ln_traldf_triad = .FALSE. !: griffies triad scheme (nam_traldf namelist) 47 LOGICAL , PUBLIC :: ln_dynldf_iso !: iso-neutral direction (nam_dynldf namelist) 48 49 LOGICAL , PUBLIC :: ln_triad_iso = .FALSE. !: pure horizontal mixing in ML (nam_traldf namelist) 50 LOGICAL , PUBLIC :: ln_botmix_triad = .FALSE. !: mixing on bottom (nam_traldf namelist) 51 REAL(wp), PUBLIC :: rn_sw_triad = 1._wp !: =1 switching triads ; =0 all four triads used (nam_traldf namelist) 52 REAL(wp), PUBLIC :: rn_slpmax = 0.01_wp !: slope limit (nam_traldf namelist) 52 53 53 54 LOGICAL , PUBLIC :: l_grad_zps = .FALSE. !: special treatment for Horz Tgradients w partial steps (triad operator) … … 749 750 ! 750 751 IF( ln_traldf_triad ) THEN ! Griffies operator : triad of slopes 751 IF(lwp) WRITE(numout,*) ' Griffies (triad) operator initialisation'752 IF(lwp) WRITE(numout,*) ' ==>>> triad) operator (Griffies)' 752 753 ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , & 753 754 & triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , & … … 757 758 ! 758 759 ELSE ! Madec operator : slopes at u-, v-, and w-points 759 IF(lwp) WRITE(numout,*) ' Madec operator initialisation'760 IF(lwp) WRITE(numout,*) ' ==>>> iso operator (Madec)' 760 761 ALLOCATE( omlmask(jpi,jpj,jpk) , & 761 762 & uslp(jpi,jpj,jpk) , uslpml(jpi,jpj) , wslpi(jpi,jpj,jpk) , wslpiml(jpi,jpj) , & -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
r9169 r9490 50 50 LOGICAL , PUBLIC :: ln_traldf_hor !: horizontal (geopotential) direction 51 51 ! LOGICAL , PUBLIC :: ln_traldf_iso !: iso-neutral direction (see ldfslp) 52 ! != iso-neutral options =! 52 53 ! LOGICAL , PUBLIC :: ln_traldf_triad !: griffies triad scheme (see ldfslp) 53 54 LOGICAL , PUBLIC :: ln_traldf_msc !: Method of Stabilizing Correction … … 58 59 ! != Coefficients =! 59 60 INTEGER , PUBLIC :: nn_aht_ijk_t !: choice of time & space variations of the lateral eddy diffusivity coef. 60 REAL(wp), PUBLIC :: rn_aht_0 !: laplacian lateral eddy diffusivity [m2/s] 61 REAL(wp), PUBLIC :: rn_bht_0 !: bilaplacian lateral eddy diffusivity [m4/s] 62 63 ! !!* Namelist namtra_ldfeiv : eddy induced velocity param. * 61 ! ! time invariant coefficients: aht_0 = 1/2 Ud*Ld (lap case) 62 ! ! bht_0 = 1/12 Ud*Ld^3 (blp case) 63 REAL(wp), PUBLIC :: rn_Ud !: lateral diffusive velocity [m/s] 64 REAL(wp), PUBLIC :: rn_Ld !: lateral diffusive length [m] 65 66 ! !!* Namelist namtra_eiv : eddy induced velocity param. * 64 67 ! != Use/diagnose eiv =! 65 68 LOGICAL , PUBLIC :: ln_ldfeiv !: eddy induced velocity flag … … 67 70 ! != Coefficients =! 68 71 INTEGER , PUBLIC :: nn_aei_ijk_t !: choice of time/space variation of the eiv coeff. 69 REAL(wp), PUBLIC :: rn_aeiv_0 !: eddy induced velocity coefficient [m2/s] 72 REAL(wp), PUBLIC :: rn_Ue !: lateral diffusive velocity [m/s] 73 REAL(wp), PUBLIC :: rn_Le !: lateral diffusive length [m] 70 74 75 ! ! Flag to control the type of lateral diffusive operator 76 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in specification of lateral diffusion 77 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 ! without operator (i.e. no lateral diffusive trend) 78 ! !! laplacian ! bilaplacian ! 79 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 ! iso-level operator 80 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 , np_blp_i = 21 ! standard iso-neutral or geopotential operator 81 INTEGER, PARAMETER, PUBLIC :: np_lap_it = 12 , np_blp_it = 22 ! triad iso-neutral or geopotential operator 82 83 INTEGER , PUBLIC :: nldf_tra = 0 !: type of lateral diffusion used defined from ln_traldf_... (namlist logicals) 71 84 LOGICAL , PUBLIC :: l_ldftra_time = .FALSE. !: flag for time variation of the lateral eddy diffusivity coef. 72 LOGICAL , PUBLIC :: l_ldfeiv_time = .FALSE. ! flag for time variation of the eiv coef.85 LOGICAL , PUBLIC :: l_ldfeiv_time = .FALSE. !: flag for time variation of the eiv coef. 73 86 74 87 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahtu, ahtv !: eddy diffusivity coef. at U- and V-points [m2/s] 75 88 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aeiu, aeiv !: eddy induced velocity coeff. [m2/s] 76 89 90 REAL(wp) :: aht0, aei0 ! constant eddy coefficients (deduced from namelist values) [m2/s] 91 REAL(wp) :: r1_2 = 0.5_wp ! =1/2 77 92 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 78 93 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 … … 108 123 !! =-30 => = F(i,j,k) = shape read in 'eddy_diffusivity.nc' file 109 124 !! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10) 110 !! = 31 = F(i,j,k,t) = F(local velocity) ( |u|e /12laplacian operator111 !! or |u|e^3/12bilaplacian operator )125 !! = 31 = F(i,j,k,t) = F(local velocity) ( 1/2 |u|e laplacian operator 126 !! or 1/12 |u|e^3 bilaplacian operator ) 112 127 !! * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init 113 128 !! 114 !! ** action : ahtu, ahtv initialized once for all or l_ldftra_time set to true 115 !! aeiu, aeiv initialized once for all or l_ldfeiv_time set to true 116 !!---------------------------------------------------------------------- 117 INTEGER :: jk ! dummy loop indices 118 INTEGER :: ierr, inum, ios ! local integer 119 REAL(wp) :: zah0 ! local scalar 129 !! ** action : ahtu, ahtv initialized one for all or l_ldftra_time set to true 130 !! aeiu, aeiv initialized one for all or l_ldfeiv_time set to true 131 !!---------------------------------------------------------------------- 132 INTEGER :: jk ! dummy loop indices 133 INTEGER :: ioptio, ierr, inum, ios, inn ! local integer 134 REAL(wp) :: zah_max, zUfac ! - - 135 CHARACTER(len=5) :: cl_Units ! units (m2/s or m4/s) 120 136 !! 121 137 NAMELIST/namtra_ldf/ ln_traldf_NONE, ln_traldf_lap , ln_traldf_blp , & ! type of operator … … 123 139 & ln_traldf_iso , ln_traldf_msc , rn_slpmax , & ! option for iso-neutral operator 124 140 & ln_triad_iso , ln_botmix_triad, rn_sw_triad , & ! option for triad operator 125 & rn_aht_0 , rn_bht_0 , nn_aht_ijk_t! lateral eddy coefficient141 & nn_aht_ijk_t , rn_Ud , rn_Ld ! lateral eddy coefficient 126 142 !!---------------------------------------------------------------------- 127 143 ! 128 144 IF(lwp) THEN ! control print 129 145 WRITE(numout,*) 130 WRITE(numout,*) 'ldf_tra_init : lateral tracer physics'146 WRITE(numout,*) 'ldf_tra_init : lateral tracer diffusion' 131 147 WRITE(numout,*) '~~~~~~~~~~~~ ' 132 148 ENDIF 149 133 150 ! 134 151 ! Choice of lateral tracer physics … … 154 171 WRITE(numout,*) ' iso-neutral Madec operator ln_traldf_iso = ', ln_traldf_iso 155 172 WRITE(numout,*) ' iso-neutral triad operator ln_traldf_triad = ', ln_traldf_triad 156 WRITE(numout,*) ' iso-neutral (Method of Stab. Corr.)ln_traldf_msc = ', ln_traldf_msc173 WRITE(numout,*) ' use the Method of Stab. Correction ln_traldf_msc = ', ln_traldf_msc 157 174 WRITE(numout,*) ' maximum isoppycnal slope rn_slpmax = ', rn_slpmax 158 175 WRITE(numout,*) ' pure lateral mixing in ML ln_triad_iso = ', ln_triad_iso … … 160 177 WRITE(numout,*) ' lateral mixing on bottom ln_botmix_triad = ', ln_botmix_triad 161 178 WRITE(numout,*) ' coefficients :' 162 WRITE(numout,*) ' lateral eddy diffusivity (lap case) rn_aht_0 = ', rn_aht_0163 WRITE(numout,*) ' lateral eddy diffusivity (bilap case) rn_bht_0 = ', rn_bht_0164 179 WRITE(numout,*) ' type of time-space variation nn_aht_ijk_t = ', nn_aht_ijk_t 165 ENDIF 166 ! 167 ! ! Parameter control 168 ! 169 IF( ln_traldf_NONE ) THEN 170 IF(lwp) WRITE(numout,*) ' ==>>> No diffusive operator selected. ahtu and ahtv are not allocated' 171 l_ldftra_time = .FALSE. 172 RETURN 173 ENDIF 180 WRITE(numout,*) ' lateral diffusive velocity (if cst) rn_Ud = ', rn_Ud, ' m/s' 181 WRITE(numout,*) ' lateral diffusive length (if cst) rn_Ld = ', rn_Ld, ' m' 182 ENDIF 183 ! 184 ! 185 ! Operator and its acting direction (set nldf_tra) 186 ! ================================= 187 ! 188 nldf_tra = np_ERROR 189 ioptio = 0 190 IF( ln_traldf_NONE ) THEN ; nldf_tra = np_no_ldf ; ioptio = ioptio + 1 ; ENDIF 191 IF( ln_traldf_lap ) THEN ; ioptio = ioptio + 1 ; ENDIF 192 IF( ln_traldf_blp ) THEN ; ioptio = ioptio + 1 ; ENDIF 193 IF( ioptio /= 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 194 ! 195 IF( .NOT.ln_traldf_NONE ) THEN !== direction ==>> type of operator ==! 196 ioptio = 0 197 IF( ln_traldf_lev ) ioptio = ioptio + 1 198 IF( ln_traldf_hor ) ioptio = ioptio + 1 199 IF( ln_traldf_iso ) ioptio = ioptio + 1 200 IF( ioptio /= 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE direction (level/hor/iso)' ) 201 ! 202 ! ! defined the type of lateral diffusion from ln_traldf_... logicals 203 ierr = 0 204 IF ( ln_traldf_lap ) THEN ! laplacian operator 205 IF ( ln_zco ) THEN ! z-coordinate 206 IF ( ln_traldf_lev ) nldf_tra = np_lap ! iso-level = horizontal (no rotation) 207 IF ( ln_traldf_hor ) nldf_tra = np_lap ! iso-level = horizontal (no rotation) 208 IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard ( rotation) 209 IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad ( rotation) 210 ENDIF 211 IF ( ln_zps ) THEN ! z-coordinate with partial step 212 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 213 IF ( ln_traldf_hor ) nldf_tra = np_lap ! horizontal (no rotation) 214 IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard (rotation) 215 IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad (rotation) 216 ENDIF 217 IF ( ln_sco ) THEN ! s-coordinate 218 IF ( ln_traldf_lev ) nldf_tra = np_lap ! iso-level (no rotation) 219 IF ( ln_traldf_hor ) nldf_tra = np_lap_i ! horizontal ( rotation) 220 IF ( ln_traldf_iso ) nldf_tra = np_lap_i ! iso-neutral: standard ( rotation) 221 IF ( ln_traldf_triad ) nldf_tra = np_lap_it ! iso-neutral: triad ( rotation) 222 ENDIF 223 ENDIF 224 ! 225 IF( ln_traldf_blp ) THEN ! bilaplacian operator 226 IF ( ln_zco ) THEN ! z-coordinate 227 IF ( ln_traldf_lev ) nldf_tra = np_blp ! iso-level = horizontal (no rotation) 228 IF ( ln_traldf_hor ) nldf_tra = np_blp ! iso-level = horizontal (no rotation) 229 IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation) 230 IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation) 231 ENDIF 232 IF ( ln_zps ) THEN ! z-coordinate with partial step 233 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 234 IF ( ln_traldf_hor ) nldf_tra = np_blp ! horizontal (no rotation) 235 IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation) 236 IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation) 237 ENDIF 238 IF ( ln_sco ) THEN ! s-coordinate 239 IF ( ln_traldf_lev ) nldf_tra = np_blp ! iso-level (no rotation) 240 IF ( ln_traldf_hor ) nldf_tra = np_blp_it ! horizontal ( rotation) 241 IF ( ln_traldf_iso ) nldf_tra = np_blp_i ! iso-neutral: standard ( rotation) 242 IF ( ln_traldf_triad ) nldf_tra = np_blp_it ! iso-neutral: triad ( rotation) 243 ENDIF 244 ENDIF 245 IF ( ierr == 1 ) CALL ctl_stop( 'iso-level in z-partial step, not allowed' ) 246 ENDIF 247 ! 248 IF( ln_ldfeiv .AND. .NOT.( ln_traldf_iso .OR. ln_traldf_triad ) ) & 249 & CALL ctl_stop( 'ln_ldfeiv=T requires iso-neutral laplacian diffusion' ) 250 IF( ln_isfcav .AND. ln_traldf_triad ) & 251 & CALL ctl_stop( ' ice shelf cavity and traldf_triad not tested' ) 252 ! 253 IF( nldf_tra == np_lap_i .OR. nldf_tra == np_lap_it .OR. & 254 & nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it ) l_ldfslp = .TRUE. ! slope of neutral surfaces required 174 255 ! 175 256 IF( ln_traldf_blp .AND. ( ln_traldf_iso .OR. ln_traldf_triad) ) THEN ! iso-neutral bilaplacian need MSC … … 177 258 ENDIF 178 259 ! 260 IF(lwp) THEN 261 WRITE(numout,*) 262 SELECT CASE( nldf_tra ) 263 CASE( np_no_ldf ) ; WRITE(numout,*) ' ==>>> NO lateral diffusion' 264 CASE( np_lap ) ; WRITE(numout,*) ' ==>>> laplacian iso-level operator' 265 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> Rotated laplacian operator (standard)' 266 CASE( np_lap_it ) ; WRITE(numout,*) ' ==>>> Rotated laplacian operator (triad)' 267 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> bilaplacian iso-level operator' 268 CASE( np_blp_i ) ; WRITE(numout,*) ' ==>>> Rotated bilaplacian operator (standard)' 269 CASE( np_blp_it ) ; WRITE(numout,*) ' ==>>> Rotated bilaplacian operator (triad)' 270 END SELECT 271 WRITE(numout,*) 272 ENDIF 273 274 ! 179 275 ! Space/time variation of eddy coefficients 180 276 ! =========================================== 181 ! ! allocate the aht arrays 182 ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr ) 183 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays') 184 ! 185 ahtu(:,:,jpk) = 0._wp ! last level always 0 186 ahtv(:,:,jpk) = 0._wp 187 ! 188 ! ! value of eddy mixing coef. 189 IF ( ln_traldf_lap ) THEN ; zah0 = rn_aht_0 ! laplacian operator 190 ELSEIF( ln_traldf_blp ) THEN ; zah0 = ABS( rn_bht_0 ) ! bilaplacian operator 191 ENDIF 192 ! 193 l_ldftra_time = .FALSE. ! no time variation except in case defined below 194 ! 195 IF( ln_traldf_lap .OR. ln_traldf_blp ) THEN ! only if a lateral diffusion operator is used 196 ! 197 SELECT CASE( nn_aht_ijk_t ) ! Specification of space time variations of ehtu, ahtv 277 ! 278 l_ldftra_time = .FALSE. ! no time variation except in case defined below 279 ! 280 IF( ln_traldf_NONE ) THEN !== no explicit diffusive operator ==! 281 ! 282 IF(lwp) WRITE(numout,*) ' ==>>> No diffusive operator selected. ahtu and ahtv are not allocated' 283 RETURN 284 ! 285 ELSE !== a lateral diffusion operator is used ==! 286 ! 287 ! ! allocate the aht arrays 288 ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr ) 289 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays') 290 ! 291 ahtu(:,:,jpk) = 0._wp ! last level always 0 292 ahtv(:,:,jpk) = 0._wp 293 !. 294 ! ! value of lap/blp eddy mixing coef. 295 IF( ln_traldf_lap ) THEN ; zUfac = r1_2 *rn_Ud ; inn = 1 ; cl_Units = ' m2/s' ! laplacian 296 ELSEIF( ln_traldf_blp ) THEN ; zUfac = r1_12*rn_Ud ; inn = 3 ; cl_Units = ' m4/s' ! bilaplacian 297 ENDIF 298 aht0 = zUfac * rn_Ld**inn ! mixing coefficient 299 zah_max = zUfac * (ra*rad)**inn ! maximum reachable coefficient (value at the Equator for e1=1 degree) 300 ! 301 ! 302 SELECT CASE( nn_aht_ijk_t ) !* Specification of space-time variations of ahtu, ahtv 198 303 ! 199 304 CASE( 0 ) !== constant ==! 200 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = constant = ', rn_aht_0201 ahtu(:,:, :) = zah0 * umask(:,:,:)202 ahtv(:,:, :) = zah0 * vmask(:,:,:)305 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = constant = ', aht0, cl_Units 306 ahtu(:,:,1:jpkm1) = aht0 307 ahtv(:,:,1:jpkm1) = aht0 203 308 ! 204 309 CASE( 10 ) !== fixed profile ==! 205 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( depth )' 206 ahtu(:,:,1) = zah0 * umask(:,:,1) ! constant surface value 207 ahtv(:,:,1) = zah0 * vmask(:,:,1) 208 CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) 209 ! 210 CASE ( -20 ) !== fixed horizontal shape read in file ==! 211 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F(i,j) read in eddy_diffusivity.nc file' 310 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( depth )' 311 IF(lwp) WRITE(numout,*) ' surface eddy diffusivity = constant = ', aht0, cl_Units 312 ahtu(:,:,1) = aht0 ! constant surface value 313 ahtv(:,:,1) = aht0 314 CALL ldf_c1d( 'TRA', ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) 315 ! 316 CASE ( -20 ) !== fixed horizontal shape and magnitude read in file ==! 317 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F(i,j) read in eddy_diffusivity.nc file' 212 318 CALL iom_open( 'eddy_diffusivity_2D.nc', inum ) 213 319 CALL iom_get ( inum, jpdom_data, 'ahtu_2D', ahtu(:,:,1) ) … … 215 321 CALL iom_close( inum ) 216 322 DO jk = 2, jpkm1 217 ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk)218 ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk)323 ahtu(:,:,jk) = ahtu(:,:,1) 324 ahtv(:,:,jk) = ahtv(:,:,1) 219 325 END DO 220 326 ! 221 327 CASE( 20 ) !== fixed horizontal shape ==! 222 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)' 223 IF( ln_traldf_lap ) CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor 224 IF( ln_traldf_blp ) CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor 328 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)' 329 IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ud,' m/s and Ld = Max(e1,e2)' 330 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)' 331 CALL ldf_c2d( 'TRA', zUfac , inn , ahtu, ahtv ) ! value proportional to scale factor^inn 225 332 ! 226 333 CASE( 21 ) !== time varying 2D field ==! 227 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( latitude, longitude, time )' 228 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 229 IF(lwp) WRITE(numout,*) ' min value = 0.1 * rn_aht_0' 230 IF(lwp) WRITE(numout,*) ' max value = rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21)' 231 IF(lwp) WRITE(numout,*) ' increased to rn_aht_0 within 20N-20S' 334 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, time )' 335 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 336 IF(lwp) WRITE(numout,*) ' min value = 0.2 * aht0 (with aht0= 1/2 rn_Ud*rn_Ld)' 337 IF(lwp) WRITE(numout,*) ' max value = aei0 (with aei0=1/2 rn_Ue*Le increased to aht0 within 20N-20S' 232 338 ! 233 339 l_ldftra_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 234 340 ! 235 IF( ln_traldf_blp ) THEN 236 CALL ctl_stop( 'ldf_tra_init: aht=F(growth rate of baroc. insta.) incompatible with bilaplacian operator' ) 237 ENDIF 341 IF( ln_traldf_blp ) CALL ctl_stop( 'ldf_tra_init: aht=F( growth rate of baroc. insta .)', & 342 & ' incompatible with bilaplacian operator' ) 238 343 ! 239 344 CASE( -30 ) !== fixed 3D shape read in file ==! 240 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef.= F(i,j,k) read in eddy_diffusivity.nc file'345 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F(i,j,k) read in eddy_diffusivity.nc file' 241 346 CALL iom_open( 'eddy_diffusivity_3D.nc', inum ) 242 347 CALL iom_get ( inum, jpdom_data, 'ahtu_3D', ahtu ) 243 348 CALL iom_get ( inum, jpdom_data, 'ahtv_3D', ahtv ) 244 349 CALL iom_close( inum ) 245 DO jk = 1, jpkm1246 ahtu(:,:,jk) = ahtu(:,:,jk) * umask(:,:,jk)247 ahtv(:,:,jk) = ahtv(:,:,jk) * vmask(:,:,jk)248 END DO249 350 ! 250 351 CASE( 30 ) !== fixed 3D shape ==! 251 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef.= F( latitude, longitude, depth )'252 IF( ln_traldf_lap ) CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor253 IF( ln_traldf_blp ) CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor254 ! ! reduction with depth255 CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv )352 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, depth )' 353 IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ud,' m/s and Ld = Max(e1,e2)' 354 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)' 355 CALL ldf_c2d( 'TRA', zUfac , inn , ahtu, ahtv ) ! surface value proportional to scale factor^inn 356 CALL ldf_c1d( 'TRA', ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) ! reduction with depth 256 357 ! 257 358 CASE( 31 ) !== time varying 3D field ==! 258 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef.= F( latitude, longitude, depth , time )'259 IF(lwp) WRITE(numout,*) ' proportional to the velocity : |u|e/12 or |u|e^3/12'359 IF(lwp) WRITE(numout,*) ' ==>>> eddy diffusivity = F( latitude, longitude, depth , time )' 360 IF(lwp) WRITE(numout,*) ' proportional to the velocity : 1/2 |u|e or 1/12 |u|e^3' 260 361 ! 261 362 l_ldftra_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 … … 265 366 END SELECT 266 367 ! 267 IF( ln_traldf_blp .AND. .NOT. l_ldftra_time ) THEN 268 ahtu(:,:,:) = SQRT( ahtu(:,:,:) ) 269 ahtv(:,:,:) = SQRT( ahtv(:,:,:) ) 368 IF( .NOT.l_ldftra_time ) THEN !* No time variation 369 IF( ln_traldf_lap ) THEN ! laplacian operator (mask only) 370 ahtu(:,:,1:jpkm1) = ahtu(:,:,1:jpkm1) * umask(:,:,1:jpkm1) 371 ahtv(:,:,1:jpkm1) = ahtv(:,:,1:jpkm1) * vmask(:,:,1:jpkm1) 372 ELSEIF( ln_traldf_blp ) THEN ! bilaplacian operator (square root + mask) 373 ahtu(:,:,1:jpkm1) = SQRT( ahtu(:,:,1:jpkm1) ) * umask(:,:,1:jpkm1) 374 ahtv(:,:,1:jpkm1) = SQRT( ahtv(:,:,1:jpkm1) ) * vmask(:,:,1:jpkm1) 375 ENDIF 270 376 ENDIF 271 377 ! … … 281 387 !! ** Purpose : update at kt the tracer lateral mixing coeff. (aht and aeiv) 282 388 !! 283 !! ** Method : 389 !! ** Method : * time varying eddy diffusivity coefficients: 284 390 !! 285 391 !! nn_aei_ijk_t = 21 aeiu, aeiv = F(i,j, t) = F(growth rate of baroclinic instability) … … 290 396 !! or |u|e^3/12 bilaplacian operator ) 291 397 !! 398 !! * time varying EIV coefficients: call to ldf_eiv routine 399 !! 292 400 !! ** action : ahtu, ahtv update at each time step 293 401 !! aeiu, aeiv - - - - (if ln_ldfeiv=T) … … 296 404 ! 297 405 INTEGER :: ji, jj, jk ! dummy loop indices 298 REAL(wp) :: zaht, zahf, zaht_min, z 1_f20! local scalar406 REAL(wp) :: zaht, zahf, zaht_min, zDaht, z1_f20 ! local scalar 299 407 !!---------------------------------------------------------------------- 300 408 ! 301 409 IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN ! eddy induced velocity coefficients 302 410 ! ! =F(growth rate of baroclinic instability) 303 ! ! max value rn_aeiv_0 ; decreased to 0 within 20N-20S304 CALL ldf_eiv( kt, rn_aeiv_0, aeiu, aeiv )411 ! ! max value aeiv_0 ; decreased to 0 within 20N-20S 412 CALL ldf_eiv( kt, aei0, aeiu, aeiv ) 305 413 ENDIF 306 414 ! … … 308 416 ! 309 417 CASE( 21 ) !== time varying 2D field ==! = F( growth rate of baroclinic instability ) 310 ! ! min value rn_aht_0 / 10311 ! ! max value rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21)312 ! ! increase to rn_aht_0 within 20N-20S418 ! ! min value 0.2*aht0 419 ! ! max value aht0 (aei0 if nn_aei_ijk_t=21) 420 ! ! increase to aht0 within 20N-20S 313 421 IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN ! use the already computed aei. 314 422 ahtu(:,:,1) = aeiu(:,:,1) 315 423 ahtv(:,:,1) = aeiv(:,:,1) 316 424 ELSE ! compute aht. 317 CALL ldf_eiv( kt, rn_aht_0, ahtu, ahtv )425 CALL ldf_eiv( kt, aht0, ahtu, ahtv ) 318 426 ENDIF 319 427 ! 320 z1_f20 = 1._wp / ( 2._wp * omega * SIN( rad * 20._wp ) ) ! 1 / ff(20 degrees) 321 zaht_min = 0.2_wp * rn_aht_0 ! minimum value for aht 428 z1_f20 = 1._wp / ( 2._wp * omega * SIN( rad * 20._wp ) ) ! 1 / ff(20 degrees) 429 zaht_min = 0.2_wp * aht0 ! minimum value for aht 430 zDaht = aht0 - zaht_min 322 431 DO jj = 1, jpj 323 432 DO ji = 1, jpi 324 433 !!gm CAUTION : here we assume lat/lon grid in 20deg N/S band (like all ORCA cfg) 325 434 !! ==>>> The Coriolis value is identical for t- & u_points, and for v- and f-points 326 zaht = ( 1._wp - MIN( 1._wp , ABS( ff_t(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min )327 zahf = ( 1._wp - MIN( 1._wp , ABS( ff_f(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min )328 ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) * umask(ji,jj,1)! min value zaht_min329 ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zahf ) * vmask(ji,jj,1)! increase within 20S-20N330 END DO 331 END DO 332 DO jk = 2, jpkm1 ! deeper value = surface value435 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 DO 440 END DO 441 DO jk = 1, jpkm1 ! deeper value = surface value + mask for all levels 333 442 ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk) 334 443 ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk) … … 338 447 IF( ln_traldf_lap ) THEN ! laplacian operator |u| e /12 339 448 DO jk = 1, jpkm1 340 ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 449 ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 ! n.b. ub,vb are masked 341 450 ahtv(:,:,jk) = ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12 342 451 END DO … … 355 464 CALL iom_put( "ahtv_3d", ahtv(:,:,:) ) ! 3D v-eddy diffusivity coeff. 356 465 ! 357 !!gm : THE IF below is to be checked (comes from Seb)358 466 IF( ln_ldfeiv ) THEN 359 467 CALL iom_put( "aeiu_2d", aeiu(:,:,1) ) ! surface u-EIV coeff. … … 372 480 !! ** Purpose : initialization of the eiv coeff. from namelist choices. 373 481 !! 374 !! ** Method : 375 !! 376 !! ** Action : aeiu , aeiv : EIV coeff. at u- & v-points 377 !! l_ldfeiv_time : =T if EIV coefficients vary with time 378 !!---------------------------------------------------------------------- 379 INTEGER :: jk ! dummy loop indices 380 INTEGER :: ierr, inum, ios ! local integer 381 ! 382 NAMELIST/namtra_ldfeiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv) 383 & nn_aei_ijk_t, rn_aeiv_0 ! eiv coefficient 482 !! ** Method : the eiv diffusivity coef. specification depends on: 483 !! nn_aei_ijk_t = 0 => = constant 484 !! ! 485 !! = 10 => = F(z) : constant with a reduction of 1/4 with depth 486 !! ! 487 !! =-20 => = F(i,j) = shape read in 'eddy_diffusivity.nc' file 488 !! = 20 = F(i,j) = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case) 489 !! = 21 = F(i,j,t) = F(growth rate of baroclinic instability) 490 !! ! 491 !! =-30 => = F(i,j,k) = shape read in 'eddy_diffusivity.nc' file 492 !! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10) 493 !! 494 !! ** Action : aeiu , aeiv : initialized one for all or l_ldftra_time set to true 495 !! l_ldfeiv_time : =T if EIV coefficients vary with time 496 !!---------------------------------------------------------------------- 497 INTEGER :: jk ! dummy loop indices 498 INTEGER :: ierr, inum, ios, inn ! local integer 499 REAL(wp) :: zah_max, zUfac ! - scalar 500 !! 501 NAMELIST/namtra_eiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv) 502 & nn_aei_ijk_t, rn_Ue, rn_Le ! eiv coefficient 384 503 !!---------------------------------------------------------------------- 385 504 ! … … 390 509 ENDIF 391 510 ! 392 REWIND( numnam_ref ) ! Namelist namtra_ ldfeiv in reference namelist : eddy induced velocity param.393 READ ( numnam_ref, namtra_ ldfeiv, IOSTAT = ios, ERR = 901)394 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ ldfeiv in reference namelist', lwp )395 ! 396 REWIND( numnam_cfg ) ! Namelist namtra_ ldfeiv in configuration namelist : eddy induced velocity param.397 READ ( numnam_cfg, namtra_ ldfeiv, IOSTAT = ios, ERR = 902 )398 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_ ldfeiv in configuration namelist', lwp )399 IF(lwm) WRITE ( numond, namtra_ ldfeiv )511 REWIND( numnam_ref ) ! Namelist namtra_eiv in reference namelist : eddy induced velocity param. 512 READ ( numnam_ref, namtra_eiv, IOSTAT = ios, ERR = 901) 513 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_eiv in reference namelist', lwp ) 514 ! 515 REWIND( numnam_cfg ) ! Namelist namtra_eiv in configuration namelist : eddy induced velocity param. 516 READ ( numnam_cfg, namtra_eiv, IOSTAT = ios, ERR = 902 ) 517 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namtra_eiv in configuration namelist', lwp ) 518 IF(lwm) WRITE ( numond, namtra_eiv ) 400 519 401 520 IF(lwp) THEN ! control print 402 WRITE(numout,*) ' Namelist namtra_ldfeiv : ' 403 WRITE(numout,*) ' Eddy Induced Velocity (eiv) param. ln_ldfeiv = ', ln_ldfeiv 404 WRITE(numout,*) ' eiv streamfunction & velocity diag. ln_ldfeiv_dia = ', ln_ldfeiv_dia 405 WRITE(numout,*) ' eddy induced velocity coef. rn_aeiv_0 = ', rn_aeiv_0 406 WRITE(numout,*) ' type of time-space variation nn_aei_ijk_t = ', nn_aei_ijk_t 521 WRITE(numout,*) ' Namelist namtra_eiv : ' 522 WRITE(numout,*) ' Eddy Induced Velocity (eiv) param. ln_ldfeiv = ', ln_ldfeiv 523 WRITE(numout,*) ' eiv streamfunction & velocity diag. ln_ldfeiv_dia = ', ln_ldfeiv_dia 524 WRITE(numout,*) ' coefficients :' 525 WRITE(numout,*) ' type of time-space variation nn_aei_ijk_t = ', nn_aht_ijk_t 526 WRITE(numout,*) ' lateral diffusive velocity (if cst) rn_Ue = ', rn_Ue, ' m/s' 527 WRITE(numout,*) ' lateral diffusive length (if cst) rn_Le = ', rn_Le, ' m' 407 528 WRITE(numout,*) 408 529 ENDIF 409 530 ! 410 IF( ln_ldfeiv .AND. ln_traldf_blp ) CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' ) 411 412 ! ! Parameter control 413 l_ldfeiv_time = .FALSE. 414 ! 415 IF( ln_ldfeiv ) THEN ! allocate the aei arrays 531 l_ldfeiv_time = .FALSE. ! no time variation except in case defined below 532 ! 533 ! 534 IF( .NOT.ln_ldfeiv ) THEN !== Parametrization not used ==! 535 ! 536 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity param is NOT used' 537 ln_ldfeiv_dia = .FALSE. 538 ! 539 ELSE !== use the parametrization ==! 540 ! 541 IF(lwp) WRITE(numout,*) ' ==>>> use eddy induced velocity parametrization' 542 IF(lwp) WRITE(numout,*) 543 ! 544 IF( ln_traldf_blp ) CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' ) 545 ! 546 ! != allocate the aei arrays 416 547 ALLOCATE( aeiu(jpi,jpj,jpk), aeiv(jpi,jpj,jpk), STAT=ierr ) 417 548 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ldf_eiv: failed to allocate arrays') 418 549 ! 419 SELECT CASE( nn_aei_ijk_t ) ! Specification of space time variations of eaiu, aeiv 420 ! 421 CASE( 0 ) !== constant ==! 422 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = constant = ', rn_aeiv_0 423 aeiu(:,:,:) = rn_aeiv_0 424 aeiv(:,:,:) = rn_aeiv_0 425 ! 426 CASE( 10 ) !== fixed profile ==! 550 ! != Specification of space-time variations of eaiu, aeiv 551 ! 552 aeiu(:,:,jpk) = 0._wp ! last level always 0 553 aeiv(:,:,jpk) = 0._wp 554 ! ! value of EIV coef. (laplacian operator) 555 zUfac = r1_2 *rn_Ue ! velocity factor 556 inn = 1 ! L-exponent 557 aei0 = zUfac * rn_Le**inn ! mixing coefficient 558 zah_max = zUfac * (ra*rad)**inn ! maximum reachable coefficient (value at the Equator) 559 560 SELECT CASE( nn_aei_ijk_t ) !* Specification of space-time variations 561 ! 562 CASE( 0 ) !-- constant --! 563 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = constant = ', aei0, ' m2/s' 564 aeiu(:,:,1:jpkm1) = aei0 565 aeiv(:,:,1:jpkm1) = aei0 566 ! 567 CASE( 10 ) !-- fixed profile --! 427 568 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( depth )' 428 aeiu(:,:,1) = rn_aeiv_0 ! constant surface value 429 aeiv(:,:,1) = rn_aeiv_0 430 CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 431 ! 432 CASE ( -20 ) !== fixed horizontal shape read in file ==! 433 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F(i,j) read in eddy_diffusivity_2D.nc file' 569 IF(lwp) WRITE(numout,*) ' surface eddy diffusivity = constant = ', aht0, ' m2/s' 570 aeiu(:,:,1) = aei0 ! constant surface value 571 aeiv(:,:,1) = aei0 572 CALL ldf_c1d( 'TRA', aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 573 ! 574 CASE ( -20 ) !-- fixed horizontal shape read in file --! 575 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F(i,j) read in eddy_diffusivity_2D.nc file' 434 576 CALL iom_open ( 'eddy_induced_velocity_2D.nc', inum ) 435 577 CALL iom_get ( inum, jpdom_data, 'aeiu', aeiu(:,:,1) ) 436 578 CALL iom_get ( inum, jpdom_data, 'aeiv', aeiv(:,:,1) ) 437 579 CALL iom_close( inum ) 438 DO jk = 2, jpk 580 DO jk = 2, jpkm1 439 581 aeiu(:,:,jk) = aeiu(:,:,1) 440 582 aeiv(:,:,jk) = aeiv(:,:,1) 441 583 END DO 442 584 ! 443 CASE( 20 ) !== fixed horizontal shape ==! 444 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or bilap case)' 445 CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv ) ! surface value proportional to scale factor 446 ! 447 CASE( 21 ) !== time varying 2D field ==! 448 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( latitude, longitude, time )' 449 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 585 CASE( 20 ) !-- fixed horizontal shape --! 586 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( e1, e2 )' 587 IF(lwp) WRITE(numout,*) ' using a fixed diffusive velocity = ', rn_Ue, ' m/s and Le = Max(e1,e2)' 588 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, ' m2/s for e1=1°)' 589 CALL ldf_c2d( 'TRA', zUfac , inn , aeiu, aeiv ) ! value proportional to scale factor^inn 590 ! 591 CASE( 21 ) !-- time varying 2D field --! 592 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( latitude, longitude, time )' 593 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 594 IF(lwp) WRITE(numout,*) ' maximum allowed value: aei0 = ', aei0, ' m2/s' 450 595 ! 451 596 l_ldfeiv_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 452 597 ! 453 CASE( -30 ) !== fixed 3D shape read in file ==!454 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixingcoef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'598 CASE( -30 ) !-- fixed 3D shape read in file --! 599 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file' 455 600 CALL iom_open ( 'eddy_induced_velocity_3D.nc', inum ) 456 601 CALL iom_get ( inum, jpdom_data, 'aeiu', aeiu ) … … 458 603 CALL iom_close( inum ) 459 604 ! 460 CASE( 30 ) !== fixed 3D shape ==! 461 IF(lwp) WRITE(numout,*) ' ==>>> tracer mixing coef. = F( latitude, longitude, depth )' 462 CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv ) ! surface value proportional to scale factor 463 ! ! reduction with depth 464 CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 605 CASE( 30 ) !-- fixed 3D shape --! 606 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity coef. = F( latitude, longitude, depth )' 607 CALL ldf_c2d( 'TRA', zUfac , inn , aeiu, aeiv ) ! surface value proportional to scale factor^inn 608 CALL ldf_c1d( 'TRA', aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) ! reduction with depth 465 609 ! 466 610 CASE DEFAULT … … 468 612 END SELECT 469 613 ! 470 ELSE 471 IF(lwp) WRITE(numout,*) ' ==>>> eddy induced velocity param is NOT used neither diagnosed' 472 ln_ldfeiv_dia = .FALSE. 614 IF( .NOT.l_ldfeiv_time ) THEN !* mask if No time variation 615 DO jk = 1, jpkm1 616 aeiu(:,:,jk) = aeiu(:,:,jk) * umask(:,:,jk) 617 ahtv(:,:,jk) = ahtv(:,:,jk) * vmask(:,:,jk) 618 END DO 619 ENDIF 620 ! 473 621 ENDIF 474 622 ! … … 493 641 INTEGER :: ji, jj, jk ! dummy loop indices 494 642 REAL(wp) :: zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei ! local scalars 495 REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, z ross, zaeiw ! 2D workspace496 !!---------------------------------------------------------------------- 497 ! 498 zn (:,:) = 0._wp! Local initialization499 zhw 500 zah 501 z ross(:,:) = 0._wp643 REAL(wp), DIMENSION(jpi,jpj) :: zn, zah, zhw, zRo, zaeiw ! 2D workspace 644 !!---------------------------------------------------------------------- 645 ! 646 zn (:,:) = 0._wp ! Local initialization 647 zhw(:,:) = 5._wp 648 zah(:,:) = 0._wp 649 zRo(:,:) = 0._wp 502 650 ! ! Compute lateral diffusive coefficient at T-point 503 651 IF( ln_traldf_triad ) THEN … … 538 686 END DO 539 687 END DO 540 END 688 ENDIF 541 689 542 690 DO jj = 2, jpjm1 543 691 DO ji = fs_2, fs_jpim1 ! vector opt. 544 692 zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 545 ! Rossby radius at w-point taken < 40km and > 2km546 z ross(ji,jj) = MAX( MIN( .4 * zn(ji,jj) / zfw, 40.e3 ), 2.e3)693 ! Rossby radius at w-point taken betwenn 2 km and 40km 694 zRo(ji,jj) = MAX( 2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 ) ) 547 695 ! Compute aeiw by multiplying Ro^2 and T^-1 548 zaeiw(ji,jj) = z ross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)696 zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 549 697 END DO 550 698 END DO … … 554 702 DO jj = 2, jpjm1 555 703 DO ji = fs_2, fs_jpim1 ! vector opt. 556 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) 704 zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease 557 705 zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0 558 706 END DO … … 620 768 DO jj = 1, jpjm1 621 769 DO ji = 1, fs_jpim1 ! vector opt. 622 zpsi_uw(ji,jj,jk) = - 0.25_wp* e2u(ji,jj) * ( wslpi(ji,jj,jk ) + wslpi(ji+1,jj,jk) ) &623 & 624 zpsi_vw(ji,jj,jk) = - 0.25_wp* e1v(ji,jj) * ( wslpj(ji,jj,jk ) + wslpj(ji,jj+1,jk) ) &625 & 770 zpsi_uw(ji,jj,jk) = - r1_4 * e2u(ji,jj) * ( wslpi(ji,jj,jk ) + wslpi(ji+1,jj,jk) ) & 771 & * ( aeiu (ji,jj,jk-1) + aeiu (ji ,jj,jk) ) * umask(ji,jj,jk) 772 zpsi_vw(ji,jj,jk) = - r1_4 * e1v(ji,jj) * ( wslpj(ji,jj,jk ) + wslpj(ji,jj+1,jk) ) & 773 & * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj ,jk) ) * vmask(ji,jj,jk) 626 774 END DO 627 775 END DO
Note: See TracChangeset
for help on using the changeset viewer.
