- Timestamp:
- 2011-10-24T11:29:46+02:00 (13 years ago)
- Location:
- branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/LDF
- Files:
-
- 3 edited
-
ldfslp.F90 (modified) (35 diffs)
-
ldftra.F90 (modified) (2 diffs)
-
ldftra_oce.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r2772 r2980 46 46 ! !! Griffies operator 47 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wslp2 !: wslp**2 from Griffies quarter cells 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi_g, triadj_g !: skew flux slopes relative to geopotentials 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi_g, triadj_g !: skew flux slopes relative to geopotentials 49 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi , triadj !: isoneutral slopes relative to model-coordinate 50 50 … … 58 58 59 59 ! Workspace arrays for ldf_slp_grif. These could be replaced by several 3D and 2D workspace 60 ! arrays from the wrk_nemo module with a bit of code re-writing. The 4D workspace 60 ! arrays from the wrk_nemo module with a bit of code re-writing. The 4D workspace 61 61 ! arrays can't be used here because of the zero-indexing of some of the ranks. ARPDBG. 62 62 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: zdzrho , zdyrho, zdxrho ! Horizontal and vertical density gradients … … 93 93 !!---------------------------------------------------------------------- 94 94 !! *** ROUTINE ldf_slp *** 95 !! 95 !! 96 96 !! ** Purpose : Compute the slopes of neutral surface (slope of isopycnal 97 97 !! surfaces referenced locally) (ln_traldf_iso=T). 98 !! 99 !! ** Method : The slope in the i-direction is computed at U- and 100 !! W-points (uslp, wslpi) and the slope in the j-direction is 98 !! 99 !! ** Method : The slope in the i-direction is computed at U- and 100 !! W-points (uslp, wslpi) and the slope in the j-direction is 101 101 !! computed at V- and W-points (vslp, wslpj). 102 102 !! They are bounded by 1/100 over the whole ocean, and within the … … 112 112 !! bottom slope (ln_sco=T) at level jpk in inildf] 113 113 !! 114 !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and j-slopes 114 !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and j-slopes 115 115 !! of now neutral surfaces at u-, w- and v- w-points, resp. 116 116 !!---------------------------------------------------------------------- … … 127 127 INTEGER :: ii0, ii1, iku ! temporary integer 128 128 INTEGER :: ij0, ij1, ikv ! temporary integer 129 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16 129 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16, zcofw ! local scalars 130 130 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 131 131 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - … … 148 148 DO jj = 1, jpjm1 149 149 DO ji = 1, fs_jpim1 ! vector opt. 150 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 151 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 150 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 151 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 152 152 END DO 153 153 END DO 154 154 END DO 155 155 IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level 156 # if defined key_vectopt_loop 156 # if defined key_vectopt_loop 157 157 DO jj = 1, 1 158 158 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) … … 161 161 DO ji = 1, jpim1 162 162 # endif 163 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 164 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 163 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 164 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 165 165 END DO 166 166 END DO … … 181 181 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml 182 182 183 183 184 184 ! I. slopes at u and v point | uslp = d/di( prd ) / d/dz( prd ) 185 185 ! =========================== | vslp = d/dj( prd ) / d/dz( prd ) 186 ! 186 ! 187 187 DO jk = 2, jpkm1 !* Slopes at u and v points 188 188 DO jj = 2, jpjm1 … … 225 225 DO jk = 2, jpkm1 226 226 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 227 DO ji = 2, jpim1 227 DO ji = 2, jpim1 228 228 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 229 229 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & … … 266 266 ! II. slopes at w point | wslpi = mij( d/di( prd ) / d/dz( prd ) 267 267 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 268 ! 268 ! 269 269 DO jk = 2, jpkm1 270 270 DO jj = 2, jpjm1 … … 308 308 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 309 309 DO ji = 2, jpim1 310 zcofw = tmask(ji,jj,jk) * z1_16 310 311 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 311 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &312 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &313 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &314 & + 4.* zwz(ji ,jj ,jk) ) * z1_16 * tmask(ji,jj,jk)312 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 313 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 314 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 315 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 315 316 316 317 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 317 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &318 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &319 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &320 & + 4.* zww(ji ,jj ,jk) ) * z1_16 * tmask(ji,jj,jk)321 END DO 322 END DO 318 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 319 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 320 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 321 & + 4.* zww(ji ,jj ,jk) ) * zcofw 322 END DO 323 END DO 323 324 DO jj = 3, jpj-2 ! other rows 324 325 DO ji = fs_2, fs_jpim1 ! vector opt. 326 zcofw = tmask(ji,jj,jk) * z1_16 325 327 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 326 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &327 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &328 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &329 & + 4.* zwz(ji ,jj ,jk) ) * z1_16 * tmask(ji,jj,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 330 332 331 333 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 332 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &333 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &334 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &335 & + 4.* zww(ji ,jj ,jk) ) * z1_16 * tmask(ji,jj,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 336 338 END DO 337 339 END DO … … 346 348 END DO 347 349 END DO 348 349 ! III. Specific grid points 350 ! =========================== 351 ! 350 351 ! III. Specific grid points 352 ! =========================== 353 ! 352 354 IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN ! ORCA_R4 configuration: horizontal diffusion in specific area 353 355 ! ! Gibraltar Strait … … 368 370 ENDIF 369 371 370 ! IV. Lateral boundary conditions 372 373 ! IV. Lateral boundary conditions 371 374 ! =============================== 372 375 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) … … 382 385 ! 383 386 END SUBROUTINE ldf_slp 384 387 385 388 386 389 SUBROUTINE ldf_slp_grif ( kt ) … … 390 393 !! ** Purpose : Compute the squared slopes of neutral surfaces (slope 391 394 !! of iso-pycnal surfaces referenced locally) (ln_traldf_grif=T) 392 !! at W-points using the Griffies quarter-cells. 393 !! 394 !! ** Method : calculates alpha and beta at T-points 395 !! at W-points using the Griffies quarter-cells. 396 !! 397 !! ** Method : calculates alpha and beta at T-points 395 398 !! 396 399 !! ** Action : - triadi_g, triadj_g T-pts i- and j-slope triads relative to geopot. (used for eiv) … … 399 402 !!---------------------------------------------------------------------- 400 403 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 401 USE oce , ONLY: zdit => ua , zdis => va ! (ua,va) used as workspace 402 USE oce , ONLY: zdjt => ta , zdjs => sa ! (ta,sa) used as workspace 403 USE wrk_nemo, ONLY: zdkt => wrk_3d_2 , zdks => wrk_3d_3 ! 3D workspace 404 USE wrk_nemo, ONLY: zalpha => wrk_3d_4 , zbeta => wrk_3d_5 ! alpha, beta at T points, at depth fsgdept 404 USE oce , ONLY: zalbet => ua ! use ua as workspace 405 405 USE wrk_nemo, ONLY: z1_mlbw => wrk_2d_1 406 ! 407 INTEGER, INTENT( in ) :: kt ! ocean time-step index408 ! 406 !! 407 INTEGER, INTENT( in ) :: kt ! ocean time-step index 408 !! 409 409 INTEGER :: ji, jj, jk, jl, ip, jp, kp ! dummy loop indices 410 INTEGER :: iku, ikv 411 REAL(wp) :: zfacti, zfactj , zatempw,zatempu,zatempv! local scalars412 REAL(wp) :: z bu, zbv, zbti, zbtj ! - -410 INTEGER :: iku, ikv ! local integer 411 REAL(wp) :: zfacti, zfactj ! local scalars 412 REAL(wp) :: zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj 413 413 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_lim2, zti_g_raw, zti_g_lim 414 414 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_lim2, ztj_g_raw, ztj_g_lim 415 415 REAL(wp) :: zdzrho_raw 416 REAL(wp) :: zbeta0 416 417 !!---------------------------------------------------------------------- 417 418 … … 424 425 !--------------------------------! 425 426 ! 426 CALL eos_alpbet( tsb, zalpha, zbeta ) !== before thermal and haline expension coeff. at T-points ==! 427 ! 428 DO jk = 1, jpkm1 !== before lateral T & S gradients at T-level jk ==! 429 DO jj = 1, jpjm1 430 DO ji = 1, fs_jpim1 ! vector opt. 431 zdit(ji,jj,jk) = ( tb(ji+1,jj,jk) - tb(ji,jj,jk) ) * umask(ji,jj,jk) ! i-gradient of T and S at jj 432 zdis(ji,jj,jk) = ( sb(ji+1,jj,jk) - sb(ji,jj,jk) ) * umask(ji,jj,jk) 433 zdjt(ji,jj,jk) = ( tb(ji,jj+1,jk) - tb(ji,jj,jk) ) * vmask(ji,jj,jk) ! j-gradient of T and S at jj 434 zdjs(ji,jj,jk) = ( sb(ji,jj+1,jk) - sb(ji,jj,jk) ) * vmask(ji,jj,jk) 435 END DO 436 END DO 437 END DO 438 IF( ln_zps ) THEN ! partial steps: correction at the last level 427 CALL eos_alpbet( tsb, zalbet, zbeta0 ) !== before local thermal/haline expension ratio at T-points ==! 428 ! 429 DO jl = 0, 1 !== unmasked before density i- j-, k-gradients ==! 430 ! 431 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln) 432 DO jk = 1, jpkm1 ! done each pair of triad 433 DO jj = 1, jpjm1 ! NB: not masked ==> a minimum value is set 434 DO ji = 1, fs_jpim1 ! vector opt. 435 zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! i-gradient of T & S at u-point 436 zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 437 zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! j-gradient of T & S at v-point 438 zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 439 zdxrho_raw = ( - zalbet(ji+ip,jj ,jk) * zdit + zbeta0*zdis ) / e1u(ji,jj) 440 zdyrho_raw = ( - zalbet(ji ,jj+jp,jk) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 441 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 442 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 443 END DO 444 END DO 445 END DO 446 ! 447 IF( ln_zps.and.l_grad_zps ) THEN ! partial steps: correction of i- & j-grad on bottom 439 448 # if defined key_vectopt_loop 440 DO jj = 1, 1441 DO ji = 1, jpij-jpi! vector opt. (forced unrolling)449 DO jj = 1, 1 450 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 442 451 # else 443 DO jj = 1, jpjm1444 DO ji = 1, jpim1452 DO jj = 1, jpjm1 453 DO ji = 1, jpim1 445 454 # endif 446 zdit(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_tem) ! i-gradient of T and S 447 zdis(ji,jj,mbku(ji,jj)) = gtsu(ji,jj,jp_sal) 448 zdjt(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_tem) ! j-gradient of T and S 449 zdjs(ji,jj,mbkv(ji,jj)) = gtsv(ji,jj,jp_sal) 450 END DO 451 END DO 452 ENDIF 453 ! 454 zdkt(:,:,1) = 0._wp !== before vertical T & S gradient at w-level ==! 455 zdks(:,:,1) = 0._wp 456 DO jk = 2, jpk 457 zdkt(:,:,jk) = ( tb(:,:,jk-1) - tb(:,:,jk) ) * tmask(:,:,jk) 458 zdks(:,:,jk) = ( sb(:,:,jk-1) - sb(:,:,jk) ) * tmask(:,:,jk) 459 END DO 460 ! 461 ! 462 DO jl = 0, 1 !== density i-, j-, and k-gradients ==! 463 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln) 464 DO jk = 1, jpkm1 ! done each pair of triad 465 DO jj = 1, jpjm1 ! NB: not masked due to the minimum value set 466 DO ji = 1, fs_jpim1 ! vector opt. 467 zdxrho_raw = ( zalpha(ji+ip,jj ,jk) * zdit(ji,jj,jk) + zbeta(ji+ip,jj ,jk) * zdis(ji,jj,jk) ) / e1u(ji,jj) 468 zdyrho_raw = ( zalpha(ji ,jj+jp,jk) * zdjt(ji,jj,jk) + zbeta(ji ,jj+jp,jk) * zdjs(ji,jj,jk) ) / e2v(ji,jj) 469 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 470 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 455 iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points) 456 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 457 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 458 zdxrho_raw = ( - zalbet(ji+ip,jj ,iku) * zdit + zbeta0*zdis ) / e1u(ji,jj) 459 zdyrho_raw = ( - zalbet(ji ,jj+jp,ikv) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 460 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 461 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 471 462 END DO 472 463 END DO 473 END DO 474 END DO 475 DO kp = 0, 1 !== density i-, j-, and k-gradients ==! 476 DO jk = 1, jpkm1 ! done each pair of triad 477 DO jj = 1, jpj ! NB: not masked due to the minimum value set 478 DO ji = 1, jpi ! vector opt. 479 zdzrho_raw = ( zalpha(ji,jj,jk) * zdkt(ji,jj,jk+kp) + zbeta(ji,jj,jk) * zdks(ji,jj,jk+kp) ) & 480 & / fse3w(ji,jj,jk+kp) 481 zdzrho(ji ,jj ,jk, kp) = - MIN( - repsln, zdzrho_raw ) ! force zdzrho >= repsln 482 END DO 483 END DO 484 END DO 485 END DO 486 ! 487 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==! 464 ENDIF 465 ! 466 END DO 467 468 DO kp = 0, 1 !== unmasked before density i- j-, k-gradients ==! 469 DO jk = 1, jpkm1 ! done each pair of triad 470 DO jj = 1, jpj ! NB: not masked ==> a minimum value is set 471 DO ji = 1, jpi ! vector opt. 472 IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp 473 zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) ) 474 zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) ) 475 ELSE 476 zdkt = 0._wp ! 1st level gradient set to zero 477 zdks = 0._wp 478 ENDIF 479 zdzrho_raw = ( - zalbet(ji ,jj ,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp) 480 zdzrho(ji ,jj ,jk, kp) = - MIN( - repsln, zdzrho_raw ) ! force zdzrho >= repsln 481 END DO 482 END DO 483 END DO 484 END DO 485 ! 486 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==! 488 487 DO ji = 1, jpi 489 488 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth … … 492 491 END DO 493 492 ! 494 ! !== intialisations to zero ==!495 ! 496 wslp2 (:,:,:) = 0._wp 497 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp 493 ! !== intialisations to zero ==! 494 ! 495 wslp2 (:,:,:) = 0._wp ! wslp2 will be cumulated 3D field set to zero 496 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 498 497 triadj_g(:,:,1,:,:) = 0._wp ; triadj_g(:,:,jpk,:,:) = 0._wp 499 !!gm _iso set to zero missing500 triadi (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp! set surface and bottom slope to zero501 triadj (:,:,1,:,:) = 0._wp ; triadj(:,:,jpk,:,:) = 0._wp502 498 !!gm _iso set to zero missing 499 triadi (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 500 triadj (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp 501 503 502 !-------------------------------------! 504 503 ! Triads just below the Mixed Layer ! 505 504 !-------------------------------------! 506 505 ! 507 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base508 DO kp = 0, 1 ! with only the slope-max limit and MASKED506 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 507 DO kp = 0, 1 ! with only the slope-max limit and MASKED 509 508 DO jj = 1, jpjm1 510 509 DO ji = 1, fs_jpim1 511 510 ip = jl ; jp = jl 512 511 jk = MIN( nmln(ji+ip,jj) , mbkt(ji+ip,jj) ) + 1 ! ML level+1 (MIN in case ML depth is the ocean depth) 512 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 513 513 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 514 & +( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj) ) * umask(ji,jj,jk)514 & - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj) ) * umask(ji,jj,jk) 515 515 jk = MIN( nmln(ji,jj+jp) , mbkt(ji,jj+jp) ) + 1 516 516 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 517 & +( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk)517 & - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 518 518 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 519 519 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) … … 527 527 !-------------------------------------! 528 528 ! 529 DO kp = 0, 1 ! k-index of triads529 DO kp = 0, 1 ! k-index of triads 530 530 DO jl = 0, 1 531 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes)531 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes) 532 532 DO jk = 1, jpkm1 533 533 DO jj = 1, jpjm1 534 DO ji = 1, fs_jpim1 ! vector opt.534 DO ji = 1, fs_jpim1 ! vector opt. 535 535 ! 536 536 ! Calculate slope relative to geopotentials used for GM skew fluxes 537 ! For s-coordinate, subtract slope at t-points (equivalent to *adding* gradient of depth)537 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 538 538 ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 539 539 ! masked by umask taken at the level of dz(rho) … … 544 544 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 545 545 zti_coord = ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 546 ztj_coord = ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) 547 ! unmasked 548 zti_g_raw = zti_raw + zti_coord ! ref to geopot surfaces 549 ztj_g_raw = ztj_raw + ztj_coord 546 ztj_coord = ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) ! unmasked 547 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 548 ztj_g_raw = ztj_raw - ztj_coord 550 549 zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 551 550 ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) … … 562 561 ztj_g_lim = zfactj * ztj_g_lim & 563 562 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 564 & * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ! masked565 ! 566 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim * umask(ji,jj,jk+kp) 563 & * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) 564 ! 565 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim * umask(ji,jj,jk+kp) ! masked 567 566 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim * vmask(ji,jj,jk+kp) 568 567 ! … … 574 573 ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 575 574 ! 576 zti_lim = zti_g_lim - zti_coord ! remove the coordinate slope ==> relative to coordinate surfaces 577 ztj_lim = ztj_g_lim - ztj_coord 578 zti_lim2 = zti_lim * zti_lim * umask(ji,jj,jk+kp) ! square of limited slopes ! masked <<== 579 ztj_lim2 = ztj_lim * ztj_lim * vmask(ji,jj,jk+kp) 580 ! 575 zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces 576 ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 577 zti_lim2 = zti_lim * zti_lim ! square of limited slopes ! masked <<== 578 ztj_lim2 = ztj_lim * ztj_lim 579 ! 580 IF( ln_triad_iso ) THEN 581 zti_raw = zti_lim2 / zti_raw 582 ztj_raw = ztj_lim2 / ztj_raw 583 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 584 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 585 zti_lim = zfacti * zti_lim & 586 & + ( 1._wp - zfacti ) * zti_raw 587 ztj_lim = zfactj * ztj_lim & 588 & + ( 1._wp - zfactj ) * ztj_raw 589 ENDIF 590 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim 591 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim 592 ! 581 593 zbu = e1u(ji ,jj) * e2u(ji ,jj) * fse3u(ji ,jj,jk ) 582 594 zbv = e1v(ji ,jj) * e2v(ji ,jj) * fse3v(ji ,jj,jk ) … … 584 596 zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 585 597 ! 586 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim2 / zti_raw ! masked 587 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim2 / ztj_raw 588 ! 589 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked 590 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2 ! masked 598 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked 599 ! agn may need to change to using ztj_g_lim**2, as wslp2 just used for eddy growth rate, needs *real* slope 600 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_lim2 ! masked 591 601 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2 592 602 END DO … … 597 607 ! 598 608 wslp2(:,:,1) = 0._wp ! force the surface wslp to zero 599 609 600 610 CALL lbc_lnk( wslp2, 'W', 1. ) ! lateral boundary confition on wslp2 only ==>>> gm : necessary ? to be checked 601 611 ! … … 610 620 !! *** ROUTINE ldf_slp_mxl *** 611 621 !! 612 !! ** Purpose : Compute the slopes of iso-neutral surface just below 622 !! ** Purpose : Compute the slopes of iso-neutral surface just below 613 623 !! the mixed layer. 614 624 !! … … 619 629 !! 620 630 !! ** Action : uslpml, wslpiml : i- & j-slopes of neutral surfaces 621 !! vslpml, wslpjml just below the mixed layer 631 !! vslpml, wslpjml just below the mixed layer 622 632 !! omlmask : mixed layer mask 623 633 !!---------------------------------------------------------------------- … … 627 637 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_dzr ! z-gradient of density (T-point) 628 638 !! 629 INTEGER :: ji , jj , jk ! dummy loop indices630 INTEGER :: iku, ikv, ik, ikm1 ! local integers639 INTEGER :: ji , jj , jk ! dummy loop indices 640 INTEGER :: iku, ikv, ik, ikm1 ! local integers 631 641 REAL(wp) :: zeps, zm1_g, zm1_2g ! local scalars 632 642 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - … … 644 654 wslpjml(1,:) = 0._wp ; wslpjml(jpi,:) = 0._wp 645 655 ! 646 ! !== surface mixed layer mask !647 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise656 ! !== surface mixed layer mask ! 657 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise 648 658 # if defined key_vectopt_loop 649 659 DO jj = 1, 1 650 DO ji = 1, jpij ! vector opt. (forced unrolling)660 DO ji = 1, jpij ! vector opt. (forced unrolling) 651 661 # else 652 662 DO jj = 1, jpj … … 679 689 DO ji = 2, jpim1 680 690 # endif 681 ! !== Slope at u- & v-points just below the Mixed Layer ==!691 ! !== Slope at u- & v-points just below the Mixed Layer ==! 682 692 ! 683 ! 693 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 684 694 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 685 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 695 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 686 696 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 687 697 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 688 ! 698 ! !- horizontal density gradient at u- & v-points 689 699 zau = p_gru(ji,jj,iku) / e1u(ji,jj) 690 700 zav = p_grv(ji,jj,ikv) / e2v(ji,jj) 691 ! 692 ! 701 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 702 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 693 703 zbu = MIN( zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) ) 694 704 zbv = MIN( zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) ) 695 ! 705 ! !- Slope at u- & v-points (uslpml, vslpml) 696 706 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 697 707 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 698 708 ! 699 ! !== i- & j-slopes at w-points just below the Mixed Layer ==!709 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! 700 710 ! 701 711 ik = MIN( nmln(ji,jj) + 1, jpk ) 702 712 ikm1 = MAX( 1, ik-1 ) 703 ! 713 ! !- vertical density gradient for w-slope (from N^2) 704 714 zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 705 ! 715 ! !- horizontal density i- & j-gradient at w-points 706 716 zci = MAX( umask(ji-1,jj,ik ) + umask(ji,jj,ik ) & 707 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 717 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 708 718 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 709 719 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) … … 712 722 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & 713 723 & + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik) 714 ! 715 ! 724 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 725 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 716 726 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zai ) ) 717 727 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj ) ) 718 ! 728 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 719 729 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 720 730 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 721 731 END DO 722 732 END DO 723 !!gm this lbc_lnk should be useless....733 !!gm this lbc_lnk should be useless.... 724 734 CALL lbc_lnk( uslpml , 'U', -1. ) ; CALL lbc_lnk( vslpml , 'V', -1. ) ! lateral boundary cond. (sign change) 725 735 CALL lbc_lnk( wslpiml, 'W', -1. ) ; CALL lbc_lnk( wslpjml, 'W', -1. ) ! lateral boundary conditions … … 734 744 !! ** Purpose : Initialization for the isopycnal slopes computation 735 745 !! 736 !! ** Method : read the nammbf namelist and check the parameter 737 !! 746 !! ** Method : read the nammbf namelist and check the parameter 747 !! values called by tra_dmp at the first timestep (nit000) 738 748 !!---------------------------------------------------------------------- 739 749 INTEGER :: ji, jj, jk ! dummy loop indices 740 750 INTEGER :: ierr ! local integer 741 751 !!---------------------------------------------------------------------- 742 743 IF(lwp) THEN 752 753 IF(lwp) THEN 744 754 WRITE(numout,*) 745 755 WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing' 746 756 WRITE(numout,*) '~~~~~~~~~~~~' 747 757 ENDIF 748 758 749 759 IF( ln_traldf_grif ) THEN ! Griffies operator : triad of slopes 750 760 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 ) … … 755 765 IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 756 766 ! 757 IF( ( ln_traldf_hor .OR. ln_dynldf_hor ) .AND. ln_sco ) &758 CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator in s-coordinate not supported' )759 !760 767 ELSE ! Madec operator : slopes at u-, v-, and w-points 761 768 ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) , & … … 770 777 wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp 771 778 772 !!gm I no longer understand this.....779 !!gm I no longer understand this..... 773 780 IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 774 781 IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' … … 803 810 LOGICAL, PUBLIC, PARAMETER :: lk_ldfslp = .FALSE. !: slopes flag 804 811 CONTAINS 805 SUBROUTINE ldf_slp( kt, prd, pn2 ) 806 INTEGER, INTENT(in) :: kt 812 SUBROUTINE ldf_slp( kt, prd, pn2 ) ! Dummy routine 813 INTEGER, INTENT(in) :: kt 807 814 REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2 808 815 WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1) … … 812 819 WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt 813 820 END SUBROUTINE ldf_slp_grif 814 SUBROUTINE ldf_slp_init ! Dummy routine821 SUBROUTINE ldf_slp_init ! Dummy routine 815 822 END SUBROUTINE ldf_slp_init 816 823 #endif -
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
r2715 r2980 67 67 & ln_traldf_level, ln_traldf_hor , ln_traldf_iso, & 68 68 & ln_traldf_grif , ln_traldf_gdia, & 69 & ln_triad_iso , ln_botmix_grif, & 69 70 & rn_aht_0 , rn_ahtb_0 , rn_aeiv_0, & 70 71 & rn_slpmax … … 94 95 WRITE(numout,*) ' maximum isoppycnal slope rn_slpmax = ', rn_slpmax 95 96 WRITE(numout,*) ' + griffies operator internal controls not set via the namelist (experimental): ' 96 WRITE(numout,*) ' calculate triads twice l _triad_iso = ', l_triad_iso97 WRITE(numout,*) ' no Shapiro filter l_no_smooth = ', l_no_smooth97 WRITE(numout,*) ' calculate triads twice ln_triad_iso = ', ln_triad_iso 98 WRITE(numout,*) ' GM -->lat mixing on bottom ln_botmix_grif = ', ln_botmix_grif 98 99 WRITE(numout,*) 99 100 ENDIF -
branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90
r2715 r2980 32 32 33 33 REAL(wp), PUBLIC :: aht0, ahtb0, aeiv0 !!: OLD namelist names 34 LOGICAL , PUBLIC :: l_triad_iso = .FALSE. !: calculate triads twice 35 LOGICAL , PUBLIC :: l_no_smooth = .FALSE. !: no Shapiro smoothing 34 LOGICAL , PUBLIC :: ln_triad_iso = .FALSE. !: calculate triads twice 35 LOGICAL , PUBLIC :: ln_botmix_grif = .FALSE. !: mixing on bottom 36 LOGICAL , PUBLIC :: l_grad_zps = .FALSE. !: special treatment for Horz Tgradients w partial steps 36 37 37 38 #if defined key_traldf_c3d
Note: See TracChangeset
for help on using the changeset viewer.
