Changeset 4488 for branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/LDF
- Timestamp:
- 2014-02-06T11:43:09+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r3848 r4488 117 117 CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 118 118 119 zeps = 1.e-20_wp !== Local constant initialization ==! 120 z1_16 = 1.0_wp / 16._wp 121 zm1_g = -1.0_wp / grav 122 zm1_2g = -0.5_wp / grav 123 ! 124 zww(:,:,:) = 0._wp 125 zwz(:,:,:) = 0._wp 126 ! 127 DO jk = 1, jpk !== i- & j-gradient of density ==! 128 DO jj = 1, jpjm1 129 DO ji = 1, fs_jpim1 ! vector opt. 130 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 131 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 132 END DO 133 END DO 134 END DO 135 IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level 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 136 138 # if defined key_vectopt_loop 137 DO jj = 1, 1138 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)139 DO jj = 1, 1 140 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 139 141 # else 140 DO jj = 1, jpjm1141 DO ji = 1, jpim1142 DO jj = 1, jpjm1 143 DO ji = 1, jpim1 142 144 # endif 143 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 144 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 145 END DO 146 END DO 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) 192 !!gm modif to suppress omlmask.... (as in Griffies case) 193 ! ! ! jk must be >= ML level for zf=1. otherwise zf=0. 194 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 195 ! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 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 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 199 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 200 !!gm end modif 201 END DO 202 END DO 203 END DO 204 CALL lbc_lnk( zwz, 'U', -1. ) ; CALL lbc_lnk( zww, 'V', -1. ) ! lateral boundary conditions 205 ! 206 ! !* horizontal Shapiro filter 207 DO jk = 2, jpkm1 208 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 209 DO ji = 2, jpim1 210 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 211 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 212 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 213 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 214 & + 4.* zwz(ji ,jj ,jk) ) 215 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 216 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 217 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 218 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 219 & + 4.* zww(ji,jj ,jk) ) 220 END DO 221 END DO 222 DO jj = 3, jpj-2 ! other rows 223 DO ji = fs_2, fs_jpim1 ! vector opt. 224 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 225 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 226 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 227 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 228 & + 4.* zwz(ji ,jj ,jk) ) 229 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 230 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 231 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 232 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 233 & + 4.* zww(ji,jj ,jk) ) 234 END DO 235 END DO 236 ! !* decrease along coastal boundaries 237 DO jj = 2, jpjm1 238 DO ji = fs_2, fs_jpim1 ! vector opt. 239 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 240 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 241 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 242 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 243 END DO 244 END DO 245 END DO 246 247 248 ! II. slopes at w point | wslpi = mij( d/di( prd ) / d/dz( prd ) 249 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 250 ! 251 DO jk = 2, jpkm1 252 DO jj = 2, jpjm1 253 DO ji = fs_2, fs_jpim1 ! vector opt. 254 ! !* Local vertical density gradient evaluated from N^2 255 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 256 ! !* Slopes at w point 257 ! ! i- & j-gradient of density at w-points 258 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 259 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 260 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 261 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 262 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 263 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * tmask (ji,jj,jk) 264 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 265 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * tmask (ji,jj,jk) 266 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 267 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 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 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 271 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 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) 275 276 !!gm modif to suppress omlmask.... (as in Griffies operator) 277 ! ! ! jk must be >= ML level for zfk=1. otherwise zfk=0. 278 ! zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp ) 279 ! zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. ) 280 ! zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk) 281 ! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 282 !!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 147 404 ENDIF 148 ! 149 zdzr(:,:,1) = 0._wp !== Local vertical density gradient at T-point == ! (evaluated from N^2) 150 DO jk = 2, jpkm1 151 ! ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point 152 ! ! trick: tmask(ik ) = 0 => all pn2 = 0 => zdzr = 0 153 ! ! else tmask(ik+1) = 0 => pn2(ik+1) = 0 => zdzr divides by 1 154 ! ! umask(ik+1) /= 0 => all pn2 /= 0 => zdzr divides by 2 155 ! ! NB: 1/(tmask+1) = (1-.5*tmask) substitute a / by a * ==> faster 156 zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp ) & 157 & * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) ) 158 END DO 159 ! 160 ! !== Slopes just below the mixed layer ==! 161 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml 162 163 164 ! I. slopes at u and v point | uslp = d/di( prd ) / d/dz( prd ) 165 ! =========================== | vslp = d/dj( prd ) / d/dz( prd ) 166 ! 167 DO jk = 2, jpkm1 !* Slopes at u and v points 168 DO jj = 2, jpjm1 169 DO ji = fs_2, fs_jpim1 ! vector opt. 170 ! ! horizontal and vertical density gradient at u- and v-points 171 zau = zgru(ji,jj,jk) / e1u(ji,jj) 172 zav = zgrv(ji,jj,jk) / e2v(ji,jj) 173 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 174 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 175 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 176 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 177 zbu = MIN( zbu, -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau ) ) 178 zbv = MIN( zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav ) ) 179 ! ! uslp and vslp output in zwz and zww, resp. 180 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 181 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 182 zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps ) & 183 & + zfi * uslpml(ji,jj) & 184 & * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) ) & 185 & / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk) 186 zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps ) & 187 & + zfj * vslpml(ji,jj) & 188 & * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) ) & 189 & / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 190 !!gm modif to suppress omlmask.... (as in Griffies case) 191 ! ! ! jk must be >= ML level for zf=1. otherwise zf=0. 192 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 193 ! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 194 ! zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 195 ! zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 196 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 197 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 198 !!gm end modif 199 END DO 200 END DO 201 END DO 202 CALL lbc_lnk( zwz, 'U', -1. ) ; CALL lbc_lnk( zww, 'V', -1. ) ! lateral boundary conditions 203 ! 204 ! !* horizontal Shapiro filter 205 DO jk = 2, jpkm1 206 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 207 DO ji = 2, jpim1 208 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 209 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 210 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 211 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 212 & + 4.* zwz(ji ,jj ,jk) ) 213 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 214 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 215 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 216 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 217 & + 4.* zww(ji,jj ,jk) ) 218 END DO 219 END DO 220 DO jj = 3, jpj-2 ! other rows 221 DO ji = fs_2, fs_jpim1 ! vector opt. 222 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 223 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 224 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 225 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 226 & + 4.* zwz(ji ,jj ,jk) ) 227 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 228 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 229 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 230 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 231 & + 4.* zww(ji,jj ,jk) ) 232 END DO 233 END DO 234 ! !* decrease along coastal boundaries 235 DO jj = 2, jpjm1 236 DO ji = fs_2, fs_jpim1 ! vector opt. 237 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 238 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 239 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 240 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 241 END DO 242 END DO 243 END DO 244 245 246 ! II. slopes at w point | wslpi = mij( d/di( prd ) / d/dz( prd ) 247 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 248 ! 249 DO jk = 2, jpkm1 250 DO jj = 2, jpjm1 251 DO ji = fs_2, fs_jpim1 ! vector opt. 252 ! !* Local vertical density gradient evaluated from N^2 253 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 254 ! !* Slopes at w point 255 ! ! i- & j-gradient of density at w-points 256 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 257 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 258 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 259 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 260 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 261 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * tmask (ji,jj,jk) 262 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 263 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * tmask (ji,jj,jk) 264 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 265 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 266 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai ) ) 267 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj ) ) 268 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 269 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 270 zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp ) 271 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * tmask(ji,jj,jk) 272 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * tmask(ji,jj,jk) 273 274 !!gm modif to suppress omlmask.... (as in Griffies operator) 275 ! ! ! jk must be >= ML level for zfk=1. otherwise zfk=0. 276 ! zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp ) 277 ! zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. ) 278 ! zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk) 279 ! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 280 !!gm end modif 281 END DO 282 END DO 283 END DO 284 CALL lbc_lnk( zwz, 'T', -1. ) ; CALL lbc_lnk( zww, 'T', -1. ) ! lateral boundary conditions 285 ! 286 ! !* horizontal Shapiro filter 287 DO jk = 2, jpkm1 288 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 289 DO ji = 2, jpim1 290 zcofw = tmask(ji,jj,jk) * z1_16 291 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 292 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 293 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 294 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 295 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 296 297 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 298 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 299 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 300 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 301 & + 4.* zww(ji ,jj ,jk) ) * zcofw 302 END DO 303 END DO 304 DO jj = 3, jpj-2 ! other rows 305 DO ji = fs_2, fs_jpim1 ! vector opt. 306 zcofw = tmask(ji,jj,jk) * z1_16 307 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 308 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 309 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 310 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 311 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 312 313 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 314 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 315 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 316 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 317 & + 4.* zww(ji ,jj ,jk) ) * zcofw 318 END DO 319 END DO 320 ! !* decrease along coastal boundaries 321 DO jj = 2, jpjm1 322 DO ji = fs_2, fs_jpim1 ! vector opt. 323 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 324 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 325 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 326 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 327 END DO 328 END DO 329 END DO 330 331 ! III. Specific grid points 332 ! =========================== 333 ! 334 IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN ! ORCA_R4 configuration: horizontal diffusion in specific area 335 ! ! Gibraltar Strait 336 ij0 = 50 ; ij1 = 53 337 ii0 = 69 ; ii1 = 71 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 338 ij0 = 51 ; ij1 = 53 339 ii0 = 68 ; ii1 = 71 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 340 ii0 = 69 ; ii1 = 71 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 341 ii0 = 69 ; ii1 = 71 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 342 ! 343 ! ! Mediterrannean Sea 344 ij0 = 49 ; ij1 = 56 345 ii0 = 71 ; ii1 = 90 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 346 ij0 = 50 ; ij1 = 56 347 ii0 = 70 ; ii1 = 90 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 348 ii0 = 71 ; ii1 = 90 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 349 ii0 = 71 ; ii1 = 90 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp 350 ENDIF 351 352 353 ! IV. Lateral boundary conditions 354 ! =============================== 355 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) 356 CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. ) 357 358 359 IF(ln_ctl) THEN 360 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 361 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) 362 ENDIF 363 ! 405 364 406 CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 365 407 ! … … 783 825 wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp 784 826 785 !!gm I no longer understand this..... 786 IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 827 IF( ln_traldf_hor .OR. ln_dynldf_hor ) THEN 787 828 IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 788 829 … … 796 837 DO jj = 2, jpjm1 797 838 DO ji = fs_2, fs_jpim1 ! vector opt. 798 uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept (ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)799 vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept (ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)800 wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw (ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5801 wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw (ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5839 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 802 843 END DO 803 844 END DO
Note: See TracChangeset
for help on using the changeset viewer.