- Timestamp:
- 2018-04-23T10:44:07+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/LDF/ldfdyn.F90
r9169 r9490 17 17 USE dom_oce ! ocean space and time domain 18 18 USE phycst ! physical constants 19 USE ldfslp ! lateral diffusion: slopes of mixing orientation 19 20 USE ldfc1d_c2d ! lateral diffusion: 1D and 2D cases 20 21 ! … … 31 32 PUBLIC ldf_dyn ! called by step.F90 32 33 33 ! 34 ! !!* Namelist namdyn_ldf : lateral mixing on momentum * 34 35 LOGICAL , PUBLIC :: ln_dynldf_NONE !: No operator (i.e. no explicit diffusion) 35 36 LOGICAL , PUBLIC :: ln_dynldf_lap !: laplacian operator … … 37 38 LOGICAL , PUBLIC :: ln_dynldf_lev !: iso-level direction 38 39 LOGICAL , PUBLIC :: ln_dynldf_hor !: horizontal (geopotential) direction 39 LOGICAL , PUBLIC :: ln_dynldf_iso !: iso-neutral direction 40 ! LOGICAL , PUBLIC :: ln_dynldf_iso !: iso-neutral direction (see ldfslp) 40 41 INTEGER , PUBLIC :: nn_ahm_ijk_t !: choice of time & space variations of the lateral eddy viscosity coef. 41 REAL(wp), PUBLIC :: rn_ahm_0 !: lateral laplacian eddy viscosity [m2/s] 42 REAL(wp), PUBLIC :: rn_ahm_b !: lateral laplacian background eddy viscosity [m2/s] 43 REAL(wp), PUBLIC :: rn_bhm_0 !: lateral bilaplacian eddy viscosity [m4/s] 44 !! If nn_ahm_ijk_t = 32 a time and space varying Smagorinsky viscosity 45 !! will be computed. 46 REAL(wp), PUBLIC :: rn_csmc !: Smagorinsky constant of proportionality 47 REAL(wp), PUBLIC :: rn_minfac !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 48 REAL(wp), PUBLIC :: rn_maxfac !: Multiplicative factor of theorectical maximum Smagorinsky viscosity 49 50 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 51 52 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy diffusivity coef. at U- and V-points [m2/s or m4/s] 42 ! ! time invariant coefficients: aht = 1/2 Ud*Ld (lap case) 43 ! ! bht = 1/12 Ud*Ld^3 (blp case) 44 REAL(wp), PUBLIC :: rn_Uv !: lateral viscous velocity [m/s] 45 REAL(wp), PUBLIC :: rn_Lv !: lateral viscous length [m] 46 ! ! Smagorinsky viscosity (nn_ahm_ijk_t = 32) 47 REAL(wp), PUBLIC :: rn_csmc !: Smagorinsky constant of proportionality 48 REAL(wp), PUBLIC :: rn_minfac !: Multiplicative factor of theorectical minimum Smagorinsky viscosity 49 REAL(wp), PUBLIC :: rn_maxfac !: Multiplicative factor of theorectical maximum Smagorinsky viscosity 50 ! ! iso-neutral laplacian (ln_dynldf_lap=ln_dynldf_iso=T) 51 REAL(wp), PUBLIC :: rn_ahm_b !: lateral laplacian background eddy viscosity [m2/s] 52 53 ! !!* Parameter to control the type of lateral viscous operator 54 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 !: error in setting the operator 55 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 !: without operator (i.e. no lateral viscous trend) 56 ! !! laplacian ! bilaplacian ! 57 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 !: iso-level operator 58 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 !: iso-neutral or geopotential operator 59 ! 60 INTEGER , PUBLIC :: nldf_dyn !: type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 61 LOGICAL , PUBLIC :: l_ldfdyn_time !: flag for time variation of the lateral eddy viscosity coef. 62 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ahmt, ahmf !: eddy viscosity coef. at T- and F-points [m2/s or m4/s] 53 64 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dtensq !: horizontal tension squared (Smagorinsky only) 54 65 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dshesq !: horizontal shearing strain squared (Smagorinsky only) 55 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: esqt, esqf !: Square of the local gridscale (e1e2/(e1+e2))**2 56 67 57 REAL(wp) :: r1_ 12 = 1._wp / 12._wp ! =1/1268 REAL(wp) :: r1_2 = 0.5_wp ! =1/2 58 69 REAL(wp) :: r1_4 = 0.25_wp ! =1/4 59 70 REAL(wp) :: r1_8 = 0.125_wp ! =1/8 71 REAL(wp) :: r1_12 = 1._wp / 12._wp ! =1/12 60 72 REAL(wp) :: r1_288 = 1._wp / 288._wp ! =1/( 12^2 * 2 ) 61 73 … … 92 104 !! or L^4|D|/8 bilaplacian operator ) 93 105 !!---------------------------------------------------------------------- 94 INTEGER :: ji, jj, jk ! dummy loop indices 95 INTEGER :: ierr, inum, ios ! local integer 96 REAL(wp) :: zah0 ! local scalar 97 ! 106 INTEGER :: ji, jj, jk ! dummy loop indices 107 INTEGER :: ioptio, ierr, inum, ios, inn ! local integer 108 REAL(wp) :: zah0, zah_max, zUfac ! local scalar 109 CHARACTER(len=5) :: cl_Units ! units (m2/s or m4/s) 110 !! 98 111 NAMELIST/namdyn_ldf/ ln_dynldf_NONE, ln_dynldf_lap, ln_dynldf_blp, & ! type of operator 99 & ln_dynldf_lev , ln_dynldf_hor, ln_dynldf_iso, & ! acting direction of the operator100 & nn_ahm_ijk_t , rn_ahm_0, rn_ahm_b, rn_bhm_0, & ! lateral eddy coefficient101 & rn_csmc , rn_minfac, rn_maxfac! Smagorinsky settings112 & ln_dynldf_lev , ln_dynldf_hor, ln_dynldf_iso, & ! acting direction of the operator 113 & nn_ahm_ijk_t , rn_Uv , rn_Lv, rn_ahm_b, & ! lateral eddy coefficient 114 & rn_csmc , rn_minfac , rn_maxfac ! Smagorinsky settings 102 115 !!---------------------------------------------------------------------- 103 116 ! … … 129 142 WRITE(numout,*) ' coefficients :' 130 143 WRITE(numout,*) ' type of time-space variation nn_ahm_ijk_t = ', nn_ahm_ijk_t 131 WRITE(numout,*) ' lateral laplacian eddy viscosity rn_ahm_0 = ', rn_ahm_0, ' m2/s' 132 WRITE(numout,*) ' background viscosity (iso case) rn_ahm_b = ', rn_ahm_b, ' m2/s' 133 WRITE(numout,*) ' lateral bilaplacian eddy viscosity rn_bhm_0 = ', rn_bhm_0, ' m4/s' 144 WRITE(numout,*) ' lateral viscous velocity (if cst) rn_Uv = ', rn_Uv, ' m/s' 145 WRITE(numout,*) ' lateral viscous length (if cst) rn_Lv = ', rn_Lv, ' m' 146 WRITE(numout,*) ' background viscosity (iso-lap case) rn_ahm_b = ', rn_ahm_b, ' m2/s' 147 ! 134 148 WRITE(numout,*) ' Smagorinsky settings (nn_ahm_ijk_t = 32) :' 135 149 WRITE(numout,*) ' Smagorinsky coefficient rn_csmc = ', rn_csmc 136 WRITE(numout,*) ' factor multiplier for theorectical lower limit for ' 137 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_minfac = ', rn_minfac 138 WRITE(numout,*) ' factor multiplier for theorectical lower upper for ' 139 WRITE(numout,*) ' Smagorinsky eddy visc (def. 1.0) rn_maxfac = ', rn_maxfac 150 WRITE(numout,*) ' factor multiplier for eddy visc.' 151 WRITE(numout,*) ' lower limit (default 1.0) rn_minfac = ', rn_minfac 152 WRITE(numout,*) ' upper limit (default 1.0) rn_maxfac = ', rn_maxfac 140 153 ENDIF 141 154 142 ! ! Parameter control 155 ! 156 ! !== type of lateral operator used ==! (set nldf_dyn) 157 ! !=====================================! 158 ! 159 nldf_dyn = np_ERROR 160 ioptio = 0 161 IF( ln_dynldf_NONE ) THEN ; nldf_dyn = np_no_ldf ; ioptio = ioptio + 1 ; ENDIF 162 IF( ln_dynldf_lap ) THEN ; ioptio = ioptio + 1 ; ENDIF 163 IF( ln_dynldf_blp ) THEN ; ioptio = ioptio + 1 ; ENDIF 164 IF( ioptio /= 1 ) CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 operator options (NONE/lap/blp)' ) 165 ! 166 IF(.NOT.ln_dynldf_NONE ) THEN !== direction ==>> type of operator ==! 167 ioptio = 0 168 IF( ln_dynldf_lev ) ioptio = ioptio + 1 169 IF( ln_dynldf_hor ) ioptio = ioptio + 1 170 IF( ln_dynldf_iso ) ioptio = ioptio + 1 171 IF( ioptio /= 1 ) CALL ctl_stop( 'dyn_ldf_init: use ONE of the 3 direction options (level/hor/iso)' ) 172 ! 173 ! ! Set nldf_dyn, the type of lateral diffusion, from ln_dynldf_... logicals 174 ierr = 0 175 IF( ln_dynldf_lap ) THEN ! laplacian operator 176 IF( ln_zco ) THEN ! z-coordinate 177 IF ( ln_dynldf_lev ) nldf_dyn = np_lap ! iso-level = horizontal (no rotation) 178 IF ( ln_dynldf_hor ) nldf_dyn = np_lap ! iso-level = horizontal (no rotation) 179 IF ( ln_dynldf_iso ) nldf_dyn = np_lap_i ! iso-neutral ( rotation) 180 ENDIF 181 IF( ln_zps ) THEN ! z-coordinate with partial step 182 IF ( ln_dynldf_lev ) nldf_dyn = np_lap ! iso-level (no rotation) 183 IF ( ln_dynldf_hor ) nldf_dyn = np_lap ! iso-level (no rotation) 184 IF ( ln_dynldf_iso ) nldf_dyn = np_lap_i ! iso-neutral ( rotation) 185 ENDIF 186 IF( ln_sco ) THEN ! s-coordinate 187 IF ( ln_dynldf_lev ) nldf_dyn = np_lap ! iso-level = horizontal (no rotation) 188 IF ( ln_dynldf_hor ) nldf_dyn = np_lap_i ! horizontal ( rotation) 189 IF ( ln_dynldf_iso ) nldf_dyn = np_lap_i ! iso-neutral ( rotation) 190 ENDIF 191 ENDIF 192 ! 193 IF( ln_dynldf_blp ) THEN ! bilaplacian operator 194 IF( ln_zco ) THEN ! z-coordinate 195 IF( ln_dynldf_lev ) nldf_dyn = np_blp ! iso-level = horizontal (no rotation) 196 IF( ln_dynldf_hor ) nldf_dyn = np_blp ! iso-level = horizontal (no rotation) 197 IF( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 198 ENDIF 199 IF( ln_zps ) THEN ! z-coordinate with partial step 200 IF( ln_dynldf_lev ) nldf_dyn = np_blp ! iso-level (no rotation) 201 IF( ln_dynldf_hor ) nldf_dyn = np_blp ! iso-level (no rotation) 202 IF( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 203 ENDIF 204 IF( ln_sco ) THEN ! s-coordinate 205 IF( ln_dynldf_lev ) nldf_dyn = np_blp ! iso-level (no rotation) 206 IF( ln_dynldf_hor ) ierr = 2 ! horizontal ( rotation) 207 IF( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 208 ENDIF 209 ENDIF 210 ! 211 IF( ierr == 2 ) CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 212 ! 213 IF( nldf_dyn == np_lap_i ) l_ldfslp = .TRUE. ! rotation require the computation of the slopes 214 ! 215 ENDIF 216 ! 217 IF(lwp) THEN 218 WRITE(numout,*) 219 SELECT CASE( nldf_dyn ) 220 CASE( np_no_ldf ) ; WRITE(numout,*) ' ==>>> NO lateral viscosity' 221 CASE( np_lap ) ; WRITE(numout,*) ' ==>>> iso-level laplacian operator' 222 CASE( np_lap_i ) ; WRITE(numout,*) ' ==>>> rotated laplacian operator with iso-level background' 223 CASE( np_blp ) ; WRITE(numout,*) ' ==>>> iso-level bi-laplacian operator' 224 END SELECT 225 WRITE(numout,*) 226 ENDIF 227 228 ! 229 ! !== Space/time variation of eddy coefficients ==! 230 ! !=================================================! 231 ! 232 l_ldfdyn_time = .FALSE. ! no time variation except in case defined below 233 ! 143 234 IF( ln_dynldf_NONE ) THEN 144 235 IF(lwp) WRITE(numout,*) ' ==>>> No viscous operator selected. ahmt and ahmf are not allocated' 145 l_ldfdyn_time = .FALSE.146 236 RETURN 147 ENDIF 148 ! 149 IF( ln_dynldf_blp .AND. ln_dynldf_iso ) THEN ! iso-neutral bilaplacian not implemented 150 CALL ctl_stop( 'dyn_ldf_init: iso-neutral bilaplacian not coded yet' ) 151 ENDIF 152 153 ! ... Space/Time variation of eddy coefficients 154 ! ! allocate the ahm arrays 155 ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 156 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 157 ! 158 ahmt(:,:,jpk) = 0._wp ! last level always 0 159 ahmf(:,:,jpk) = 0._wp 160 ! 161 ! ! value of eddy mixing coef. 162 IF ( ln_dynldf_lap ) THEN ; zah0 = rn_ahm_0 ! laplacian operator 163 ELSEIF( ln_dynldf_blp ) THEN ; zah0 = ABS( rn_bhm_0 ) ! bilaplacian operator 164 ELSE ! NO viscous operator 165 CALL ctl_warn( 'ldf_dyn_init: No lateral viscous operator used ' ) 166 ENDIF 167 ! 168 l_ldfdyn_time = .FALSE. ! no time variation except in case defined below 169 ! 170 IF( ln_dynldf_lap .OR. ln_dynldf_blp ) THEN ! only if a lateral diffusion operator is used 171 ! 172 SELECT CASE( nn_ahm_ijk_t ) ! Specification of space time variations of ahmt, ahmf 237 ! 238 ELSE !== a lateral diffusion operator is used ==! 239 ! 240 ! ! allocate the ahm arrays 241 ALLOCATE( ahmt(jpi,jpj,jpk) , ahmf(jpi,jpj,jpk) , STAT=ierr ) 242 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate arrays') 243 ! 244 ahmt(:,:,jpk) = 0._wp ! last level always 0 245 ahmf(:,:,jpk) = 0._wp 246 ! 247 ! ! value of lap/blp eddy mixing coef. 248 IF( ln_dynldf_lap ) THEN ; zUfac = r1_2 *rn_Uv ; inn = 1 ; cl_Units = ' m2/s' ! laplacian 249 ELSEIF( ln_dynldf_blp ) THEN ; zUfac = r1_12*rn_Uv ; inn = 3 ; cl_Units = ' m4/s' ! bilaplacian 250 ENDIF 251 zah0 = zUfac * rn_Lv**inn ! mixing coefficient 252 zah_max = zUfac * (ra*rad)**inn ! maximum reachable coefficient (value at the Equator) 253 ! 254 SELECT CASE( nn_ahm_ijk_t ) !* Specification of space-time variations of ahmt, ahmf 173 255 ! 174 256 CASE( 0 ) !== constant ==! 175 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = constant '176 ahmt(:,:, :) = zah0 * tmask(:,:,:)177 ahmf(:,:, :) = zah0 * fmask(:,:,:)257 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity. = constant = ', zah0, cl_Units 258 ahmt(:,:,1:jpkm1) = zah0 259 ahmf(:,:,1:jpkm1) = zah0 178 260 ! 179 261 CASE( 10 ) !== fixed profile ==! 180 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = F( depth )' 181 ahmt(:,:,1) = zah0 * tmask(:,:,1) ! constant surface value 182 ahmf(:,:,1) = zah0 * fmask(:,:,1) 183 CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 262 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( depth )' 263 IF(lwp) WRITE(numout,*) ' surface viscous coef. = constant = ', zah0, cl_Units 264 ahmt(:,:,1) = zah0 ! constant surface value 265 ahmf(:,:,1) = zah0 266 CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) 184 267 ! 185 268 CASE ( -20 ) !== fixed horizontal shape read in file ==! 186 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef.= F(i,j) read in eddy_viscosity.nc file'269 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F(i,j) read in eddy_viscosity.nc file' 187 270 CALL iom_open( 'eddy_viscosity_2D.nc', inum ) 188 271 CALL iom_get ( inum, jpdom_data, 'ahmt_2d', ahmt(:,:,1) ) 189 272 CALL iom_get ( inum, jpdom_data, 'ahmf_2d', ahmf(:,:,1) ) 190 273 CALL iom_close( inum ) 191 !!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ???192 !! do we introduce a scaling by the max value of the array, and then multiply by zah0 ????193 !! better: check that the max is <=1 i.e. it is a shape from 0 to 1, not a coef that has physical dimension194 274 DO jk = 2, jpkm1 195 ahmt(:,:,jk) = ahmt(:,:,1) * tmask(:,:,jk)196 ahmf(:,:,jk) = ahmf(:,:,1) * fmask(:,:,jk)275 ahmt(:,:,jk) = ahmt(:,:,1) 276 ahmf(:,:,jk) = ahmf(:,:,1) 197 277 END DO 198 278 ! 199 279 CASE( 20 ) !== fixed horizontal shape ==! 200 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)' 201 IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor 202 IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor^3 280 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( e1, e2 ) or F( e1^3, e2^3 ) (lap. or blp. case)' 281 IF(lwp) WRITE(numout,*) ' using a fixed viscous velocity = ', rn_Uv ,' m/s and Lv = Max(e1,e2)' 282 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)' 283 CALL ldf_c2d( 'DYN', zUfac , inn , ahmt, ahmf ) ! surface value proportional to scale factor^inn 203 284 ! 204 285 CASE( -30 ) !== fixed 3D shape read in file ==! 205 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = F(i,j,k) read in eddy_diffusivity_3D.nc file'286 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F(i,j,k) read in eddy_viscosity_3D.nc file' 206 287 CALL iom_open( 'eddy_viscosity_3D.nc', inum ) 207 288 CALL iom_get ( inum, jpdom_data, 'ahmt_3d', ahmt ) 208 289 CALL iom_get ( inum, jpdom_data, 'ahmf_3d', ahmf ) 209 290 CALL iom_close( inum ) 210 !!gm Question : info for LAP or BLP case to take into account the SQRT in the bilaplacian case ????211 !! do we introduce a scaling by the max value of the array, and then multiply by zah0 ????212 DO jk = 1, jpkm1213 ahmt(:,:,jk) = ahmt(:,:,jk) * tmask(:,:,jk)214 ahmf(:,:,jk) = ahmf(:,:,jk) * fmask(:,:,jk)215 END DO216 291 ! 217 292 CASE( 30 ) !== fixed 3D shape ==! 218 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef.= F( latitude, longitude, depth )'219 IF( ln_dynldf_lap ) CALL ldf_c2d( 'DYN', 'LAP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor220 IF( ln_dynldf_blp ) CALL ldf_c2d( 'DYN', 'BLP', zah0, ahmt, ahmf ) ! surface value proportional to scale factor221 ! ! reduction with depth222 CALL ldf_c1d( 'DYN', r1_4, ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf )293 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( latitude, longitude, depth )' 294 IF(lwp) WRITE(numout,*) ' using a fixed viscous velocity = ', rn_Uv ,' m/s and Ld = Max(e1,e2)' 295 IF(lwp) WRITE(numout,*) ' maximum reachable coefficient (at the Equator) = ', zah_max, cl_Units, ' for e1=1°)' 296 CALL ldf_c2d( 'DYN', zUfac , inn , ahmt, ahmf ) ! surface value proportional to scale factor^inn 297 CALL ldf_c1d( 'DYN', ahmt(:,:,1), ahmf(:,:,1), ahmt, ahmf ) ! reduction with depth 223 298 ! 224 299 CASE( 31 ) !== time varying 3D field ==! 225 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef.= F( latitude, longitude, depth , time )'226 IF(lwp) WRITE(numout,*) ' proportional to the velocity : |u|e/12 or |u|e^3/12'300 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( latitude, longitude, depth , time )' 301 IF(lwp) WRITE(numout,*) ' proportional to the local velocity : 1/2 |u|e (lap) or 1/12 |u|e^3 (blp)' 227 302 ! 228 303 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 229 304 ! 230 305 CASE( 32 ) !== time varying 3D field ==! 231 IF(lwp) WRITE(numout,*) ' ==>>> momentum mixing coef. = F( latitude, longitude, depth , time )' 232 IF(lwp) WRITE(numout,*) ' proportional to the local deformation rate and gridscale (Smagorinsky)' 233 IF(lwp) WRITE(numout,*) ' : L^2|D| or L^4|D|/8' 306 IF(lwp) WRITE(numout,*) ' ==>>> eddy viscosity = F( latitude, longitude, depth , time )' 307 IF(lwp) WRITE(numout,*) ' proportional to the local deformation rate and gridscale (Smagorinsky)' 234 308 ! 235 309 l_ldfdyn_time = .TRUE. ! will be calculated by call to ldf_dyn routine in step.F90 236 310 ! 237 ! allocate arrays used in ldf_dyn.311 ! ! allocate arrays used in ldf_dyn. 238 312 ALLOCATE( dtensq(jpi,jpj) , dshesq(jpi,jpj) , esqt(jpi,jpj) , esqf(jpi,jpj) , STAT=ierr ) 239 313 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'ldf_dyn_init: failed to allocate Smagorinsky arrays') 240 314 ! 241 ! Set local gridscale values 242 DO jj = 2, jpjm1 315 DO jj = 2, jpjm1 ! Set local gridscale values 243 316 DO ji = fs_2, fs_jpim1 244 317 esqt(ji,jj) = ( e1e2t(ji,jj) /( e1t(ji,jj) + e2t(ji,jj) ) )**2 … … 251 324 END SELECT 252 325 ! 253 IF( ln_dynldf_blp .AND. .NOT. l_ldfdyn_time ) THEN ! bilapcian and no time variation: 254 ahmt(:,:,:) = SQRT( ahmt(:,:,:) ) ! take the square root of the coefficient 255 ahmf(:,:,:) = SQRT( ahmf(:,:,:) ) 326 IF( .NOT.l_ldfdyn_time ) THEN !* No time variation 327 IF( ln_dynldf_lap ) THEN ! laplacian operator (mask only) 328 ahmt(:,:,1:jpkm1) = ahmt(:,:,1:jpkm1) * tmask(:,:,1:jpkm1) 329 ahmf(:,:,1:jpkm1) = ahmf(:,:,1:jpkm1) * fmask(:,:,1:jpkm1) 330 ELSEIF( ln_dynldf_blp ) THEN ! bilaplacian operator (square root + mask) 331 ahmt(:,:,1:jpkm1) = SQRT( ahmt(:,:,1:jpkm1) ) * tmask(:,:,1:jpkm1) 332 ahmf(:,:,1:jpkm1) = SQRT( ahmf(:,:,1:jpkm1) ) * fmask(:,:,1:jpkm1) 333 ENDIF 256 334 ENDIF 257 335 ! … … 341 419 & - ( vb(ji,jj,jk) * r1_e1v(ji,jj) - vb(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) & 342 420 & * r1_e2t(ji,jj) * e1t(ji,jj) ) * tmask(ji,jj,jk) 343 dtensq(ji,jj) = zdb *zdb421 dtensq(ji,jj) = zdb * zdb 344 422 END DO 345 423 END DO … … 351 429 & + ( vb(ji+1,jj,jk) * r1_e2v(ji+1,jj) - vb(ji,jj,jk) * r1_e2v(ji,jj) ) & 352 430 & * r1_e1f(ji,jj) * e2f(ji,jj) ) * fmask(ji,jj,jk) 353 dshesq(ji,jj) = zdb *zdb431 dshesq(ji,jj) = zdb * zdb 354 432 END DO 355 433 END DO … … 385 463 ! 386 464 IF( ln_dynldf_blp ) THEN ! bilaplacian operator : sqrt( (C_smag/pi)^2 L^4 |D|/8) 387 ! = sqrt( A_lap_smag L^2/8 ) 388 ! stability limits already applied to laplacian values 389 ! effective default limits are |U|L^3/12 < B_hm < L^4/(32*2dt) 390 ! 465 ! ! = sqrt( A_lap_smag L^2/8 ) 466 ! ! stability limits already applied to laplacian values 467 ! ! effective default limits are 1/12 |U|L^3 < B_hm < 1//(32*2dt) L^4 391 468 DO jk = 1, jpkm1 392 !393 469 DO jj = 2, jpjm1 394 470 DO ji = fs_2, fs_jpim1 395 ! 396 ahmt(ji,jj,jk) = sqrt( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 397 ahmf(ji,jj,jk) = sqrt( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 398 ! 471 ahmt(ji,jj,jk) = SQRT( r1_8 * esqt(ji,jj) * ahmt(ji,jj,jk) ) 472 ahmf(ji,jj,jk) = SQRT( r1_8 * esqf(ji,jj) * ahmf(ji,jj,jk) ) 399 473 END DO 400 474 END DO
Note: See TracChangeset
for help on using the changeset viewer.
