- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- Location:
- branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/LDF
- Files:
-
- 14 deleted
- 3 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
r4624 r6225 6 6 !! History : OPA ! 1997-07 (G. Madec) multi dimensional coefficients 7 7 !! NEMO 1.0 ! 2002-09 (G. Madec) F90: Free form and module 8 !! ----------------------------------------------------------------------9 10 !!---------------------------------------------------------------------- 11 !! ldf_dyn_init : initialization, namelist read, and parameters control 12 !! ldf_dyn_c3d : 3D eddy viscosity coefficient initialization13 !! ldf_dyn_ c2d : 2D eddy viscosity coefficient initialization14 !! ldf_dyn _c1d : 1D eddy viscosity coefficient initialization8 !! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification, 9 !! ! add velocity dependent coefficient and optional read in file 10 !!---------------------------------------------------------------------- 11 12 !!---------------------------------------------------------------------- 13 !! ldf_dyn_init : initialization, namelist read, and parameters control 14 !! ldf_dyn : update lateral eddy viscosity coefficients at each time step 15 15 !!---------------------------------------------------------------------- 16 16 USE oce ! ocean dynamics and tracers 17 17 USE dom_oce ! ocean space and time domain 18 USE ldfdyn_oce ! ocean dynamics lateral physics19 18 USE phycst ! physical constants 20 USE ldf slp ! ???21 USE ioipsl19 USE ldfc1d_c2d ! lateral diffusion: 1D and 2D cases 20 ! 22 21 USE in_out_manager ! I/O manager 22 USE iom ! I/O module for ehanced bottom friction file 23 USE timing ! Timing 23 24 USE lib_mpp ! distribued memory computing library 24 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 28 29 PRIVATE 29 30 30 PUBLIC ldf_dyn_init ! called by opa.F90 31 32 INTERFACE ldf_zpf 33 MODULE PROCEDURE ldf_zpf_1d, ldf_zpf_1d_3d, ldf_zpf_3d 34 END INTERFACE 31 PUBLIC ldf_dyn_init ! called by nemogcm.F90 32 PUBLIC ldf_dyn ! called by step.F90 33 34 ! !!* Namelist namdyn_ldf : lateral mixing on momentum * 35 LOGICAL , PUBLIC :: ln_dynldf_lap !: laplacian operator 36 LOGICAL , PUBLIC :: ln_dynldf_blp !: bilaplacian operator 37 LOGICAL , PUBLIC :: ln_dynldf_lev !: iso-level direction 38 LOGICAL , PUBLIC :: ln_dynldf_hor !: horizontal (geopotential) direction 39 LOGICAL , PUBLIC :: ln_dynldf_iso !: iso-neutral direction 40 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 45 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 46 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy diffusivity coef. at U- and V-points [m2/s or m4/s] 48 49 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 50 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 51 REAL(wp) :: r1_288 = 1._wp / 288._wp ! =1/( 12^2 * 2 ) 35 52 36 53 !! * Substitutions 37 # include " domzgr_substitute.h90"38 !!---------------------------------------------------------------------- 39 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)54 # include "vectopt_loop_substitute.h90" 55 !!---------------------------------------------------------------------- 56 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 40 57 !! $Id$ 41 58 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 49 66 !! ** Purpose : set the horizontal ocean dynamics physics 50 67 !! 51 !! ** Method : 52 !! - default option : ahm = constant coef. = rn_ahm_0 (namelist) 53 !! - 'key_dynldf_c1d': ahm = F(depth) see ldf_dyn_c1d.h90 54 !! - 'key_dynldf_c2d': ahm = F(latitude,longitude) see ldf_dyn_c2d.h90 55 !! - 'key_dynldf_c3d': ahm = F(latitude,longitude,depth) see ldf_dyn_c3d.h90 56 !! 57 !! N.B. User defined include files. By default, 3d and 2d coef. 58 !! are set to a constant value given in the namelist and the 1d 59 !! coefficients are initialized to a hyperbolic tangent vertical 60 !! profile. 61 !! 62 !! Reference : Madec, G. and M. Imbard, 1996: Climate Dynamics, 12, 381-388. 63 !!---------------------------------------------------------------------- 64 INTEGER :: ioptio ! ??? 65 INTEGER :: ios ! Local : output status for namelist read 66 LOGICAL :: ll_print = .FALSE. ! Logical flag for printing viscosity coef. 67 !! 68 NAMELIST/namdyn_ldf/ ln_dynldf_lap , ln_dynldf_bilap, & 69 & ln_dynldf_level, ln_dynldf_hor , ln_dynldf_iso, & 70 & rn_ahm_0_lap , rn_ahmb_0 , rn_ahm_0_blp , & 71 & rn_cmsmag_1 , rn_cmsmag_2 , rn_cmsh, & 72 & rn_ahm_m_lap , rn_ahm_m_blp 73 74 !!---------------------------------------------------------------------- 75 68 !! ** Method : the eddy viscosity coef. specification depends on: 69 !! - the operator: 70 !! ln_dynldf_lap = T laplacian operator 71 !! ln_dynldf_blp = T bilaplacian operator 72 !! - the parameter nn_ahm_ijk_t: 73 !! nn_ahm_ijk_t = 0 => = constant 74 !! = 10 => = F(z) : = constant with a reduction of 1/4 with depth 75 !! =-20 => = F(i,j) = shape read in 'eddy_viscosity.nc' file 76 !! = 20 = F(i,j) = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case) 77 !! =-30 => = F(i,j,k) = shape read in 'eddy_viscosity.nc' file 78 !! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10) 79 !! = 31 = F(i,j,k,t) = F(local velocity) ( |u|e /12 laplacian operator 80 !! or |u|e^3/12 bilaplacian operator ) 81 !!---------------------------------------------------------------------- 82 INTEGER :: jk ! dummy loop indices 83 INTEGER :: ierr, inum, ios ! local integer 84 REAL(wp) :: zah0 ! local scalar 85 ! 86 NAMELIST/namdyn_ldf/ ln_dynldf_lap, ln_dynldf_blp, & 87 & ln_dynldf_lev, ln_dynldf_hor, ln_dynldf_iso, & 88 & nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0 89 !!---------------------------------------------------------------------- 90 ! 76 91 REWIND( numnam_ref ) ! Namelist namdyn_ldf in reference namelist : Lateral physics 77 92 READ ( numnam_ref, namdyn_ldf, IOSTAT = ios, ERR = 901) … … 88 103 WRITE(numout,*) '~~~~~~~' 89 104 WRITE(numout,*) ' Namelist namdyn_ldf : set lateral mixing parameters' 90 WRITE(numout,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap 91 WRITE(numout,*) ' bilaplacian operator ln_dynldf_bilap = ', ln_dynldf_bilap 92 WRITE(numout,*) ' iso-level ln_dynldf_level = ', ln_dynldf_level 93 WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor = ', ln_dynldf_hor 94 WRITE(numout,*) ' iso-neutral ln_dynldf_iso = ', ln_dynldf_iso 95 WRITE(numout,*) ' horizontal laplacian eddy viscosity rn_ahm_0_lap = ', rn_ahm_0_lap 96 WRITE(numout,*) ' background viscosity rn_ahmb_0 = ', rn_ahmb_0 97 WRITE(numout,*) ' horizontal bilaplacian eddy viscosity rn_ahm_0_blp = ', rn_ahm_0_blp 98 WRITE(numout,*) ' upper limit for laplacian eddy visc rn_ahm_m_lap = ', rn_ahm_m_lap 99 WRITE(numout,*) ' upper limit for bilap eddy viscosity rn_ahm_m_blp = ', rn_ahm_m_blp 100 101 ENDIF 102 103 ahm0 = rn_ahm_0_lap ! OLD namelist variables defined from DOCTOR namelist variables 104 ahmb0 = rn_ahmb_0 105 ahm0_blp = rn_ahm_0_blp 106 107 ! ... check of lateral diffusive operator on tracers 108 ! ==> will be done in trazdf module 109 110 ! ... Space variation of eddy coefficients 111 ioptio = 0 112 #if defined key_dynldf_c3d 113 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth)' 114 ioptio = ioptio+1 115 #endif 116 #if defined key_dynldf_c2d 117 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude)' 118 ioptio = ioptio+1 119 #endif 120 #if defined key_dynldf_c1d 121 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( depth )' 122 ioptio = ioptio+1 123 IF( ln_sco ) CALL ctl_stop( 'key_dynldf_c1d cannot be used in s-coordinate (ln_sco)' ) 124 #endif 125 IF( ioptio == 0 ) THEN 126 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = constant (default option)' 127 ELSEIF( ioptio > 1 ) THEN 128 CALL ctl_stop( 'use only one of the following keys: key_dynldf_c3d, key_dynldf_c2d, key_dynldf_c1d' ) 129 ENDIF 130 131 132 IF( ln_dynldf_bilap ) THEN 133 IF(lwp) WRITE(numout,*) ' biharmonic momentum diffusion' 134 IF( .NOT. ln_dynldf_lap ) ahm0 = ahm0_blp ! Allow spatially varying coefs, which use ahm0 as input 135 IF( ahm0_blp > 0 .AND. .NOT. lk_esopa ) CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be negative' ) 136 ELSE 137 IF(lwp) WRITE(numout,*) ' harmonic momentum diff. (default)' 138 IF( ahm0 < 0 .AND. .NOT. lk_esopa ) CALL ctl_stop( 'The horizontal viscosity coef. ahm0 must be positive' ) 139 ENDIF 140 141 142 ! Lateral eddy viscosity 143 ! ====================== 144 #if defined key_dynldf_c3d 145 CALL ldf_dyn_c3d( ll_print ) ! ahm = 3D coef. = F( longitude, latitude, depth ) 146 #elif defined key_dynldf_c2d 147 CALL ldf_dyn_c2d( ll_print ) ! ahm = 1D coef. = F( longitude, latitude ) 148 #elif defined key_dynldf_c1d 149 CALL ldf_dyn_c1d( ll_print ) ! ahm = 1D coef. = F( depth ) 150 #else 151 ! Constant coefficients 152 IF(lwp) WRITE(numout,*) 153 IF(lwp) WRITE(numout,*) 'inildf: constant eddy viscosity coef. ' 154 IF(lwp) WRITE(numout,*) '~~~~~~' 155 IF(lwp) WRITE(numout,*) ' ahm1 = ahm2 = ahm0 = ',ahm0 156 #endif 157 nkahm_smag = 0 158 #if defined key_dynldf_smag 159 nkahm_smag = 1 160 #endif 161 105 ! 106 WRITE(numout,*) ' type :' 107 WRITE(numout,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap 108 WRITE(numout,*) ' bilaplacian operator ln_dynldf_blp = ', ln_dynldf_blp 109 ! 110 WRITE(numout,*) ' direction of action :' 111 WRITE(numout,*) ' iso-level ln_dynldf_lev = ', ln_dynldf_lev 112 WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor = ', ln_dynldf_hor 113 WRITE(numout,*) ' iso-neutral ln_dynldf_iso = ', ln_dynldf_iso 114 ! 115 WRITE(numout,*) ' coefficients :' 116 WRITE(numout,*) ' type of time-space variation nn_ahm_ijk_t = ', nn_ahm_ijk_t 117 WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0_lap = ', rn_ahm_0, ' m2/s' 118 WRITE(numout,*) ' background viscosity (iso case) rn_ahm_b = ', rn_ahm_b, ' m2/s' 119 WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_ahm_0_blp = ', rn_bhm_0, ' m4/s' 120 ENDIF 121 122 ! ! Parameter control 123 IF( .NOT.ln_dynldf_lap .AND. .NOT.ln_dynldf_blp ) THEN 124 IF(lwp) WRITE(numout,*) ' No viscous operator selected. ahmt and ahmf are not allocated' 125 l_ldfdyn_time = .FALSE. 126 RETURN 127 ENDIF 128 ! 129 IF( ln_dynldf_blp .AND. ln_dynldf_iso ) THEN ! iso-neutral bilaplacian not implemented 130 CALL ctl_stop( 'dyn_ldf_init: iso-neutral bilaplacian not coded yet' ) 131 ENDIF 132 133 ! ... Space/Time variation of eddy coefficients 134 ! ! allocate the ahm arrays 135 ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 136 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 137 ! 138 ahmt(:,:,jpk) = 0._wp ! last level always 0 139 ahmf(:,:,jpk) = 0._wp 140 ! 141 ! ! value of eddy mixing coef. 142 IF ( ln_dynldf_lap ) THEN ; zah0 = rn_ahm_0 ! laplacian operator 143 ELSEIF( ln_dynldf_blp ) THEN ; zah0 = ABS( rn_bhm_0 ) ! bilaplacian operator 144 ELSE ! NO viscous operator 145 CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' ) 146 ENDIF 147 ! 148 l_ldfdyn_time = .FALSE. ! no time variation except in case defined below 149 ! 150 IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! only if a lateral diffusion operator is used 151 ! 152 SELECT CASE( nn_ahm_ijk_t ) ! Specification of space time variations of ahmt, ahmf 153 ! 154 CASE( 0 ) !== constant ==! 155 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = constant ' 156 ahmt(:,:,:) = zah0 * tmask(:,:,:) 157 ahmf(:,:,:) = zah0 * fmask(:,:,:) 158 ! 159 CASE( 10 ) !== fixed profile ==! 160 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( depth )' 161 ahmt(:,:,1) = zah0 * tmask(:,:,1) ! constant surface value 162 ahmf(:,:,1) = zah0 * fmask(:,:,1) 163 CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 164 ! 165 CASE ( -20 ) !== fixed horizontal shape read in file ==! 166 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F(i,j) read in eddy_viscosity.nc file' 167 CALL iom_open( 'eddy_viscosity_2D.nc', inum ) 168 CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) ) 169 CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) ) 170 CALL iom_close( inum ) 171 !!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ??? 172 !! do we introduce a scaling by the max value of the array, and then multiply by zah0 ???? 173 !! better: check that the max is <=1 i.e. it is a shape from 0 to 1, not a coef that has physical dimension 174 DO jk = 2, jpkm1 175 ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk) 176 ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk) 177 END DO 178 ! 179 CASE( 20 ) !== fixed horizontal shape ==! 180 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)' 181 IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor 182 IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor^3 183 ! 184 CASE( -30 ) !== fixed 3D shape read in file ==! 185 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file' 186 CALL iom_open( 'eddy_viscosity_3D.nc', inum ) 187 CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt ) 188 CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf ) 189 CALL iom_close( inum ) 190 !!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ???? 191 !! do we introduce a scaling by the max value of the array, and then multiply by zah0 ???? 192 DO jk = 1, jpkm1 193 ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk) 194 ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk) 195 END DO 196 ! 197 CASE( 30 ) !== fixed 3D shape ==! 198 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth )' 199 IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor 200 IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor 201 ! ! reduction with depth 202 CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 203 ! 204 CASE( 31 ) !== time varying 3D field ==! 205 IF(lwp) WRITE(numout,*) ' momentum mixing coef. = F( latitude, longitude, depth , time )' 206 IF(lwp) WRITE(numout,*) ' proportional to the velocity : |u|e/12 or |u|e^3/12' 207 ! 208 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 209 ! 210 CASE DEFAULT 211 CALL ctl_stop('ldf_dyn_init: wrong choice for nn_ahm_ijk_t, the type of space-time variation of ahm') 212 END SELECT 213 ! 214 IF( ln_dynldf_blp .AND. .NOT. l_ldfdyn_time ) THEN ! bilapcian and no time variation: 215 ahmt(:,:,:) = SQRT( ahmt(:,:,:) ) ! take the square root of the coefficient 216 ahmf(:,:,:) = SQRT( ahmf(:,:,:) ) 217 ENDIF 218 ! 219 ENDIF 162 220 ! 163 221 END SUBROUTINE ldf_dyn_init 164 222 165 #if defined key_dynldf_c3d 166 # include "ldfdyn_c3d.h90" 167 #elif defined key_dynldf_c2d 168 # include "ldfdyn_c2d.h90" 169 #elif defined key_dynldf_c1d 170 # include "ldfdyn_c1d.h90" 171 #endif 172 173 174 SUBROUTINE ldf_zpf_1d( ld_print, pdam, pwam, pbot, pdep, pah ) 175 !!---------------------------------------------------------------------- 176 !! *** ROUTINE ldf_zpf *** 177 !! 178 !! ** Purpose : vertical adimensional profile for eddy coefficient 179 !! 180 !! ** Method : 1D eddy viscosity coefficients ( depth ) 181 !!---------------------------------------------------------------------- 182 LOGICAL , INTENT(in ) :: ld_print ! If true, output arrays on numout 183 REAL(wp), INTENT(in ) :: pdam ! depth of the inflection point 184 REAL(wp), INTENT(in ) :: pwam ! width of inflection 185 REAL(wp), INTENT(in ) :: pbot ! bottom value (0<pbot<= 1) 186 REAL(wp), INTENT(in ), DIMENSION(jpk) :: pdep ! depth of the gridpoint (T, U, V, F) 187 REAL(wp), INTENT(inout), DIMENSION(jpk) :: pah ! adimensional vertical profile 188 !! 189 INTEGER :: jk ! dummy loop indices 190 REAL(wp) :: zm00, zm01, zmhb, zmhs ! temporary scalars 191 !!---------------------------------------------------------------------- 192 193 zm00 = TANH( ( pdam - gdept_1d(1 ) ) / pwam ) 194 zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam ) 195 zmhs = zm00 / zm01 196 zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01 197 198 DO jk = 1, jpk 199 pah(jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam ) ) 200 END DO 201 202 IF(lwp .AND. ld_print ) THEN ! Control print 203 WRITE(numout,*) 204 WRITE(numout,*) ' ahm profile : ' 205 WRITE(numout,*) 206 WRITE(numout,'(" jk ahm "," depth t-level " )') 207 DO jk = 1, jpk 208 WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(jk), pdep(jk) 209 END DO 210 ENDIF 211 ! 212 END SUBROUTINE ldf_zpf_1d 213 214 215 SUBROUTINE ldf_zpf_1d_3d( ld_print, pdam, pwam, pbot, pdep, pah ) 216 !!---------------------------------------------------------------------- 217 !! *** ROUTINE ldf_zpf *** 218 !! 219 !! ** Purpose : vertical adimensional profile for eddy coefficient 220 !! 221 !! ** Method : 1D eddy viscosity coefficients ( depth ) 222 !!---------------------------------------------------------------------- 223 LOGICAL , INTENT(in ) :: ld_print ! If true, output arrays on numout 224 REAL(wp), INTENT(in ) :: pdam ! depth of the inflection point 225 REAL(wp), INTENT(in ) :: pwam ! width of inflection 226 REAL(wp), INTENT(in ) :: pbot ! bottom value (0<pbot<= 1) 227 REAL(wp), INTENT(in ), DIMENSION (:) :: pdep ! depth of the gridpoint (T, U, V, F) 228 REAL(wp), INTENT(inout), DIMENSION (:,:,:) :: pah ! adimensional vertical profile 229 !! 230 INTEGER :: jk ! dummy loop indices 231 REAL(wp) :: zm00, zm01, zmhb, zmhs, zcf ! temporary scalars 232 !!---------------------------------------------------------------------- 233 234 zm00 = TANH( ( pdam - gdept_1d(1 ) ) / pwam ) 235 zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam ) 236 zmhs = zm00 / zm01 237 zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01 238 239 DO jk = 1, jpk 240 zcf = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(jk) ) / pwam ) ) 241 pah(:,:,jk) = zcf 242 END DO 243 244 IF(lwp .AND. ld_print ) THEN ! Control print 245 WRITE(numout,*) 246 WRITE(numout,*) ' ahm profile : ' 247 WRITE(numout,*) 248 WRITE(numout,'(" jk ahm "," depth t-level " )') 249 DO jk = 1, jpk 250 WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(jk) 251 END DO 252 ENDIF 253 ! 254 END SUBROUTINE ldf_zpf_1d_3d 255 256 257 SUBROUTINE ldf_zpf_3d( ld_print, pdam, pwam, pbot, pdep, pah ) 258 !!---------------------------------------------------------------------- 259 !! *** ROUTINE ldf_zpf *** 260 !! 261 !! ** Purpose : vertical adimensional profile for eddy coefficient 262 !! 263 !! ** Method : 3D for partial step or s-coordinate 264 !!---------------------------------------------------------------------- 265 LOGICAL , INTENT(in ) :: ld_print ! If true, output arrays on numout 266 REAL(wp), INTENT(in ) :: pdam ! depth of the inflection point 267 REAL(wp), INTENT(in ) :: pwam ! width of inflection 268 REAL(wp), INTENT(in ) :: pbot ! bottom value (0<pbot<= 1) 269 REAL(wp), INTENT(in ), DIMENSION (:,:,:) :: pdep ! dep of the gridpoint (T, U, V, F) 270 REAL(wp), INTENT(inout), DIMENSION (:,:,:) :: pah ! adimensional vertical profile 271 !! 272 INTEGER :: jk ! dummy loop indices 273 REAL(wp) :: zm00, zm01, zmhb, zmhs ! temporary scalars 274 !!---------------------------------------------------------------------- 275 276 zm00 = TANH( ( pdam - gdept_1d(1 ) ) / pwam ) 277 zm01 = TANH( ( pdam - gdept_1d(jpkm1) ) / pwam ) 278 zmhs = zm00 / zm01 279 zmhb = ( 1.e0 - pbot ) / ( 1.e0 - zmhs ) / zm01 280 281 DO jk = 1, jpk 282 pah(:,:,jk) = 1.e0 + zmhb * ( zm00 - TANH( ( pdam - pdep(:,:,jk) ) / pwam ) ) 283 END DO 284 285 IF(lwp .AND. ld_print ) THEN ! Control print 286 WRITE(numout,*) 287 WRITE(numout,*) ' ahm profile : ' 288 WRITE(numout,*) 289 WRITE(numout,'(" jk ahm "," depth t-level " )') 290 DO jk = 1, jpk 291 WRITE(numout,'(i6,2f12.4,3x,2f12.4)') jk, pah(1,1,jk), pdep(1,1,jk) 292 END DO 293 ENDIF 294 ! 295 END SUBROUTINE ldf_zpf_3d 223 224 SUBROUTINE ldf_dyn( kt ) 225 !!---------------------------------------------------------------------- 226 !! *** ROUTINE ldf_dyn *** 227 !! 228 !! ** Purpose : update at kt the momentum lateral mixing coeff. (ahmt and ahmf) 229 !! 230 !! ** Method : time varying eddy viscosity coefficients: 231 !! 232 !! nn_ahm_ijk_t = 31 ahmt, ahmf = F(i,j,k,t) = F(local velocity) 233 !! ( |u|e /12 or |u|e^3/12 for laplacian or bilaplacian operator ) 234 !! BLP case : sqrt of the eddy coef, since bilaplacian is en re-entrant laplacian 235 !! 236 !! ** action : ahmt, ahmf update at each time step 237 !!---------------------------------------------------------------------- 238 INTEGER, INTENT(in) :: kt ! time step index 239 ! 240 INTEGER :: ji, jj, jk ! dummy loop indices 241 REAL(wp) :: zu2pv2_ij_p1, zu2pv2_ij, zu2pv2_ij_m1, zetmax, zefmax ! local scalar 242 !!---------------------------------------------------------------------- 243 ! 244 IF( nn_timing == 1 ) CALL timing_start('ldf_dyn') 245 ! 246 SELECT CASE( nn_ahm_ijk_t ) !== Eddy vicosity coefficients ==! 247 ! 248 CASE( 31 ) !== time varying 3D field ==! = F( local velocity ) 249 ! 250 IF( ln_dynldf_lap ) THEN ! laplacian operator : |u| e /12 = |u/144| e 251 DO jk = 1, jpkm1 252 DO jj = 2, jpjm1 253 DO ji = fs_2, fs_jpim1 254 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 255 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 256 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 257 zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 258 zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 259 ahmt(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax * tmask(ji,jj,jk) ! 288= 12*12 * 2 260 ahmf(ji,jj,jk) = SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax * fmask(ji,jj,jk) 261 END DO 262 END DO 263 END DO 264 ELSEIF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( |u| e^3 /12 ) = sqrt( |u/144| e ) * e 265 DO jk = 1, jpkm1 266 DO jj = 2, jpjm1 267 DO ji = fs_2, fs_jpim1 268 zu2pv2_ij_p1 = ub(ji ,jj+1,jk) * ub(ji ,jj+1,jk) + vb(ji+1,jj ,jk) * vb(ji+1,jj ,jk) 269 zu2pv2_ij = ub(ji ,jj ,jk) * ub(ji ,jj ,jk) + vb(ji ,jj ,jk) * vb(ji ,jj ,jk) 270 zu2pv2_ij_m1 = ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) 271 zetmax = MAX( e1t(ji,jj) , e2t(ji,jj) ) 272 zefmax = MAX( e1f(ji,jj) , e2f(ji,jj) ) 273 ahmt(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_m1) * r1_288 ) * zetmax ) * zetmax * tmask(ji,jj,jk) 274 ahmf(ji,jj,jk) = SQRT( SQRT( (zu2pv2_ij + zu2pv2_ij_p1) * r1_288 ) * zefmax ) * zefmax * fmask(ji,jj,jk) 275 END DO 276 END DO 277 END DO 278 ENDIF 279 ! 280 CALL lbc_lnk( ahmt, 'T', 1. ) ; CALL lbc_lnk( ahmf, 'F', 1. ) 281 ! 282 END SELECT 283 ! 284 CALL iom_put( "ahmt_2d", ahmt(:,:,1) ) ! surface u-eddy diffusivity coeff. 285 CALL iom_put( "ahmf_2d", ahmf(:,:,1) ) ! surface v-eddy diffusivity coeff. 286 CALL iom_put( "ahmt_3d", ahmt(:,:,:) ) ! 3D u-eddy diffusivity coeff. 287 CALL iom_put( "ahmf_3d", ahmf(:,:,:) ) ! 3D v-eddy diffusivity coeff. 288 ! 289 IF( nn_timing == 1 ) CALL timing_stop('ldf_dyn') 290 ! 291 END SUBROUTINE ldf_dyn 296 292 297 293 !!====================================================================== -
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r4488 r6225 11 11 !! 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) add Griffies operator 12 12 !! - ! 2010-11 (F. Dupond, G. Madec) bug correction in slopes just below the ML 13 !! 3.7 ! 2013-12 (F. Lemarie, G. Madec) add limiter on triad slopes 13 14 !!---------------------------------------------------------------------- 14 #if defined key_ldfslp || defined key_esopa 15 15 16 !!---------------------------------------------------------------------- 16 !! 'key_ldfslp' Rotation of lateral mixing tensor17 !!----------------------------------------------------------------------18 !! ldf_slp_grif : calculates the triads of isoneutral slopes (Griffies operator)19 17 !! ldf_slp : calculates the slopes of neutral surface (Madec operator) 18 !! ldf_slp_triad : calculates the triads of isoneutral slopes (Griffies operator) 20 19 !! ldf_slp_mxl : calculates the slopes at the base of the mixed layer (Madec operator) 21 20 !! ldf_slp_init : initialization of the slopes computation … … 23 22 USE oce ! ocean dynamics and tracers 24 23 USE dom_oce ! ocean space and time domain 25 USE ldftra_oce ! lateral diffusion: traceur 26 USE ldfdyn_oce ! lateral diffusion: dynamics 24 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 27 25 USE phycst ! physical constants 28 26 USE zdfmxl ! mixed layer depth 29 27 USE eosbn2 ! equation of states 30 USE lbclnk ! ocean lateral boundary conditions (or mpp link)28 ! 31 29 USE in_out_manager ! I/O manager 32 30 USE prtctl ! Print control 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 32 USE lib_mpp ! distribued memory computing library 33 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 33 34 USE wrk_nemo ! work arrays 34 35 USE timing ! Timing 35 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)36 36 37 37 IMPLICIT NONE 38 38 PRIVATE 39 39 40 PUBLIC ldf_slp ! routine called by step.F90 41 PUBLIC ldf_slp_grif ! routine called by step.F90 42 PUBLIC ldf_slp_init ! routine called by opa.F90 43 44 LOGICAL , PUBLIC, PARAMETER :: lk_ldfslp = .TRUE. !: slopes flag 45 ! !! Madec operator 46 ! Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator 40 PUBLIC ldf_slp ! routine called by step.F90 41 PUBLIC ldf_slp_triad ! routine called by step.F90 42 PUBLIC ldf_slp_init ! routine called by nemogcm.F90 43 44 LOGICAL , PUBLIC :: l_ldfslp = .FALSE. !: slopes flag 45 46 LOGICAL , PUBLIC :: ln_traldf_iso = .TRUE. !: iso-neutral direction 47 LOGICAL , PUBLIC :: ln_traldf_triad = .FALSE. !: griffies triad scheme 48 49 LOGICAL , PUBLIC :: ln_triad_iso = .FALSE. !: pure horizontal mixing in ML 50 LOGICAL , PUBLIC :: ln_botmix_triad = .FALSE. !: mixing on bottom 51 REAL(wp), PUBLIC :: rn_sw_triad = 1._wp !: =1 switching triads ; =0 all four triads used 52 REAL(wp), PUBLIC :: rn_slpmax = 0.01_wp !: slope limit 53 54 LOGICAL , PUBLIC :: l_grad_zps = .FALSE. !: special treatment for Horz Tgradients w partial steps (triad operator) 55 56 ! !! Classic operator (Madec) 47 57 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: uslp, wslpi !: i_slope at U- and W-points 48 58 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vslp, wslpj !: j-slope at V- and W-points 49 ! !! Griffies operator59 ! !! triad operator (Griffies) 50 60 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wslp2 !: wslp**2 from Griffies quarter cells 51 61 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi_g, triadj_g !: skew flux slopes relative to geopotentials 52 62 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi , triadj !: isoneutral slopes relative to model-coordinate 53 54 ! !! Madec operator 55 ! Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator 63 ! !! both operators 64 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ah_wslp2 !: ah * slope^2 at w-point 65 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akz !: stabilizing vertical diffusivity 66 67 ! !! Madec operator 56 68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: omlmask ! mask of the surface mixed layer at T-pt 57 69 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: uslpml, wslpiml ! i_slope at U- and W-points just below the mixed layer … … 61 73 62 74 !! * Substitutions 63 # include "domzgr_substitute.h90"64 # include "ldftra_substitute.h90"65 # include "ldfeiv_substitute.h90"66 75 # include "vectopt_loop_substitute.h90" 67 76 !!---------------------------------------------------------------------- 68 !! NEMO/OPA 4.0 , NEMO Consortium (201 1)77 !! NEMO/OPA 4.0 , NEMO Consortium (2014) 69 78 !! $Id$ 70 79 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 102 111 !! 103 112 INTEGER :: ji , jj , jk ! dummy loop indices 104 INTEGER :: ii0, ii1 , iku! temporary integer105 INTEGER :: ij0, ij1 , ikv! temporary integer106 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16, zcofw ! local scalars113 INTEGER :: ii0, ii1 ! temporary integer 114 INTEGER :: ij0, ij1 ! temporary integer 115 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16, zcofw, z1_slpmax ! local scalars 107 116 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 108 117 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - 109 118 REAL(wp) :: zck, zfk, zbw ! - - 119 REAL(wp) :: zdepu, zdepv ! - - 120 REAL(wp), POINTER, DIMENSION(:,: ) :: zslpml_hmlpu, zslpml_hmlpv 110 121 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zww 111 122 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdzr … … 116 127 ! 117 128 CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 118 119 IF ( ln_traldf_iso .OR. ln_dynldf_iso ) THEN 120 121 zeps = 1.e-20_wp !== Local constant initialization ==! 122 z1_16 = 1.0_wp / 16._wp 123 zm1_g = -1.0_wp / grav 124 zm1_2g = -0.5_wp / grav 125 ! 126 zww(:,:,:) = 0._wp 127 zwz(:,:,:) = 0._wp 128 ! 129 DO jk = 1, jpk !== i- & j-gradient of density ==! 130 DO jj = 1, jpjm1 131 DO ji = 1, fs_jpim1 ! vector opt. 132 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 133 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 134 END DO 135 END DO 136 END DO 137 IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level 138 # if defined key_vectopt_loop 139 DO jj = 1, 1 140 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 141 # else 142 DO jj = 1, jpjm1 143 DO ji = 1, jpim1 144 # endif 145 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 146 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 147 END DO 148 END DO 149 ENDIF 150 ! 151 zdzr(:,:,1) = 0._wp !== Local vertical density gradient at T-point == ! (evaluated from N^2) 152 DO jk = 2, jpkm1 153 ! ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 154 ! ! trick: tmask(ik ) = 0 => all pn2 = 0 => zdzr = 0 155 ! ! else tmask(ik+1) = 0 => pn2(ik+1) = 0 => zdzr divides by 1 156 ! ! umask(ik+1) /= 0 => all pn2 /= 0 => zdzr divides by 2 157 ! ! NB: 1/(tmask+1) = (1-.5*tmask) substitute a / by a * ==> faster 158 zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp ) & 159 & * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 160 END DO 161 ! 162 ! !== Slopes just below the mixed layer ==! 163 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml 164 165 166 ! I. slopes at u and v point | uslp = d/di( prd ) / d/dz( prd ) 167 ! =========================== | vslp = d/dj( prd ) / d/dz( prd ) 168 ! 169 DO jk = 2, jpkm1 !* Slopes at u and v points 170 DO jj = 2, jpjm1 171 DO ji = fs_2, fs_jpim1 ! vector opt. 172 ! ! horizontal and vertical density gradient at u- and v-points 173 zau = zgru(ji,jj,jk) / e1u(ji,jj) 174 zav = zgrv(ji,jj,jk) / e2v(ji,jj) 175 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 176 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 177 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 178 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 179 zbu = MIN( zbu, -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau ) ) 180 zbv = MIN( zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav ) ) 181 ! ! uslp and vslp output in zwz and zww, resp. 182 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 183 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 184 zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps ) & 185 & + zfi * uslpml(ji,jj) & 186 & * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) ) & 187 & / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk) 188 zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps ) & 189 & + zfj * vslpml(ji,jj) & 190 & * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) ) & 191 & / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 129 CALL wrk_alloc( jpi,jpj, zslpml_hmlpu, zslpml_hmlpv ) 130 131 zeps = 1.e-20_wp !== Local constant initialization ==! 132 z1_16 = 1.0_wp / 16._wp 133 zm1_g = -1.0_wp / grav 134 zm1_2g = -0.5_wp / grav 135 z1_slpmax = 1._wp / rn_slpmax 136 ! 137 zww(:,:,:) = 0._wp 138 zwz(:,:,:) = 0._wp 139 ! 140 DO jk = 1, jpk !== i- & j-gradient of density ==! 141 DO jj = 1, jpjm1 142 DO ji = 1, fs_jpim1 ! vector opt. 143 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 144 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 145 END DO 146 END DO 147 END DO 148 IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level 149 DO jj = 1, jpjm1 150 DO ji = 1, jpim1 151 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 152 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 153 END DO 154 END DO 155 ENDIF 156 IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the bottom ocean level 157 DO jj = 1, jpjm1 158 DO ji = 1, jpim1 159 IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 160 IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 161 END DO 162 END DO 163 ENDIF 164 ! 165 zdzr(:,:,1) = 0._wp !== Local vertical density gradient at T-point == ! (evaluated from N^2) 166 DO jk = 2, jpkm1 167 ! ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 168 ! ! trick: tmask(ik ) = 0 => all pn2 = 0 => zdzr = 0 169 ! ! else tmask(ik+1) = 0 => pn2(ik+1) = 0 => zdzr divides by 1 170 ! ! umask(ik+1) /= 0 => all pn2 /= 0 => zdzr divides by 2 171 ! ! NB: 1/(tmask+1) = (1-.5*tmask) substitute a / by a * ==> faster 172 zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp ) & 173 & * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 174 END DO 175 ! 176 ! !== Slopes just below the mixed layer ==! 177 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml 178 179 180 ! I. slopes at u and v point | uslp = d/di( prd ) / d/dz( prd ) 181 ! =========================== | vslp = d/dj( prd ) / d/dz( prd ) 182 ! 183 IF ( ln_isfcav ) THEN 184 DO jj = 2, jpjm1 185 DO ji = fs_2, fs_jpim1 ! vector opt. 186 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) & 187 & - 0.5_wp * ( risfdep(ji,jj) + risfdep(ji+1,jj ) ) ) 188 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) & 189 & - 0.5_wp * ( risfdep(ji,jj) + risfdep(ji ,jj+1) ) ) 190 END DO 191 END DO 192 ELSE 193 DO jj = 2, jpjm1 194 DO ji = fs_2, fs_jpim1 ! vector opt. 195 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) 196 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) 197 END DO 198 END DO 199 END IF 200 201 DO jk = 2, jpkm1 !* Slopes at u and v points 202 DO jj = 2, jpjm1 203 DO ji = fs_2, fs_jpim1 ! vector opt. 204 ! ! horizontal and vertical density gradient at u- and v-points 205 zau = zgru(ji,jj,jk) * r1_e1u(ji,jj) 206 zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj) 207 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 208 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 209 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 210 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 211 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,jk)* ABS( zau ) ) 212 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,jk)* ABS( zav ) ) 213 ! ! uslp and vslp output in zwz and zww, resp. 214 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 215 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 216 ! thickness of water column between surface and level k at u/v point 217 zdepu = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji+1,jj ,jk) ) & 218 - ( risfdep(ji,jj) + risfdep(ji+1,jj) ) - e3u_n(ji,jj,miku(ji,jj)) ) 219 zdepv = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji,jj+1,jk) ) & 220 - ( risfdep(ji,jj) + risfdep(ji,jj+1) ) - e3v_n(ji,jj,mikv(ji,jj)) ) 221 ! 222 zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps ) & 223 & + zfi * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk) 224 zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps ) & 225 & + zfj * zdepv * zslpml_hmlpv(ji,jj) ) * vmask(ji,jj,jk) 192 226 !!gm modif to suppress omlmask.... (as in Griffies case) 193 ! 194 ! 195 ! 196 ! zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )197 ! zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )198 ! 199 ! 227 ! ! ! jk must be >= ML level for zf=1. otherwise zf=0. 228 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 229 ! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 230 ! zci = 0.5 * ( gdept_n(ji+1,jj,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 231 ! zcj = 0.5 * ( gdept_n(ji,jj+1,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 232 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 233 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 200 234 !!gm end modif 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk)&232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) *e2t(ji,jj)262 263 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * tmask (ji,jj,jk)264 265 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * tmask (ji,jj,jk)266 267 268 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai ) )269 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj ) )270 271 272 zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp )273 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * tmask(ji,jj,jk)274 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * tmask(ji,jj,jk)235 END DO 236 END DO 237 END DO 238 CALL lbc_lnk( zwz, 'U', -1. ) ; CALL lbc_lnk( zww, 'V', -1. ) ! lateral boundary conditions 239 ! 240 ! !* horizontal Shapiro filter 241 DO jk = 2, jpkm1 242 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 243 DO ji = 2, jpim1 244 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 245 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 246 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 247 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 248 & + 4.* zwz(ji ,jj ,jk) ) 249 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 250 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 251 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 252 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 253 & + 4.* zww(ji,jj ,jk) ) 254 END DO 255 END DO 256 DO jj = 3, jpj-2 ! other rows 257 DO ji = fs_2, fs_jpim1 ! vector opt. 258 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 259 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 260 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 261 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 262 & + 4.* zwz(ji ,jj ,jk) ) 263 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 264 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 265 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 266 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 267 & + 4.* zww(ji,jj ,jk) ) 268 END DO 269 END DO 270 ! !* decrease along coastal boundaries 271 DO jj = 2, jpjm1 272 DO ji = fs_2, fs_jpim1 ! vector opt. 273 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 274 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 275 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 276 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 277 END DO 278 END DO 279 END DO 280 281 282 ! II. slopes at w point | wslpi = mij( d/di( prd ) / d/dz( prd ) 283 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 284 ! 285 DO jk = 2, jpkm1 286 DO jj = 2, jpjm1 287 DO ji = fs_2, fs_jpim1 ! vector opt. 288 ! !* Local vertical density gradient evaluated from N^2 289 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 290 ! !* Slopes at w point 291 ! ! i- & j-gradient of density at w-points 292 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 293 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 294 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 295 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 296 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 297 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * wmask (ji,jj,jk) 298 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 299 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * wmask (ji,jj,jk) 300 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 301 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 302 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zai ) ) 303 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zaj ) ) 304 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 305 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 306 zck = ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj) - gdepw_n(ji,jj,mikt(ji,jj)), 10._wp ) 307 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 308 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk) 275 309 276 310 !!gm modif to suppress omlmask.... (as in Griffies operator) 277 ! 278 ! 279 ! zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. )280 ! 281 ! 311 ! ! ! jk must be >= ML level for zfk=1. otherwise zfk=0. 312 ! zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp ) 313 ! zck = gdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. ) 314 ! zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk) 315 ! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 282 316 !!gm end modif 283 END DO 284 END DO 285 END DO 286 CALL lbc_lnk( zwz, 'T', -1. ) ; CALL lbc_lnk( zww, 'T', -1. ) ! lateral boundary conditions 287 ! 288 ! !* horizontal Shapiro filter 289 DO jk = 2, jpkm1 290 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 291 DO ji = 2, jpim1 292 zcofw = tmask(ji,jj,jk) * z1_16 293 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 294 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 295 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 296 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 297 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 298 299 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 300 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 301 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 302 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 303 & + 4.* zww(ji ,jj ,jk) ) * zcofw 304 END DO 305 END DO 306 DO jj = 3, jpj-2 ! other rows 307 DO ji = fs_2, fs_jpim1 ! vector opt. 308 zcofw = tmask(ji,jj,jk) * z1_16 309 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 310 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 311 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 312 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 313 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 314 315 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 316 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 317 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 318 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 319 & + 4.* zww(ji ,jj ,jk) ) * zcofw 320 END DO 321 END DO 322 ! !* decrease along coastal boundaries 323 DO jj = 2, jpjm1 324 DO ji = fs_2, fs_jpim1 ! vector opt. 325 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 326 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 327 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 328 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 329 END DO 330 END DO 331 END DO 332 333 ! III. Specific grid points 334 ! =========================== 335 ! 336 IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN ! ORCA_R4 configuration: horizontal diffusion in specific area 337 ! ! Gibraltar Strait 338 ij0 = 50 ; ij1 = 53 339 ii0 = 69 ; ii1 = 71 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 340 ij0 = 51 ; ij1 = 53 341 ii0 = 68 ; ii1 = 71 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 342 ii0 = 69 ; ii1 = 71 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 343 ii0 = 69 ; ii1 = 71 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 344 ! 345 ! ! Mediterrannean Sea 346 ij0 = 49 ; ij1 = 56 347 ii0 = 71 ; ii1 = 90 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 348 ij0 = 50 ; ij1 = 56 349 ii0 = 70 ; ii1 = 90 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 350 ii0 = 71 ; ii1 = 90 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 351 ii0 = 71 ; ii1 = 90 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 352 ENDIF 353 354 355 ! IV. Lateral boundary conditions 356 ! =============================== 357 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 358 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 359 360 361 IF(ln_ctl) THEN 362 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 363 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 364 ENDIF 365 ! 366 367 ELSEIF ( lk_vvl ) THEN 368 369 IF(lwp) THEN 370 WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 371 ENDIF 372 373 ! geopotential diffusion in s-coordinates on tracers and/or momentum 374 ! The slopes of s-surfaces are computed at each time step due to vvl 375 ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 376 377 ! set the slope of diffusion to the slope of s-surfaces 378 ! ( c a u t i o n : minus sign as fsdep has positive value ) 379 DO jk = 1, jpk 380 DO jj = 2, jpjm1 381 DO ji = fs_2, fs_jpim1 ! vector opt. 382 uslp(ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 383 vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 384 wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 385 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 386 END DO 387 END DO 388 END DO 389 390 ! Lateral boundary conditions on the slopes 391 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 392 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 393 394 if( kt == nit000 ) then 395 IF(lwp) WRITE(numout,*) ' max slop: u',SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)), & 396 & ' wi', sqrt(MAXVAL(wslpi)), ' wj', sqrt(MAXVAL(wslpj)) 397 endif 398 399 IF(ln_ctl) THEN 400 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 401 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 402 ENDIF 403 317 END DO 318 END DO 319 END DO 320 CALL lbc_lnk( zwz, 'T', -1. ) ; CALL lbc_lnk( zww, 'T', -1. ) ! lateral boundary conditions 321 ! 322 ! !* horizontal Shapiro filter 323 DO jk = 2, jpkm1 324 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 325 DO ji = 2, jpim1 326 zcofw = wmask(ji,jj,jk) * z1_16 327 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 328 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 329 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 330 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 331 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 332 333 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 334 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 335 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 336 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 337 & + 4.* zww(ji ,jj ,jk) ) * zcofw 338 END DO 339 END DO 340 DO jj = 3, jpj-2 ! other rows 341 DO ji = fs_2, fs_jpim1 ! vector opt. 342 zcofw = wmask(ji,jj,jk) * z1_16 343 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 344 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 345 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 346 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 347 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 348 349 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 350 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 351 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 352 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 353 & + 4.* zww(ji ,jj ,jk) ) * zcofw 354 END DO 355 END DO 356 ! !* decrease in vicinity of topography 357 DO jj = 2, jpjm1 358 DO ji = fs_2, fs_jpim1 ! vector opt. 359 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 360 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 361 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 362 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 363 END DO 364 END DO 365 END DO 366 367 ! IV. Lateral boundary conditions 368 ! =============================== 369 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 370 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 371 372 IF(ln_ctl) THEN 373 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 374 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 404 375 ENDIF 405 376 ! 406 377 CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 378 CALL wrk_dealloc( jpi,jpj, zslpml_hmlpu, zslpml_hmlpv ) 407 379 ! 408 380 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp') … … 411 383 412 384 413 SUBROUTINE ldf_slp_ grif( kt )414 !!---------------------------------------------------------------------- 415 !! *** ROUTINE ldf_slp_ grif***385 SUBROUTINE ldf_slp_triad ( kt ) 386 !!---------------------------------------------------------------------- 387 !! *** ROUTINE ldf_slp_triad *** 416 388 !! 417 389 !! ** Purpose : Compute the squared slopes of neutral surfaces (slope 418 !! of iso-pycnal surfaces referenced locally) (ln_traldf_ grif=T)390 !! of iso-pycnal surfaces referenced locally) (ln_traldf_triad=T) 419 391 !! at W-points using the Griffies quarter-cells. 420 392 !! … … 431 403 REAL(wp) :: zfacti, zfactj ! local scalars 432 404 REAL(wp) :: znot_thru_surface ! local scalars 433 REAL(wp) :: zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj 405 REAL(wp) :: zdit, zdis, zdkt, zbu, zbti, zisw 406 REAL(wp) :: zdjt, zdjs, zdks, zbv, zbtj, zjsw 434 407 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_g_raw, zti_g_lim 435 408 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim 436 409 REAL(wp) :: zdzrho_raw 437 REAL(wp) :: zbeta0 410 REAL(wp) :: zbeta0, ze3_e1, ze3_e2 438 411 REAL(wp), POINTER, DIMENSION(:,:) :: z1_mlbw 439 412 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalbet … … 442 415 !!---------------------------------------------------------------------- 443 416 ! 444 IF( nn_timing == 1 ) CALL timing_start('ldf_slp_ grif')417 IF( nn_timing == 1 ) CALL timing_start('ldf_slp_triad') 445 418 ! 446 419 CALL wrk_alloc( jpi,jpj, z1_mlbw ) … … 452 425 ! Some preliminary calculation ! 453 426 !--------------------------------! 454 !455 CALL eos_alpbet( tsb, zalbet, zbeta0 ) !== before local thermal/haline expension ratio at T-points ==!456 427 ! 457 428 DO jl = 0, 1 !== unmasked before density i- j-, k-gradients ==! … … 465 436 zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! j-gradient of T & S at v-point 466 437 zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 467 zdxrho_raw = ( - zalbet(ji+ip,jj ,jk) * zdit + zbeta0*zdis ) /e1u(ji,jj)468 zdyrho_raw = ( - zalbet(ji ,jj+jp,jk) * zdjt + zbeta0*zdjs ) /e2v(ji,jj)469 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw) ! keep the sign470 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw)438 zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj) 439 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj) 440 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 441 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 471 442 END DO 472 443 END DO 473 444 END DO 474 445 ! 475 IF( ln_zps.and.l_grad_zps ) THEN ! partial steps: correction of i- & j-grad on bottom 476 # if defined key_vectopt_loop 477 DO jj = 1, 1 478 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 479 # else 446 IF( ln_zps .AND. l_grad_zps ) THEN ! partial steps: correction of i- & j-grad on bottom 480 447 DO jj = 1, jpjm1 481 448 DO ji = 1, jpim1 482 # endif483 449 iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points) 484 450 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 485 451 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 486 zdxrho_raw = ( - zalbet(ji+ip,jj ,iku) * zdit + zbeta0*zdis ) /e1u(ji,jj)487 zdyrho_raw = ( - zalbet(ji ,jj+jp,ikv) * zdjt + zbeta0*zdjs ) /e2v(ji,jj)452 zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj) 453 zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj) 488 454 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 489 455 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) … … 505 471 zdks = 0._wp 506 472 ENDIF 507 zdzrho_raw = ( - zalbet(ji ,jj ,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp)508 zdzrho(ji ,jj ,jk, kp) = - MIN( - repsln,zdzrho_raw ) ! force zdzrho >= repsln473 zdzrho_raw = ( - zalbet(ji,jj,jk) * zdkt + zbeta0*zdks ) / e3w_n(ji,jj,jk+kp) 474 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln 509 475 END DO 510 476 END DO … … 515 481 DO ji = 1, jpi 516 482 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth 517 z1_mlbw(ji,jj) = 1._wp / fsdepw(ji,jj,jk)483 z1_mlbw(ji,jj) = 1._wp / gdepw_n(ji,jj,jk) 518 484 END DO 519 485 END DO … … 539 505 ! 540 506 jk = nmln(ji+ip,jj) + 1 541 IF( jk .GT. mbkt(ji+ip,jj) ) THEN !ML reaches bottom 542 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 543 ELSE 544 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 545 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 546 & - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj) ) * umask(ji,jj,jk) 547 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 507 IF( jk > mbkt(ji+ip,jj) ) THEN ! ML reaches bottom 508 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 509 ELSE 510 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 511 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 512 & - ( gdept_n(ji+1,jj,jk-kp) - gdept_n(ji,jj,jk-kp) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk) 513 ze3_e1 = e3w_n(ji+ip,jj,jk-kp) * r1_e1u(ji,jj) 514 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1 , ABS( zti_g_raw ) ), zti_g_raw ) 548 515 ENDIF 549 516 ! 550 517 jk = nmln(ji,jj+jp) + 1 551 IF( jk .GT.mbkt(ji,jj+jp) ) THEN !ML reaches bottom552 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp518 IF( jk > mbkt(ji,jj+jp) ) THEN !ML reaches bottom 519 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp 553 520 ELSE 554 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 555 & - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 556 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 521 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 522 & - ( gdept_n(ji,jj+1,jk-kp) - gdept_n(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 523 ze3_e2 = e3w_n(ji,jj+jp,jk-kp) / e2v(ji,jj) 524 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2 , ABS( ztj_g_raw ) ), ztj_g_raw ) 557 525 ENDIF 558 526 END DO … … 583 551 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 584 552 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 585 553 ! 586 554 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 587 zti_coord = znot_thru_surface * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) /e1u(ji,jj)588 ztj_coord = znot_thru_surface * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)! unmasked555 zti_coord = znot_thru_surface * ( gdept_n(ji+1,jj ,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) 556 ztj_coord = znot_thru_surface * ( gdept_n(ji ,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) ! unmasked 589 557 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 590 558 ztj_g_raw = ztj_raw - ztj_coord 591 zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 592 ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 559 ! additional limit required in bilaplacian case 560 ze3_e1 = e3w_n(ji+ip,jj ,jk+kp) * r1_e1u(ji,jj) 561 ze3_e2 = e3w_n(ji ,jj+jp,jk+kp) * r1_e2v(ji,jj) 562 ! NB: hard coded factor 5 (can be a namelist parameter...) 563 zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) 564 ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw ) 593 565 ! 594 566 ! Below ML use limited zti_g as is & mask … … 600 572 zti_g_lim = ( zfacti * zti_g_lim & 601 573 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 602 & * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp)574 & * gdepw_n(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 603 575 ztj_g_lim = ( zfactj * ztj_g_lim & 604 576 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 605 & * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp)577 & * gdepw_n(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 606 578 ! 607 579 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim … … 619 591 ! 620 592 IF( ln_triad_iso ) THEN 621 zti_raw = zti_lim* *2/ zti_raw622 ztj_raw = ztj_lim* *2/ ztj_raw593 zti_raw = zti_lim*zti_lim / zti_raw 594 ztj_raw = ztj_lim*ztj_lim / ztj_raw 623 595 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 624 596 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 625 zti_lim = zfacti * zti_lim & 626 & + ( 1._wp - zfacti ) * zti_raw 627 ztj_lim = zfactj * ztj_lim & 628 & + ( 1._wp - zfactj ) * ztj_raw 597 zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw 598 ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw 629 599 ENDIF 630 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim 631 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim 632 ! 633 zbu = e1u(ji ,jj) * e2u(ji ,jj) * fse3u(ji ,jj,jk ) 634 zbv = e1v(ji ,jj) * e2v(ji ,jj) * fse3v(ji ,jj,jk ) 635 zbti = e1t(ji+ip,jj) * e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp) 636 zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 637 ! 638 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked 639 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim**2 ! masked 640 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim**2 600 ! ! switching triad scheme 601 zisw = (rn_sw_triad - 1._wp ) + rn_sw_triad & 602 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) ) ) 603 zjsw = (rn_sw_triad - 1._wp ) + rn_sw_triad & 604 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) ) ) 605 ! 606 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim * zisw 607 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 608 ! 609 zbu = e1e2u(ji ,jj ) * e3u_n(ji ,jj ,jk ) 610 zbv = e1e2v(ji ,jj ) * e3v_n(ji ,jj ,jk ) 611 zbti = e1e2t(ji+ip,jj ) * e3w_n(ji+ip,jj ,jk+kp) 612 zbtj = e1e2t(ji ,jj+jp) * e3w_n(ji ,jj+jp,jk+kp) 613 ! 614 wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked 615 wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim 641 616 END DO 642 617 END DO … … 654 629 CALL wrk_dealloc( jpi,jpj, 2,2, zti_mlb, ztj_mlb, kkstart = 0, klstart = 0 ) 655 630 ! 656 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_ grif')657 ! 658 END SUBROUTINE ldf_slp_ grif631 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_triad') 632 ! 633 END SUBROUTINE ldf_slp_triad 659 634 660 635 … … 682 657 INTEGER :: ji , jj , jk ! dummy loop indices 683 658 INTEGER :: iku, ikv, ik, ikm1 ! local integers 684 REAL(wp) :: zeps, zm1_g, zm1_2g 659 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_slpmax ! local scalars 685 660 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 686 661 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - … … 693 668 zm1_g = -1.0_wp / grav 694 669 zm1_2g = -0.5_wp / grav 670 z1_slpmax = 1._wp / rn_slpmax 695 671 ! 696 672 uslpml (1,:) = 0._wp ; uslpml (jpi,:) = 0._wp … … 701 677 ! !== surface mixed layer mask ! 702 678 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise 703 # if defined key_vectopt_loop704 DO jj = 1, 1705 DO ji = 1, jpij ! vector opt. (forced unrolling)706 # else707 679 DO jj = 1, jpj 708 680 DO ji = 1, jpi 709 # endif710 681 ik = nmln(ji,jj) - 1 711 682 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp … … 727 698 !----------------------------------------------------------------------- 728 699 ! 729 # if defined key_vectopt_loop730 DO jj = 1, 1731 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)732 # else733 700 DO jj = 2, jpjm1 734 701 DO ji = 2, jpim1 735 # endif736 702 ! !== Slope at u- & v-points just below the Mixed Layer ==! 737 703 ! … … 742 708 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 743 709 ! !- horizontal density gradient at u- & v-points 744 zau = p_gru(ji,jj,iku) /e1u(ji,jj)745 zav = p_grv(ji,jj,ikv) /e2v(ji,jj)710 zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj) 711 zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj) 746 712 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 747 713 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 748 zbu = MIN( zbu , - 100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) )749 zbv = MIN( zbv , - 100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) )714 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,iku)* ABS( zau ) ) 715 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,ikv)* ABS( zav ) ) 750 716 ! !- Slope at u- & v-points (uslpml, vslpml) 751 717 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) … … 763 729 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 764 730 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) 765 zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) &731 zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) & 766 732 & + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1 ) ) / zci * tmask(ji,jj,ik) 767 733 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & … … 769 735 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 770 736 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 771 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/ fse3w(ji,jj,ik)* ABS( zai ) )772 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/ fse3w(ji,jj,ik)* ABS( zaj ) )737 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zai ) ) 738 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zaj ) ) 773 739 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 774 740 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) … … 791 757 !! ** Purpose : Initialization for the isopycnal slopes computation 792 758 !! 793 !! ** Method : read the nammbf namelist and check the parameter 794 !! values called by tra_dmp at the first timestep (nit000) 759 !! ** Method : 795 760 !!---------------------------------------------------------------------- 796 761 INTEGER :: ji, jj, jk ! dummy loop indices … … 805 770 WRITE(numout,*) '~~~~~~~~~~~~' 806 771 ENDIF 807 808 IF( ln_traldf_grif ) THEN ! Griffies operator : triad of slopes 809 ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , wslp2(jpi,jpj,jpk) , STAT=ierr ) 810 ALLOCATE( triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , STAT=ierr ) 811 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' ) 812 ! 772 ! 773 ALLOCATE( ah_wslp2(jpi,jpj,jpk) , akz(jpi,jpj,jpk) , STAT=ierr ) 774 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate ah_slp2 or akz' ) 775 ! 776 IF( ln_traldf_triad ) THEN ! Griffies operator : triad of slopes 777 IF(lwp) WRITE(numout,*) ' Griffies (triad) operator initialisation' 778 ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , & 779 & triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , & 780 & wslp2 (jpi,jpj,jpk) , STAT=ierr ) 781 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' ) 813 782 IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 814 783 ! 815 784 ELSE ! Madec operator : slopes at u-, v-, and w-points 816 ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) , & 817 & omlmask(jpi,jpj,jpk) , uslpml(jpi,jpj) , vslpml(jpi,jpj) , wslpiml(jpi,jpj) , wslpjml(jpi,jpj) , STAT=ierr ) 785 IF(lwp) WRITE(numout,*) ' Madec operator initialisation' 786 ALLOCATE( omlmask(jpi,jpj,jpk) , & 787 & uslp(jpi,jpj,jpk) , uslpml(jpi,jpj) , wslpi(jpi,jpj,jpk) , wslpiml(jpi,jpj) , & 788 & vslp(jpi,jpj,jpk) , vslpml(jpi,jpj) , wslpj(jpi,jpj,jpk) , wslpjml(jpi,jpj) , STAT=ierr ) 818 789 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Madec operator slope ' ) 819 790 … … 825 796 wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp 826 797 827 IF( ln_traldf_hor .OR. ln_dynldf_hor ) THEN 828 IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 829 830 ! geopotential diffusion in s-coordinates on tracers and/or momentum 831 ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 832 ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 833 834 ! set the slope of diffusion to the slope of s-surfaces 835 ! ( c a u t i o n : minus sign as fsdep has positive value ) 836 DO jk = 1, jpk 837 DO jj = 2, jpjm1 838 DO ji = fs_2, fs_jpim1 ! vector opt. 839 uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk) 840 vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk) 841 wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5 842 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5 843 END DO 844 END DO 845 END DO 846 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) ! Lateral boundary conditions 847 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 848 ENDIF 798 !!gm I no longer understand this..... 799 !!gm IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (.NOT.ln_linssh .AND. ln_rstart) ) THEN 800 ! IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 801 ! 802 ! ! geopotential diffusion in s-coordinates on tracers and/or momentum 803 ! ! The slopes of s-surfaces are computed once (no call to ldfslp in step) 804 ! ! The slopes for momentum diffusion are i- or j- averaged of those on tracers 805 ! 806 ! ! set the slope of diffusion to the slope of s-surfaces 807 ! ! ( c a u t i o n : minus sign as dep has positive value ) 808 ! DO jk = 1, jpk 809 ! DO jj = 2, jpjm1 810 ! DO ji = fs_2, fs_jpim1 ! vector opt. 811 ! uslp (ji,jj,jk) = - ( gdept_n(ji+1,jj,jk) - gdept_n(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 812 ! vslp (ji,jj,jk) = - ( gdept_n(ji,jj+1,jk) - gdept_n(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 813 ! wslpi(ji,jj,jk) = - ( gdepw_n(ji+1,jj,jk) - gdepw_n(ji-1,jj,jk) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5 814 ! wslpj(ji,jj,jk) = - ( gdepw_n(ji,jj+1,jk) - gdepw_n(ji,jj-1,jk) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5 815 ! END DO 816 ! END DO 817 ! END DO 818 ! CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) ! Lateral boundary conditions 819 ! CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 820 !!gm ENDIF 849 821 ENDIF 850 822 ! … … 852 824 ! 853 825 END SUBROUTINE ldf_slp_init 854 855 #else856 !!------------------------------------------------------------------------857 !! Dummy module : NO Rotation of lateral mixing tensor858 !!------------------------------------------------------------------------859 LOGICAL, PUBLIC, PARAMETER :: lk_ldfslp = .FALSE. !: slopes flag860 CONTAINS861 SUBROUTINE ldf_slp( kt, prd, pn2 ) ! Dummy routine862 INTEGER, INTENT(in) :: kt863 REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2864 WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)865 END SUBROUTINE ldf_slp866 SUBROUTINE ldf_slp_grif( kt ) ! Dummy routine867 INTEGER, INTENT(in) :: kt868 WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt869 END SUBROUTINE ldf_slp_grif870 SUBROUTINE ldf_slp_init ! Dummy routine871 END SUBROUTINE ldf_slp_init872 #endif873 826 874 827 !!====================================================================== -
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
r4624 r6225 2 2 !!====================================================================== 3 3 !! *** MODULE ldftra *** 4 !! Ocean physics: lateral diffusivity coefficient 4 !! Ocean physics: lateral diffusivity coefficients 5 5 !!===================================================================== 6 !! History : ! 1997-07 (G. Madec) from inimix.F split in 2 routines 7 !! NEMO 1.0 ! 2002-09 (G. Madec) F90: Free form and module 8 !! 2.0 ! 2005-11 (G. Madec) 6 !! History : ! 1997-07 (G. Madec) from inimix.F split in 2 routines 7 !! NEMO 1.0 ! 2002-09 (G. Madec) F90: Free form and module 8 !! 2.0 ! 2005-11 (G. Madec) 9 !! 3.7 ! 2013-12 (F. Lemarie, G. Madec) restructuration/simplification of aht/aeiv specification, 10 !! ! add velocity dependent coefficient and optional read in file 9 11 !!---------------------------------------------------------------------- 10 12 11 13 !!---------------------------------------------------------------------- 12 14 !! ldf_tra_init : initialization, namelist read, and parameters control 13 !! ldf_tra_c3d : 3D eddy viscosity coefficient initialization 14 !! ldf_tra_c2d : 2D eddy viscosity coefficient initialization 15 !! ldf_tra_c1d : 1D eddy viscosity coefficient initialization 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 !! ldf_eiv : time evolution of the eiv coefficients (function of the growth rate of baroclinic instability) 18 !! ldf_eiv_trp : add to the input ocean transport the contribution of the EIV parametrization 19 !! ldf_eiv_dia : diagnose the eddy induced velocity from the eiv streamfunction 16 20 !!---------------------------------------------------------------------- 17 21 USE oce ! ocean dynamics and tracers 18 22 USE dom_oce ! ocean space and time domain 19 23 USE phycst ! physical constants 20 USE ldftra_oce ! ocean tracer lateral physics 21 USE ldfslp ! ??? 24 USE ldfslp ! lateral diffusion: slope of iso-neutral surfaces 25 USE ldfc1d_c2d ! lateral diffusion: 1D & 2D cases 26 USE diaar5, ONLY: lk_diaar5 27 ! 28 USE trc_oce, ONLY: lk_offline ! offline flag 22 29 USE in_out_manager ! I/O manager 23 USE io ipsl30 USE iom ! I/O module for ehanced bottom friction file 24 31 USE lib_mpp ! distribued memory computing library 25 32 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 33 USE wrk_nemo ! work arrays 34 USE timing ! timing 26 35 27 36 IMPLICIT NONE 28 37 PRIVATE 29 38 30 PUBLIC ldf_tra_init ! called by opa.F90 39 PUBLIC ldf_tra_init ! called by nemogcm.F90 40 PUBLIC ldf_tra ! called by step.F90 41 PUBLIC ldf_eiv_init ! called by nemogcm.F90 42 PUBLIC ldf_eiv ! called by step.F90 43 PUBLIC ldf_eiv_trp ! called by traadv.F90 44 PUBLIC ldf_eiv_dia ! called by traldf_iso and traldf_iso_triad.F90 45 46 ! !!* Namelist namtra_ldf : lateral mixing on tracers * 47 ! != Operator type =! 48 LOGICAL , PUBLIC :: ln_traldf_lap !: laplacian operator 49 LOGICAL , PUBLIC :: ln_traldf_blp !: bilaplacian operator 50 ! != Direction of action =! 51 LOGICAL , PUBLIC :: ln_traldf_lev !: iso-level direction 52 LOGICAL , PUBLIC :: ln_traldf_hor !: horizontal (geopotential) direction 53 ! LOGICAL , PUBLIC :: ln_traldf_iso !: iso-neutral direction (see ldfslp) 54 ! LOGICAL , PUBLIC :: ln_traldf_triad !: griffies triad scheme (see ldfslp) 55 LOGICAL , PUBLIC :: ln_traldf_msc !: Method of Stabilizing Correction 56 ! LOGICAL , PUBLIC :: ln_triad_iso !: pure horizontal mixing in ML (see ldfslp) 57 ! LOGICAL , PUBLIC :: ln_botmix_triad !: mixing on bottom (see ldfslp) 58 ! REAL(wp), PUBLIC :: rn_sw_triad !: =1/0 switching triad / all 4 triads used (see ldfslp) 59 ! REAL(wp), PUBLIC :: rn_slpmax !: slope limit (see ldfslp) 60 ! != Coefficients =! 61 INTEGER , PUBLIC :: nn_aht_ijk_t !: choice of time & space variations of the lateral eddy diffusivity coef. 62 REAL(wp), PUBLIC :: rn_aht_0 !: laplacian lateral eddy diffusivity [m2/s] 63 REAL(wp), PUBLIC :: rn_bht_0 !: bilaplacian lateral eddy diffusivity [m4/s] 64 65 ! !!* Namelist namtra_ldfeiv : eddy induced velocity param. * 66 ! != Use/diagnose eiv =! 67 LOGICAL , PUBLIC :: ln_ldfeiv !: eddy induced velocity flag 68 LOGICAL , PUBLIC :: ln_ldfeiv_dia !: diagnose & output eiv streamfunction and velocity (IOM) 69 ! != Coefficients =! 70 INTEGER , PUBLIC :: nn_aei_ijk_t !: choice of time/space variation of the eiv coeff. 71 REAL(wp), PUBLIC :: rn_aeiv_0 !: eddy induced velocity coefficient [m2/s] 72 73 LOGICAL , PUBLIC :: l_ldftra_time = .FALSE. !: flag for time variation of the lateral eddy diffusivity coef. 74 LOGICAL , PUBLIC :: l_ldfeiv_time = .FALSE. ! flag for time variation of the eiv coef. 75 76 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahtu, ahtv !: eddy diffusivity coef. at U- and V-points [m2/s] 77 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aeiu, aeiv !: eddy induced velocity coeff. [m2/s] 78 79 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 80 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 31 81 32 82 !! * Substitutions 33 # include "domzgr_substitute.h90"34 83 # include "vectopt_loop_substitute.h90" 35 84 !!---------------------------------------------------------------------- 36 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)85 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 37 86 !! $Id$ 38 87 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 46 95 !! ** Purpose : initializations of the tracer lateral mixing coeff. 47 96 !! 48 !! ** Method : the Eddy diffusivity and eddy induced velocity ceoff. 49 !! are defined as follows: 50 !! default option : constant coef. aht0, aeiv0 (namelist) 51 !! 'key_traldf_c1d': depth dependent coef. defined in 52 !! in ldf_tra_c1d routine 53 !! 'key_traldf_c2d': latitude and longitude dependent coef. 54 !! defined in ldf_tra_c2d routine 55 !! 'key_traldf_c3d': latitude, longitude, depth dependent coef. 56 !! defined in ldf_tra_c3d routine 57 !! 58 !! N.B. User defined include files. By default, 3d and 2d coef. 59 !! are set to a constant value given in the namelist and the 1d 60 !! coefficients are initialized to a hyperbolic tangent vertical 61 !! profile. 62 !!---------------------------------------------------------------------- 63 INTEGER :: ioptio ! temporary integer 64 INTEGER :: ios ! temporary integer 65 LOGICAL :: ll_print = .FALSE. ! =T print eddy coef. in numout 66 !! 67 NAMELIST/namtra_ldf/ ln_traldf_lap , ln_traldf_bilap, & 68 & ln_traldf_level, ln_traldf_hor , ln_traldf_iso, & 69 & ln_traldf_grif , ln_traldf_gdia , & 70 & ln_triad_iso , ln_botmix_grif , & 71 & rn_aht_0 , rn_ahtb_0 , rn_aeiv_0, & 72 & rn_slpmax , rn_chsmag , rn_smsh, & 73 & rn_aht_m 74 !!---------------------------------------------------------------------- 75 76 ! Define the lateral tracer physics parameters 77 ! ============================================= 78 79 97 !! ** Method : * the eddy diffusivity coef. specification depends on: 98 !! 99 !! ln_traldf_lap = T laplacian operator 100 !! ln_traldf_blp = T bilaplacian operator 101 !! 102 !! nn_aht_ijk_t = 0 => = constant 103 !! ! 104 !! = 10 => = F(z) : constant with a reduction of 1/4 with depth 105 !! ! 106 !! =-20 => = F(i,j) = shape read in 'eddy_diffusivity.nc' file 107 !! = 20 = F(i,j) = F(e1,e2) or F(e1^3,e2^3) (lap or bilap case) 108 !! = 21 = F(i,j,t) = F(growth rate of baroclinic instability) 109 !! ! 110 !! =-30 => = F(i,j,k) = shape read in 'eddy_diffusivity.nc' file 111 !! = 30 = F(i,j,k) = 2D (case 20) + decrease with depth (case 10) 112 !! = 31 = F(i,j,k,t) = F(local velocity) ( |u|e /12 laplacian operator 113 !! or |u|e^3/12 bilaplacian operator ) 114 !! * initialisation of the eddy induced velocity coefficient by a call to ldf_eiv_init 115 !! 116 !! ** action : ahtu, ahtv initialized once for all or l_ldftra_time set to true 117 !! aeiu, aeiv initialized once for all or l_ldfeiv_time set to true 118 !!---------------------------------------------------------------------- 119 INTEGER :: jk ! dummy loop indices 120 INTEGER :: ierr, inum, ios ! local integer 121 REAL(wp) :: zah0 ! local scalar 122 ! 123 NAMELIST/namtra_ldf/ ln_traldf_lap, ln_traldf_blp , & ! type of operator 124 & ln_traldf_lev, ln_traldf_hor , ln_traldf_triad, & ! acting direction of the operator 125 & ln_traldf_iso, ln_traldf_msc , rn_slpmax , & ! option for iso-neutral operator 126 & ln_triad_iso , ln_botmix_triad, rn_sw_triad , & ! option for triad operator 127 & rn_aht_0 , rn_bht_0 , nn_aht_ijk_t ! lateral eddy coefficient 128 !!---------------------------------------------------------------------- 129 ! 130 ! Choice of lateral tracer physics 131 ! ================================= 132 ! 80 133 REWIND( numnam_ref ) ! Namelist namtra_ldf in reference namelist : Lateral physics on tracers 81 134 READ ( numnam_ref, namtra_ldf, IOSTAT = ios, ERR = 901) 82 135 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldf in reference namelist', lwp ) 83 136 ! 84 137 REWIND( numnam_cfg ) ! Namelist namtra_ldf in configuration namelist : Lateral physics on tracers 85 138 READ ( numnam_cfg, namtra_ldf, IOSTAT = ios, ERR = 902 ) 86 139 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldf in configuration namelist', lwp ) 87 140 IF(lwm) WRITE ( numond, namtra_ldf ) 88 141 ! 89 142 IF(lwp) THEN ! control print 90 143 WRITE(numout,*) … … 92 145 WRITE(numout,*) '~~~~~~~~~~~~ ' 93 146 WRITE(numout,*) ' Namelist namtra_ldf : lateral mixing parameters (type, direction, coefficients)' 94 WRITE(numout,*) ' laplacian operator ln_traldf_lap = ', ln_traldf_lap 95 WRITE(numout,*) ' bilaplacian operator ln_traldf_bilap = ', ln_traldf_bilap 96 WRITE(numout,*) ' iso-level ln_traldf_level = ', ln_traldf_level 97 WRITE(numout,*) ' horizontal (geopotential) ln_traldf_hor = ', ln_traldf_hor 98 WRITE(numout,*) ' iso-neutral ln_traldf_iso = ', ln_traldf_iso 99 WRITE(numout,*) ' iso-neutral (Griffies) ln_traldf_grif = ', ln_traldf_grif 100 WRITE(numout,*) ' Griffies strmfn diagnostics ln_traldf_gdia = ', ln_traldf_gdia 101 WRITE(numout,*) ' lateral eddy diffusivity rn_aht_0 = ', rn_aht_0 102 WRITE(numout,*) ' background hor. diffusivity rn_ahtb_0 = ', rn_ahtb_0 103 WRITE(numout,*) ' eddy induced velocity coef. rn_aeiv_0 = ', rn_aeiv_0 104 WRITE(numout,*) ' maximum isoppycnal slope rn_slpmax = ', rn_slpmax 105 WRITE(numout,*) ' pure lateral mixing in ML ln_triad_iso = ', ln_triad_iso 106 WRITE(numout,*) ' lateral mixing on bottom ln_botmix_grif = ', ln_botmix_grif 147 ! 148 WRITE(numout,*) ' type :' 149 WRITE(numout,*) ' laplacian operator ln_traldf_lap = ', ln_traldf_lap 150 WRITE(numout,*) ' bilaplacian operator ln_traldf_blp = ', ln_traldf_blp 151 ! 152 WRITE(numout,*) ' direction of action :' 153 WRITE(numout,*) ' iso-level ln_traldf_lev = ', ln_traldf_lev 154 WRITE(numout,*) ' horizontal (geopotential) ln_traldf_hor = ', ln_traldf_hor 155 WRITE(numout,*) ' iso-neutral Madec operator ln_traldf_iso = ', ln_traldf_iso 156 WRITE(numout,*) ' iso-neutral triad operator ln_traldf_triad = ', ln_traldf_triad 157 WRITE(numout,*) ' iso-neutral (Method of Stab. Corr.) ln_traldf_msc = ', ln_traldf_msc 158 WRITE(numout,*) ' maximum isoppycnal slope rn_slpmax = ', rn_slpmax 159 WRITE(numout,*) ' pure lateral mixing in ML ln_triad_iso = ', ln_triad_iso 160 WRITE(numout,*) ' switching triad or not rn_sw_triad = ', rn_sw_triad 161 WRITE(numout,*) ' lateral mixing on bottom ln_botmix_triad = ', ln_botmix_triad 162 ! 163 WRITE(numout,*) ' coefficients :' 164 WRITE(numout,*) ' lateral eddy diffusivity (lap case) rn_aht_0 = ', rn_aht_0 165 WRITE(numout,*) ' lateral eddy diffusivity (bilap case) rn_bht_0 = ', rn_bht_0 166 WRITE(numout,*) ' type of time-space variation nn_aht_ijk_t = ', nn_aht_ijk_t 167 ENDIF 168 ! 169 ! ! Parameter control 170 ! 171 IF( .NOT.ln_traldf_lap .AND. .NOT.ln_traldf_blp ) THEN 172 IF(lwp) WRITE(numout,*) ' No diffusive operator selected. ahtu and ahtv are not allocated' 173 l_ldftra_time = .FALSE. 174 RETURN 175 ENDIF 176 ! 177 IF( ln_traldf_blp .AND. ( ln_traldf_iso .OR. ln_traldf_triad) ) THEN ! iso-neutral bilaplacian need MSC 178 IF( .NOT.ln_traldf_msc ) CALL ctl_stop( 'tra_ldf_init: iso-neutral bilaplacian requires ln_traldf_msc=.true.' ) 179 ENDIF 180 ! 181 ! Space/time variation of eddy coefficients 182 ! =========================================== 183 ! ! allocate the aht arrays 184 ALLOCATE( ahtu(jpi,jpj,jpk) , ahtv(jpi,jpj,jpk) , STAT=ierr ) 185 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_tra_init: failed to allocate arrays') 186 ! 187 ahtu(:,:,jpk) = 0._wp ! last level always 0 188 ahtv(:,:,jpk) = 0._wp 189 ! 190 ! ! value of eddy mixing coef. 191 IF ( ln_traldf_lap ) THEN ; zah0 = rn_aht_0 ! laplacian operator 192 ELSEIF( ln_traldf_blp ) THEN ; zah0 = ABS( rn_bht_0 ) ! bilaplacian operator 193 ENDIF 194 ! 195 l_ldftra_time = .FALSE. ! no time variation except in case defined below 196 ! 197 IF( ln_traldf_lap .OR. ln_traldf_blp ) THEN ! only if a lateral diffusion operator is used 198 ! 199 SELECT CASE( nn_aht_ijk_t ) ! Specification of space time variations of ehtu, ahtv 200 ! 201 CASE( 0 ) !== constant ==! 202 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = constant = ', rn_aht_0 203 ahtu(:,:,:) = zah0 * umask(:,:,:) 204 ahtv(:,:,:) = zah0 * vmask(:,:,:) 205 ! 206 CASE( 10 ) !== fixed profile ==! 207 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F( depth )' 208 ahtu(:,:,1) = zah0 * umask(:,:,1) ! constant surface value 209 ahtv(:,:,1) = zah0 * vmask(:,:,1) 210 CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) 211 ! 212 CASE ( -20 ) !== fixed horizontal shape read in file ==! 213 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F(i,j) read in eddy_diffusivity.nc file' 214 CALL iom_open( 'eddy_diffusivity_2D.nc', inum ) 215 CALL iom_get ( inum, jpdom_data, 'ahtu_2D', ahtu(:,:,1) ) 216 CALL iom_get ( inum, jpdom_data, 'ahtv_2D', ahtv(:,:,1) ) 217 CALL iom_close( inum ) 218 DO jk = 2, jpkm1 219 ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk) 220 ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk) 221 END DO 222 ! 223 CASE( 20 ) !== fixed horizontal shape ==! 224 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or blp case)' 225 IF( ln_traldf_lap ) CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor 226 IF( ln_traldf_blp ) CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor 227 ! 228 CASE( 21 ) !== time varying 2D field ==! 229 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F( latitude, longitude, time )' 230 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 231 IF(lwp) WRITE(numout,*) ' min value = 0.1 * rn_aht_0' 232 IF(lwp) WRITE(numout,*) ' max value = rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21)' 233 IF(lwp) WRITE(numout,*) ' increased to rn_aht_0 within 20N-20S' 234 ! 235 l_ldftra_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 236 ! 237 IF( ln_traldf_blp ) THEN 238 CALL ctl_stop( 'ldf_tra_init: aht=F(growth rate of baroc. insta.) incompatible with bilaplacian operator' ) 239 ENDIF 240 ! 241 CASE( -30 ) !== fixed 3D shape read in file ==! 242 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F(i,j,k) read in eddy_diffusivity.nc file' 243 CALL iom_open( 'eddy_diffusivity_3D.nc', inum ) 244 CALL iom_get ( inum, jpdom_data, 'ahtu_3D', ahtu ) 245 CALL iom_get ( inum, jpdom_data, 'ahtv_3D', ahtv ) 246 CALL iom_close( inum ) 247 DO jk = 1, jpkm1 248 ahtu(:,:,jk) = ahtu(:,:,jk) * umask(:,:,jk) 249 ahtv(:,:,jk) = ahtv(:,:,jk) * vmask(:,:,jk) 250 END DO 251 ! 252 CASE( 30 ) !== fixed 3D shape ==! 253 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F( latitude, longitude, depth )' 254 IF( ln_traldf_lap ) CALL ldf_c2d( 'TRA', 'LAP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor 255 IF( ln_traldf_blp ) CALL ldf_c2d( 'TRA', 'BLP', zah0, ahtu, ahtv ) ! surface value proportional to scale factor 256 ! ! reduction with depth 257 CALL ldf_c1d( 'TRA', r1_4, ahtu(:,:,1), ahtv(:,:,1), ahtu, ahtv ) 258 ! 259 CASE( 31 ) !== time varying 3D field ==! 260 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F( latitude, longitude, depth , time )' 261 IF(lwp) WRITE(numout,*) ' proportional to the velocity : |u|e/12 or |u|e^3/12' 262 ! 263 l_ldftra_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 264 ! 265 CASE DEFAULT 266 CALL ctl_stop('ldf_tra_init: wrong choice for nn_aht_ijk_t, the type of space-time variation of aht') 267 END SELECT 268 ! 269 IF( ln_traldf_blp .AND. .NOT. l_ldftra_time ) THEN 270 ahtu(:,:,:) = SQRT( ahtu(:,:,:) ) 271 ahtv(:,:,:) = SQRT( ahtv(:,:,:) ) 272 ENDIF 273 ! 274 ENDIF 275 ! 276 END SUBROUTINE ldf_tra_init 277 278 279 SUBROUTINE ldf_tra( kt ) 280 !!---------------------------------------------------------------------- 281 !! *** ROUTINE ldf_tra *** 282 !! 283 !! ** Purpose : update at kt the tracer lateral mixing coeff. (aht and aeiv) 284 !! 285 !! ** Method : time varying eddy diffusivity coefficients: 286 !! 287 !! nn_aei_ijk_t = 21 aeiu, aeiv = F(i,j, t) = F(growth rate of baroclinic instability) 288 !! with a reduction to 0 in vicinity of the Equator 289 !! nn_aht_ijk_t = 21 ahtu, ahtv = F(i,j, t) = F(growth rate of baroclinic instability) 290 !! 291 !! = 31 ahtu, ahtv = F(i,j,k,t) = F(local velocity) ( |u|e /12 laplacian operator 292 !! or |u|e^3/12 bilaplacian operator ) 293 !! 294 !! ** action : ahtu, ahtv update at each time step 295 !! aeiu, aeiv - - - - (if ln_ldfeiv=T) 296 !!---------------------------------------------------------------------- 297 INTEGER, INTENT(in) :: kt ! time step 298 ! 299 INTEGER :: ji, jj, jk ! dummy loop indices 300 REAL(wp) :: zaht, zaht_min, z1_f20 ! local scalar 301 !!---------------------------------------------------------------------- 302 ! 303 IF( nn_aei_ijk_t == 21 ) THEN ! eddy induced velocity coefficients 304 ! ! =F(growth rate of baroclinic instability) 305 ! ! max value rn_aeiv_0 ; decreased to 0 within 20N-20S 306 CALL ldf_eiv( kt, rn_aeiv_0, aeiu, aeiv ) 307 IF(lwp .AND. kt<=nit000+20 ) WRITE(numout,*) ' kt , ldf_eiv appel', kt 308 ENDIF 309 ! 310 SELECT CASE( nn_aht_ijk_t ) ! Eddy diffusivity coefficients 311 ! 312 CASE( 21 ) !== time varying 2D field ==! = F( growth rate of baroclinic instability ) 313 ! ! min value rn_aht_0 / 10 314 ! ! max value rn_aht_0 (rn_aeiv_0 if nn_aei_ijk_t=21) 315 ! ! increase to rn_aht_0 within 20N-20S 316 IF( nn_aei_ijk_t /= 21 ) THEN 317 CALL ldf_eiv( kt, rn_aht_0, ahtu, ahtv ) 318 IF(lwp .AND. kt<=nit000+20 ) WRITE(numout,*) ' kt , ldf_eiv appel 2', kt 319 ELSE 320 ahtu(:,:,1) = aeiu(:,:,1) 321 ahtv(:,:,1) = aeiv(:,:,1) 322 IF(lwp .AND. kt<=nit000+20 ) WRITE(numout,*) ' kt , ahtu=aeiu', kt 323 ENDIF 324 ! 325 z1_f20 = 1._wp / ( 2._wp * omega * SIN( rad * 20._wp ) ) ! 1 / ff(20 degrees) 326 zaht_min = 0.2_wp * rn_aht_0 ! minimum value for aht 327 DO jj = 1, jpj 328 DO ji = 1, jpi 329 zaht = ( 1._wp - MIN( 1._wp , ABS( ff(ji,jj) * z1_f20 ) ) ) * ( rn_aht_0 - zaht_min ) 330 ahtu(ji,jj,1) = ( MAX( zaht_min, ahtu(ji,jj,1) ) + zaht ) * umask(ji,jj,1) ! min value zaht_min 331 ahtv(ji,jj,1) = ( MAX( zaht_min, ahtv(ji,jj,1) ) + zaht ) * vmask(ji,jj,1) ! increase within 20S-20N 332 END DO 333 END DO 334 DO jk = 2, jpkm1 ! deeper value = surface value 335 ahtu(:,:,jk) = ahtu(:,:,1) * umask(:,:,jk) 336 ahtv(:,:,jk) = ahtv(:,:,1) * vmask(:,:,jk) 337 END DO 338 ! 339 CASE( 31 ) !== time varying 3D field ==! = F( local velocity ) 340 IF( ln_traldf_lap ) THEN ! laplacian operator |u| e /12 341 DO jk = 1, jpkm1 342 ahtu(:,:,jk) = ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 343 ahtv(:,:,jk) = ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12 344 END DO 345 ELSEIF( ln_traldf_blp ) THEN ! bilaplacian operator sqrt( |u| e^3 /12 ) = sqrt( |u| e /12 ) * e 346 DO jk = 1, jpkm1 347 ahtu(:,:,jk) = SQRT( ABS( ub(:,:,jk) ) * e1u(:,:) * r1_12 ) * e1u(:,:) 348 ahtv(:,:,jk) = SQRT( ABS( vb(:,:,jk) ) * e2v(:,:) * r1_12 ) * e2v(:,:) 349 END DO 350 ENDIF 351 ! 352 END SELECT 353 ! 354 IF( .NOT.lk_offline ) THEN 355 CALL iom_put( "ahtu_2d", ahtu(:,:,1) ) ! surface u-eddy diffusivity coeff. 356 CALL iom_put( "ahtv_2d", ahtv(:,:,1) ) ! surface v-eddy diffusivity coeff. 357 CALL iom_put( "ahtu_3d", ahtu(:,:,:) ) ! 3D u-eddy diffusivity coeff. 358 CALL iom_put( "ahtv_3d", ahtv(:,:,:) ) ! 3D v-eddy diffusivity coeff. 359 ! 360 !!gm : THE IF below is to be checked (comes from Seb) 361 IF( ln_ldfeiv ) THEN 362 CALL iom_put( "aeiu_2d", aeiu(:,:,1) ) ! surface u-EIV coeff. 363 CALL iom_put( "aeiv_2d", aeiv(:,:,1) ) ! surface v-EIV coeff. 364 CALL iom_put( "aeiu_3d", aeiu(:,:,:) ) ! 3D u-EIV coeff. 365 CALL iom_put( "aeiv_3d", aeiv(:,:,:) ) ! 3D v-EIV coeff. 366 ENDIF 367 ENDIF 368 ! 369 END SUBROUTINE ldf_tra 370 371 372 SUBROUTINE ldf_eiv_init 373 !!---------------------------------------------------------------------- 374 !! *** ROUTINE ldf_eiv_init *** 375 !! 376 !! ** Purpose : initialization of the eiv coeff. from namelist choices. 377 !! 378 !! ** Method : 379 !! 380 !! ** Action : aeiu , aeiv : EIV coeff. at u- & v-points 381 !! l_ldfeiv_time : =T if EIV coefficients vary with time 382 !!---------------------------------------------------------------------- 383 INTEGER :: jk ! dummy loop indices 384 INTEGER :: ierr, inum, ios ! local integer 385 ! 386 NAMELIST/namtra_ldfeiv/ ln_ldfeiv , ln_ldfeiv_dia, & ! eddy induced velocity (eiv) 387 & nn_aei_ijk_t, rn_aeiv_0 ! eiv coefficient 388 !!---------------------------------------------------------------------- 389 ! 390 REWIND( numnam_ref ) ! Namelist namtra_ldfeiv in reference namelist : eddy induced velocity param. 391 READ ( numnam_ref, namtra_ldfeiv, IOSTAT = ios, ERR = 901) 392 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldfeiv in reference namelist', lwp ) 393 ! 394 REWIND( numnam_cfg ) ! Namelist namtra_ldfeiv in configuration namelist : eddy induced velocity param. 395 READ ( numnam_cfg, namtra_ldfeiv, IOSTAT = ios, ERR = 902 ) 396 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_ldfeiv in configuration namelist', lwp ) 397 IF(lwm) WRITE ( numond, namtra_ldfeiv ) 398 399 IF(lwp) THEN ! control print 107 400 WRITE(numout,*) 108 ENDIF 109 110 ! ! convert DOCTOR namelist names into OLD names 111 aht0 = rn_aht_0 112 ahtb0 = rn_ahtb_0 113 aeiv0 = rn_aeiv_0 114 115 ! ! Parameter control 116 117 ! ... Check consistency for type and direction : 118 ! ==> will be done in traldf module 119 120 ! ... Space variation of eddy coefficients 121 ioptio = 0 122 #if defined key_traldf_c3d 123 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F( latitude, longitude, depth)' 124 ioptio = ioptio + 1 125 #endif 126 #if defined key_traldf_c2d 127 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F( latitude, longitude)' 128 ioptio = ioptio + 1 129 #endif 130 #if defined key_traldf_c1d 131 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F( depth )' 132 ioptio = ioptio + 1 133 IF( .NOT. ln_zco ) CALL ctl_stop( 'key_traldf_c1d can only be used in z-coordinate - full step' ) 134 #endif 135 IF( ioptio == 0 ) THEN 136 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = constant (default option)' 137 ELSEIF( ioptio > 1 ) THEN 138 CALL ctl_stop(' use only one of the following keys:', & 139 & ' key_traldf_c3d, key_traldf_c2d, key_traldf_c1d' ) 140 ENDIF 141 142 IF( ln_traldf_bilap ) THEN 143 IF(lwp) WRITE(numout,*) ' biharmonic tracer diffusion' 144 IF( aht0 > 0 .AND. .NOT. lk_esopa ) CALL ctl_stop( 'The horizontal diffusivity coef. aht0 must be negative' ) 401 WRITE(numout,*) 'ldf_eiv_init : eddy induced velocity parametrization' 402 WRITE(numout,*) '~~~~~~~~~~~~ ' 403 WRITE(numout,*) ' Namelist namtra_ldfeiv : ' 404 WRITE(numout,*) ' Eddy Induced Velocity (eiv) param. ln_ldfeiv = ', ln_ldfeiv 405 WRITE(numout,*) ' eiv streamfunction & velocity diag. ln_ldfeiv_dia = ', ln_ldfeiv_dia 406 WRITE(numout,*) ' eddy induced velocity coef. rn_aeiv_0 = ', rn_aeiv_0 407 WRITE(numout,*) ' type of time-space variation nn_aei_ijk_t = ', nn_aei_ijk_t 408 WRITE(numout,*) 409 ENDIF 410 ! 411 IF( ln_ldfeiv .AND. ln_traldf_blp ) CALL ctl_stop( 'ldf_eiv_init: eddy induced velocity ONLY with laplacian diffusivity' ) 412 413 ! ! Parameter control 414 l_ldfeiv_time = .FALSE. 415 ! 416 IF( ln_ldfeiv ) THEN ! allocate the aei arrays 417 ALLOCATE( aeiu(jpi,jpj,jpk), aeiv(jpi,jpj,jpk), STAT=ierr ) 418 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'ldf_eiv: failed to allocate arrays') 419 ! 420 SELECT CASE( nn_aei_ijk_t ) ! Specification of space time variations of eaiu, aeiv 421 ! 422 CASE( 0 ) !== constant ==! 423 IF(lwp) WRITE(numout,*) ' eddy induced velocity coef. = constant = ', rn_aeiv_0 424 aeiu(:,:,:) = rn_aeiv_0 425 aeiv(:,:,:) = rn_aeiv_0 426 ! 427 CASE( 10 ) !== fixed profile ==! 428 IF(lwp) WRITE(numout,*) ' eddy induced velocity coef. = F( depth )' 429 aeiu(:,:,1) = rn_aeiv_0 ! constant surface value 430 aeiv(:,:,1) = rn_aeiv_0 431 CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 432 ! 433 CASE ( -20 ) !== fixed horizontal shape read in file ==! 434 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F(i,j) read in eddy_diffusivity_2D.nc file' 435 CALL iom_open ( 'eddy_induced_velocity_2D.nc', inum ) 436 CALL iom_get ( inum, jpdom_data, 'aeiu', aeiu(:,:,1) ) 437 CALL iom_get ( inum, jpdom_data, 'aeiv', aeiv(:,:,1) ) 438 CALL iom_close( inum ) 439 DO jk = 2, jpk 440 aeiu(:,:,jk) = aeiu(:,:,1) 441 aeiv(:,:,jk) = aeiv(:,:,1) 442 END DO 443 ! 444 CASE( 20 ) !== fixed horizontal shape ==! 445 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap or bilap case)' 446 CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv ) ! surface value proportional to scale factor 447 ! 448 CASE( 21 ) !== time varying 2D field ==! 449 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F( latitude, longitude, time )' 450 IF(lwp) WRITE(numout,*) ' = F( growth rate of baroclinic instability )' 451 ! 452 l_ldfeiv_time = .TRUE. ! will be calculated by call to ldf_tra routine in step.F90 453 ! 454 CASE( -30 ) !== fixed 3D shape read in file ==! 455 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file' 456 CALL iom_open ( 'eddy_induced_velocity_3D.nc', inum ) 457 CALL iom_get ( inum, jpdom_data, 'aeiu', aeiu ) 458 CALL iom_get ( inum, jpdom_data, 'aeiv', aeiv ) 459 CALL iom_close( inum ) 460 ! 461 CASE( 30 ) !== fixed 3D shape ==! 462 IF(lwp) WRITE(numout,*) ' tracer mixing coef. = F( latitude, longitude, depth )' 463 CALL ldf_c2d( 'TRA', 'LAP', rn_aeiv_0, aeiu, aeiv ) ! surface value proportional to scale factor 464 ! ! reduction with depth 465 CALL ldf_c1d( 'TRA', r1_4, aeiu(:,:,1), aeiv(:,:,1), aeiu, aeiv ) 466 ! 467 CASE DEFAULT 468 CALL ctl_stop('ldf_tra_init: wrong choice for nn_aei_ijk_t, the type of space-time variation of aei') 469 END SELECT 470 ! 145 471 ELSE 146 IF(lwp) WRITE(numout,*) ' harmonic tracer diffusion (default)' 147 IF( aht0 < 0 .AND. .NOT. lk_esopa ) CALL ctl_stop('The horizontal diffusivity coef. aht0 must be positive' ) 148 ENDIF 149 150 151 ! Lateral eddy diffusivity and eddy induced velocity coefficients 152 ! ================================================================ 153 #if defined key_traldf_c3d 154 CALL ldf_tra_c3d( ll_print ) ! aht = 3D coef. = F( longitude, latitude, depth ) 155 #elif defined key_traldf_c2d 156 CALL ldf_tra_c2d( ll_print ) ! aht = 2D coef. = F( longitude, latitude ) 157 #elif defined key_traldf_c1d 158 CALL ldf_tra_c1d( ll_print ) ! aht = 1D coef. = F( depth ) 159 #else 160 ! Constant coefficients 161 IF(lwp)WRITE(numout,*) 162 IF(lwp)WRITE(numout,*) ' constant eddy diffusivity coef. ahtu = ahtv = ahtw = aht0 = ', aht0 163 IF( lk_traldf_eiv ) THEN 164 IF(lwp)WRITE(numout,*) ' constant eddy induced velocity coef. aeiu = aeiv = aeiw = aeiv0 = ', aeiv0 472 IF(lwp) WRITE(numout,*) ' eddy induced velocity param is NOT used neither diagnosed' 473 ln_ldfeiv_dia = .FALSE. 474 ENDIF 475 ! 476 END SUBROUTINE ldf_eiv_init 477 478 479 SUBROUTINE ldf_eiv( kt, paei0, paeiu, paeiv ) 480 !!---------------------------------------------------------------------- 481 !! *** ROUTINE ldf_eiv *** 482 !! 483 !! ** Purpose : Compute the eddy induced velocity coefficient from the 484 !! growth rate of baroclinic instability. 485 !! 486 !! ** Method : coefficient function of the growth rate of baroclinic instability 487 !! 488 !! Reference : Treguier et al. JPO 1997 ; Held and Larichev JAS 1996 489 !!---------------------------------------------------------------------- 490 INTEGER , INTENT(in ) :: kt ! ocean time-step index 491 REAL(wp) , INTENT(inout) :: paei0 ! max value [m2/s] 492 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: paeiu, paeiv ! eiv coefficient [m2/s] 493 ! 494 INTEGER :: ji, jj, jk ! dummy loop indices 495 REAL(wp) :: zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei ! local scalars 496 REAL(wp), DIMENSION(:,:), POINTER :: zn, zah, zhw, zross, zaeiw ! 2D workspace 497 !!---------------------------------------------------------------------- 498 ! 499 IF( nn_timing == 1 ) CALL timing_start('ldf_eiv') 500 ! 501 CALL wrk_alloc( jpi,jpj, zn, zah, zhw, zross, zaeiw ) 502 ! 503 zn (:,:) = 0._wp ! Local initialization 504 zhw (:,:) = 5._wp 505 zah (:,:) = 0._wp 506 zross(:,:) = 0._wp 507 ! ! Compute lateral diffusive coefficient at T-point 508 IF( ln_traldf_triad ) THEN 509 DO jk = 1, jpk 510 DO jj = 2, jpjm1 511 DO ji = 2, jpim1 512 ! Take the max of N^2 and zero then take the vertical sum 513 ! of the square root of the resulting N^2 ( required to compute 514 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 515 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 516 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk) 517 ! Compute elements required for the inverse time scale of baroclinic 518 ! eddies using the isopycnal slopes calculated in ldfslp.F : 519 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 520 ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk) 521 zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w 522 zhw(ji,jj) = zhw(ji,jj) + ze3w 523 END DO 524 END DO 525 END DO 526 ELSE 527 DO jk = 1, jpk 528 DO jj = 2, jpjm1 529 DO ji = 2, jpim1 530 ! Take the max of N^2 and zero then take the vertical sum 531 ! of the square root of the resulting N^2 ( required to compute 532 ! internal Rossby radius Ro = .5 * sum_jpk(N) / f 533 zn2 = MAX( rn2b(ji,jj,jk), 0._wp ) 534 zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * e3w_n(ji,jj,jk) 535 ! Compute elements required for the inverse time scale of baroclinic 536 ! eddies using the isopycnal slopes calculated in ldfslp.F : 537 ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 538 ze3w = e3w_n(ji,jj,jk) * tmask(ji,jj,jk) 539 zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 540 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w 541 zhw(ji,jj) = zhw(ji,jj) + ze3w 542 END DO 543 END DO 544 END DO 545 END IF 546 547 DO jj = 2, jpjm1 548 DO ji = fs_2, fs_jpim1 ! vector opt. 549 zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 550 ! Rossby radius at w-point taken < 40km and > 2km 551 zross(ji,jj) = MAX( MIN( .4 * zn(ji,jj) / zfw, 40.e3 ), 2.e3 ) 552 ! Compute aeiw by multiplying Ro^2 and T^-1 553 zaeiw(ji,jj) = zross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 554 END DO 555 END DO 556 557 !!gm IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 558 !!gm DO jj = 2, jpjm1 559 !!gm DO ji = fs_2, fs_jpim1 ! vector opt. 560 !!gm ! Take the minimum between aeiw and 1000 m2/s over shelves (depth shallower than 650 m) 561 !!gm IF( mbkt(ji,jj) <= 20 ) zaeiw(ji,jj) = MIN( zaeiw(ji,jj), 1000. ) 562 !!gm END DO 563 !!gm END DO 564 !!gm ENDIF 565 566 ! !== Bound on eiv coeff. ==! 567 z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) ) 568 DO jj = 2, jpjm1 569 DO ji = fs_2, fs_jpim1 ! vector opt. 570 zzaei = MIN( 1._wp, ABS( ff(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj) ! tropical decrease 571 zaeiw(ji,jj) = MIN( zzaei , paei0 ) ! Max value = paei0 572 END DO 573 END DO 574 CALL lbc_lnk( zaeiw(:,:), 'W', 1. ) ! lateral boundary condition 575 ! 576 DO jj = 2, jpjm1 !== aei at u- and v-points ==! 577 DO ji = fs_2, fs_jpim1 ! vector opt. 578 paeiu(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji+1,jj ) ) * umask(ji,jj,1) 579 paeiv(ji,jj,1) = 0.5_wp * ( zaeiw(ji,jj) + zaeiw(ji ,jj+1) ) * vmask(ji,jj,1) 580 END DO 581 END DO 582 CALL lbc_lnk( paeiu(:,:,1), 'U', 1. ) ; CALL lbc_lnk( paeiv(:,:,1), 'V', 1. ) ! lateral boundary condition 583 584 DO jk = 2, jpkm1 !== deeper values equal the surface one ==! 585 paeiu(:,:,jk) = paeiu(:,:,1) * umask(:,:,jk) 586 paeiv(:,:,jk) = paeiv(:,:,1) * vmask(:,:,jk) 587 END DO 588 ! 589 CALL wrk_dealloc( jpi,jpj, zn, zah, zhw, zross, zaeiw ) 590 ! 591 IF( nn_timing == 1 ) CALL timing_stop('ldf_eiv') 592 ! 593 END SUBROUTINE ldf_eiv 594 595 596 SUBROUTINE ldf_eiv_trp( kt, kit000, pun, pvn, pwn, cdtype ) 597 !!---------------------------------------------------------------------- 598 !! *** ROUTINE ldf_eiv_trp *** 599 !! 600 !! ** Purpose : add to the input ocean transport the contribution of 601 !! the eddy induced velocity parametrization. 602 !! 603 !! ** Method : The eddy induced transport is computed from a flux stream- 604 !! function which depends on the slope of iso-neutral surfaces 605 !! (see ldf_slp). For example, in the i-k plan : 606 !! psi_uw = mk(aeiu) e2u mi(wslpi) [in m3/s] 607 !! Utr_eiv = - dk[psi_uw] 608 !! Vtr_eiv = + di[psi_uw] 609 !! ln_ldfeiv_dia = T : output the associated streamfunction, 610 !! velocity and heat transport (call ldf_eiv_dia) 611 !! 612 !! ** Action : pun, pvn increased by the eiv transport 613 !!---------------------------------------------------------------------- 614 INTEGER , INTENT(in ) :: kt ! ocean time-step index 615 INTEGER , INTENT(in ) :: kit000 ! first time step index 616 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 617 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pun ! in : 3 ocean transport components [m3/s] 618 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pvn ! out: 3 ocean transport components [m3/s] 619 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pwn ! increased by the eiv [m3/s] 620 !! 621 INTEGER :: ji, jj, jk ! dummy loop indices 622 REAL(wp) :: zuwk, zuwk1, zuwi, zuwi1 ! local scalars 623 REAL(wp) :: zvwk, zvwk1, zvwj, zvwj1 ! - - 624 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpsi_uw, zpsi_vw 625 !!---------------------------------------------------------------------- 626 ! 627 IF( nn_timing == 1 ) CALL timing_start( 'ldf_eiv_trp') 628 ! 629 CALL wrk_alloc( jpi,jpj,jpk, zpsi_uw, zpsi_vw ) 630 631 IF( kt == kit000 ) THEN 632 IF(lwp) WRITE(numout,*) 633 IF(lwp) WRITE(numout,*) 'ldf_eiv_trp : eddy induced advection on ', cdtype,' :' 634 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ add to velocity fields the eiv component' 635 ENDIF 636 165 637 166 ENDIF 167 #endif 168 169 #if defined key_traldf_smag && ! defined key_traldf_c3d 170 CALL ctl_stop( 'key_traldf_smag can only be used with key_traldf_c3d' ) 171 #endif 172 #if defined key_traldf_smag 173 IF(lwp) WRITE(numout,*)' SMAGORINSKY DIFFUSION' 174 IF(lwp .AND. rn_smsh < 1) WRITE(numout,*)' only shear is used ' 175 IF(lwp.and.ln_traldf_bilap) CALL ctl_stop(' SMAGORINSKY + BILAPLACIAN - UNSTABLE OR NON_CONSERVATIVE' ) 176 #endif 177 178 ! 179 END SUBROUTINE ldf_tra_init 180 181 #if defined key_traldf_c3d 182 # include "ldftra_c3d.h90" 183 #elif defined key_traldf_c2d 184 # include "ldftra_c2d.h90" 185 #elif defined key_traldf_c1d 186 # include "ldftra_c1d.h90" 187 #endif 638 zpsi_uw(:,:, 1 ) = 0._wp ; zpsi_vw(:,:, 1 ) = 0._wp 639 zpsi_uw(:,:,jpk) = 0._wp ; zpsi_vw(:,:,jpk) = 0._wp 640 ! 641 DO jk = 2, jpkm1 642 DO jj = 1, jpjm1 643 DO ji = 1, fs_jpim1 ! vector opt. 644 zpsi_uw(ji,jj,jk) = - 0.25_wp * e2u(ji,jj) * ( wslpi(ji,jj,jk ) + wslpi(ji+1,jj,jk) ) & 645 & * ( aeiu (ji,jj,jk-1) + aeiu (ji ,jj,jk) ) * umask(ji,jj,jk) 646 zpsi_vw(ji,jj,jk) = - 0.25_wp * e1v(ji,jj) * ( wslpj(ji,jj,jk ) + wslpj(ji,jj+1,jk) ) & 647 & * ( aeiv (ji,jj,jk-1) + aeiv (ji,jj ,jk) ) * vmask(ji,jj,jk) 648 END DO 649 END DO 650 END DO 651 ! 652 DO jk = 1, jpkm1 653 DO jj = 1, jpjm1 654 DO ji = 1, fs_jpim1 ! vector opt. 655 pun(ji,jj,jk) = pun(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 656 pvn(ji,jj,jk) = pvn(ji,jj,jk) - ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 657 END DO 658 END DO 659 END DO 660 DO jk = 1, jpkm1 661 DO jj = 2, jpjm1 662 DO ji = fs_2, fs_jpim1 ! vector opt. 663 pwn(ji,jj,jk) = pwn(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj ,jk) & 664 & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji ,jj-1,jk) ) 665 END DO 666 END DO 667 END DO 668 ! 669 ! ! diagnose the eddy induced velocity and associated heat transport 670 IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) CALL ldf_eiv_dia( zpsi_uw, zpsi_vw ) 671 ! 672 CALL wrk_dealloc( jpi,jpj,jpk, zpsi_uw, zpsi_vw ) 673 ! 674 IF( nn_timing == 1 ) CALL timing_stop( 'ldf_eiv_trp') 675 ! 676 END SUBROUTINE ldf_eiv_trp 677 678 679 SUBROUTINE ldf_eiv_dia( psi_uw, psi_vw ) 680 !!---------------------------------------------------------------------- 681 !! *** ROUTINE ldf_eiv_dia *** 682 !! 683 !! ** Purpose : diagnose the eddy induced velocity and its associated 684 !! vertically integrated heat transport. 685 !! 686 !! ** Method : 687 !! 688 !!---------------------------------------------------------------------- 689 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: psi_uw, psi_vw ! streamfunction [m3/s] 690 ! 691 INTEGER :: ji, jj, jk ! dummy loop indices 692 REAL(wp) :: zztmp ! local scalar 693 REAL(wp), DIMENSION(:,:) , POINTER :: zw2d ! 2D workspace 694 REAL(wp), DIMENSION(:,:,:), POINTER :: zw3d ! 3D workspace 695 !!---------------------------------------------------------------------- 696 ! 697 IF( nn_timing == 1 ) CALL timing_start( 'ldf_eiv_dia') 698 ! 699 ! !== eiv stream function: output ==! 700 CALL lbc_lnk( psi_uw, 'U', -1. ) ! lateral boundary condition 701 CALL lbc_lnk( psi_vw, 'V', -1. ) 702 ! 703 !!gm CALL iom_put( "psi_eiv_uw", psi_uw ) ! output 704 !!gm CALL iom_put( "psi_eiv_vw", psi_vw ) 705 ! 706 ! !== eiv velocities: calculate and output ==! 707 CALL wrk_alloc( jpi,jpj,jpk, zw3d ) 708 ! 709 zw3d(:,:,jpk) = 0._wp ! bottom value always 0 710 ! 711 DO jk = 1, jpkm1 ! e2u e3u u_eiv = -dk[psi_uw] 712 zw3d(:,:,jk) = ( psi_uw(:,:,jk+1) - psi_uw(:,:,jk) ) / ( e2u(:,:) * e3u_n(:,:,jk) ) 713 END DO 714 CALL iom_put( "uoce_eiv", zw3d ) 715 ! 716 DO jk = 1, jpkm1 ! e1v e3v v_eiv = -dk[psi_vw] 717 zw3d(:,:,jk) = ( psi_vw(:,:,jk+1) - psi_vw(:,:,jk) ) / ( e1v(:,:) * e3v_n(:,:,jk) ) 718 END DO 719 CALL iom_put( "voce_eiv", zw3d ) 720 ! 721 DO jk = 1, jpkm1 ! e1 e2 w_eiv = dk[psix] + dk[psix] 722 DO jj = 2, jpjm1 723 DO ji = fs_2, fs_jpim1 ! vector opt. 724 zw3d(ji,jj,jk) = ( psi_vw(ji,jj,jk) - psi_vw(ji ,jj-1,jk) & 725 & + psi_uw(ji,jj,jk) - psi_uw(ji-1,jj ,jk) ) / e1e2t(ji,jj) 726 END DO 727 END DO 728 END DO 729 CALL lbc_lnk( zw3d, 'T', 1. ) ! lateral boundary condition 730 CALL iom_put( "woce_eiv", zw3d ) 731 ! 732 CALL wrk_dealloc( jpi,jpj,jpk, zw3d ) 733 ! 734 ! 735 IF( lk_diaar5 ) THEN !== eiv heat transport: calculate and output ==! 736 CALL wrk_alloc( jpi,jpj, zw2d ) 737 ! 738 zztmp = 0.5_wp * rau0 * rcp 739 zw2d(:,:) = 0._wp 740 DO jk = 1, jpkm1 741 DO jj = 2, jpjm1 742 DO ji = fs_2, fs_jpim1 ! vector opt. 743 zw2d(ji,jj) = zw2d(ji,jj) + zztmp * ( psi_uw(ji,jj,jk+1) - psi_uw(ji,jj,jk) ) & 744 & * ( tsn (ji,jj,jk,jp_tem) + tsn (ji+1,jj,jk,jp_tem) ) 745 END DO 746 END DO 747 END DO 748 CALL lbc_lnk( zw2d, 'U', -1. ) 749 CALL iom_put( "ueiv_heattr", zw2d ) ! heat transport in i-direction 750 zw2d(:,:) = 0._wp 751 DO jk = 1, jpkm1 752 DO jj = 2, jpjm1 753 DO ji = fs_2, fs_jpim1 ! vector opt. 754 zw2d(ji,jj) = zw2d(ji,jj) + zztmp * ( psi_vw(ji,jj,jk+1) - psi_vw(ji,jj,jk) ) & 755 & * ( tsn (ji,jj,jk,jp_tem) + tsn (ji,jj+1,jk,jp_tem) ) 756 END DO 757 END DO 758 END DO 759 CALL lbc_lnk( zw2d, 'V', -1. ) 760 CALL iom_put( "veiv_heattr", zw2d ) ! heat transport in i-direction 761 ! 762 CALL wrk_dealloc( jpi,jpj, zw2d ) 763 ENDIF 764 ! 765 IF( nn_timing == 1 ) CALL timing_stop( 'ldf_eiv_dia') 766 ! 767 END SUBROUTINE ldf_eiv_dia 188 768 189 769 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.