Changeset 2840
- Timestamp:
- 2011-09-15T14:38:15+02:00 (13 years ago)
- Location:
- branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO
- Files:
-
- 6 edited
-
OPA_SRC/LDF/ldfslp.F90 (modified) (32 diffs)
-
OPA_SRC/LDF/ldftra.F90 (modified) (2 diffs)
-
OPA_SRC/LDF/ldftra_oce.F90 (modified) (1 diff)
-
OPA_SRC/TRA/eosbn2.F90 (modified) (21 diffs)
-
OPA_SRC/TRA/traldf_iso_grif.F90 (modified) (14 diffs)
-
TOP_SRC/TRP/trcldf.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r2772 r2840 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.ln_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) 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 ) 462 END DO 463 END DO 464 ENDIF 465 ! 466 END DO 467 468 DO kp = 0, 1 !== unmasked before density i- j-, k-gradients ==! 464 469 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 ) 471 END DO 472 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 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 483 482 END DO 484 483 END DO … … 494 493 ! !== intialisations to zero ==! 495 494 ! 496 wslp2 (:,:,:) = 0._wp 497 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp 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 498 !!gm _iso set to zero missing 500 triadi (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp! set surface and bottom slope to zero501 triadj (:,:,1,:,:) = 0._wp ; triadj(:,:,jpk,:,:) = 0._wp502 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 ! … … 506 505 ! 507 506 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 508 DO kp = 0, 1 ! with only the slope-max limit and MASKED 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 ) … … 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 ! masked587 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim2 / ztj_raw588 !589 598 !!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 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 !!---------------------------------------------------------------------- … … 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) ) … … 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) … … 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' )767 ! IF( ( ln_traldf_hor .OR. ln_dynldf_hor ) .AND. ln_sco ) & 768 ! CALL ctl_stop( 'ldf_slp_init: horizontal Griffies operator in s-coordinate not supported' ) 759 769 ! 760 770 ELSE ! Madec operator : slopes at u-, v-, and w-points … … 804 814 CONTAINS 805 815 SUBROUTINE ldf_slp( kt, prd, pn2 ) ! Dummy routine 806 INTEGER, INTENT(in) :: kt 816 INTEGER, INTENT(in) :: kt 807 817 REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2 808 818 WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1) -
branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra.F90
r2715 r2840 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, ln_grad_zps, & 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_iso 97 WRITE(numout,*) ' no Shapiro filter l_no_smooth = ', l_no_smooth 97 WRITE(numout,*) ' calculate triads twice ln_triad_iso = ', ln_triad_iso 98 WRITE(numout,*) ' special Tgradts w ptl steps ln_grad_zps = ', ln_grad_zps 99 WRITE(numout,*) ' GM -->lat mixing on bottom ln_botmix = ', ln_botmix 98 100 WRITE(numout,*) 99 101 ENDIF -
branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_oce.F90
r2715 r2840 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_grad_zps = .FALSE. !: special treatment for Horz Tgradients w partial steps 36 LOGICAL , PUBLIC :: ln_botmix = .FALSE. !: mixing on bottom 36 37 37 38 #if defined key_traldf_c3d -
branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
r2715 r2840 3 3 !! *** MODULE eosbn2 *** 4 4 !! Ocean diagnostic variable : equation of state - in situ and potential density 5 !! - Brunt-Vaisala frequency 5 !! - Brunt-Vaisala frequency 6 6 !!============================================================================== 7 7 !! History : OPA ! 1989-03 (O. Marti) Original code … … 27 27 !! eos_insitu_2d : Compute the in situ density for 2d fields 28 28 !! eos_bn2 : Compute the Brunt-Vaisala frequency 29 !! eos_alpbet : calculates the in situ thermal and haline expansion coeff.29 !! eos_alpbet : calculates the in situ thermal/haline expansion ratio 30 30 !! tfreez : Compute the surface freezing temperature 31 31 !! eos_init : set eos parameters (namelist) … … 41 41 PRIVATE 42 42 43 ! !! * Interface 43 ! !! * Interface 44 44 INTERFACE eos 45 45 MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d 46 END INTERFACE 46 END INTERFACE 47 47 INTERFACE bn2 48 48 MODULE PROCEDURE eos_bn2 49 END INTERFACE 49 END INTERFACE 50 50 51 51 PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules … … 61 61 62 62 REAL(wp), PUBLIC :: ralpbet !: alpha / beta ratio 63 63 64 64 !! * Substitutions 65 65 # include "domzgr_substitute.h90" … … 75 75 !!---------------------------------------------------------------------- 76 76 !! *** ROUTINE eos_insitu *** 77 !! 78 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 77 !! 78 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 79 79 !! potential temperature and salinity using an equation of state 80 80 !! defined through the namelist parameter nn_eos. … … 134 134 !CDIR NOVERRCHK 135 135 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 136 ! 136 ! 137 137 DO jk = 1, jpkm1 138 138 DO jj = 1, jpj … … 199 199 !!---------------------------------------------------------------------- 200 200 !! *** ROUTINE eos_insitu_pot *** 201 !! 201 !! 202 202 !! ** Purpose : Compute the in situ density (ratio rho/rau0) and the 203 203 !! potential volumic mass (Kg/m3) from potential temperature and 204 !! salinity fields using an equation of state defined through the 204 !! salinity fields using an equation of state defined through the 205 205 !! namelist parameter nn_eos. 206 206 !! … … 230 230 !! nn_eos = 2 : linear equation of state function of temperature and 231 231 !! salinity 232 !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0 232 !! prd(t,s) = ( rho(t,s) - rau0 ) / rau0 233 233 !! = rn_beta * s - rn_alpha * tn - 1. 234 234 !! rhop(t,s) = rho(t,s) … … 265 265 !CDIR NOVERRCHK 266 266 zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) ) 267 ! 267 ! 268 268 DO jk = 1, jpkm1 269 269 DO jj = 1, jpj … … 336 336 !! *** ROUTINE eos_insitu_2d *** 337 337 !! 338 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 338 !! ** Purpose : Compute the in situ density (ratio rho/rau0) from 339 339 !! potential temperature and salinity using an equation of state 340 340 !! defined through the namelist parameter nn_eos. * 2D field case … … 374 374 ! ! 2 : salinity [psu] 375 375 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m] 376 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 376 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density 377 377 !! 378 378 INTEGER :: ji, jj ! dummy loop indices … … 449 449 DO jj = 1, jpjm1 450 450 DO ji = 1, fs_jpim1 ! vector opt. 451 prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 451 prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1) 452 452 END DO 453 453 END DO … … 468 468 !! ** Purpose : Compute the local Brunt-Vaisala frequency at the time- 469 469 !! step of the input arguments 470 !! 470 !! 471 471 !! ** Method : 472 472 !! * nn_eos = 0 : UNESCO sea water properties … … 482 482 !! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 483 483 !! The use of potential density to compute N^2 introduces e r r o r 484 !! in the sign of N^2 at great depths. We recommand the use of 484 !! in the sign of N^2 at great depths. We recommand the use of 485 485 !! nn_eos = 0, except for academical studies. 486 486 !! Macro-tasked on horizontal slab (jk-loop) … … 497 497 !! 498 498 INTEGER :: ji, jj, jk ! dummy loop indices 499 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars 499 REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars 500 500 #if defined key_zdfddm 501 501 REAL(wp) :: zds ! local scalars … … 504 504 505 505 ! pn2 : interior points only (2=< jk =< jpkm1 ) 506 ! -------------------------- 506 ! -------------------------- 507 507 ! 508 508 SELECT CASE( nn_eos ) … … 542 542 & - 0.121555e-07_wp ) * zh 543 543 ! 544 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 544 pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! N^2 545 545 & * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) & 546 546 & - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) … … 565 565 & - rn_beta * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) ) ) & 566 566 & / fse3w(:,:,jk) * tmask(:,:,jk) 567 END DO 567 END DO 568 568 #if defined key_zdfddm 569 569 DO jk = 2, jpkm1 ! Rrau = (alpha / beta) (dk[t] / dk[s]) 570 570 DO jj = 1, jpj 571 571 DO ji = 1, jpi 572 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 572 zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) 573 573 IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp 574 574 rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds … … 587 587 588 588 589 SUBROUTINE eos_alpbet( pts, palp h, pbeta)590 !!---------------------------------------------------------------------- 591 !! *** ROUTINE ldf_slp_grif***592 !! 593 !! ** Purpose : Calculates the thermal and haline expansion coefficients at T-points.594 !! 595 !! ** Method : calculates alpha and beta at T-points589 SUBROUTINE eos_alpbet( pts, palpbet, beta0 ) 590 !!---------------------------------------------------------------------- 591 !! *** ROUTINE eos_alpbet *** 592 !! 593 !! ** Purpose : Calculates the in situ thermal/haline expansion ratio at T-points 594 !! 595 !! ** Method : calculates alpha / beta ratio at T-points 596 596 !! * nn_eos = 0 : UNESCO sea water properties 597 !! The brunt-vaisala frequency is computed using the polynomial 598 !! polynomial expression of McDougall (1987): 599 !! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w 600 !! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is 601 !! computed and used in zdfddm module : 602 !! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] ) 597 !! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial 598 !! polynomial expression of McDougall (1987). 599 !! Scalar beta0 is returned = 1. 603 600 !! * nn_eos = 1 : linear equation of state (temperature only) 604 !! N^2 = grav * rn_alpha * dk[ t ]/e3w 601 !! The ratio is undefined, so we return alpha as palpbet 602 !! Scalar beta0 is returned = 0. 605 603 !! * nn_eos = 2 : linear equation of state (temperature & salinity) 606 !! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w 607 !! * nn_eos = 3 : Jackett JAOT 2003 ??? 608 !! 609 !! ** Action : - palph, pbeta : thermal and haline expansion coeff. at T-point 610 !!---------------------------------------------------------------------- 611 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 612 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palph, pbeta ! thermal & haline expansion coeff. 613 ! 604 !! The alpha/beta ratio is returned as ralpbet 605 !! Scalar beta0 is returned = 1. 606 !! 607 !! ** Action : - palpbet : thermal/haline expansion ratio at T-points 608 !! : beta0 : 1. or 0. 609 !!---------------------------------------------------------------------- 610 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity 611 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palpbet ! thermal/haline expansion ratio 612 REAL(wp), INTENT( out) :: beta0 ! set = 1 except with case 1 eos, rho=rho(T) 613 !! 614 614 INTEGER :: ji, jj, jk ! dummy loop indices 615 REAL(wp) :: zt, zs, zh ! local scalars 615 REAL(wp) :: zt, zs, zh ! local scalars 616 616 !!---------------------------------------------------------------------- 617 617 ! … … 624 624 zt = pts(ji,jj,jk,jp_tem) ! potential temperature 625 625 zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35) 626 zh = fsdept(ji,jj,jk) ! depth in meters 627 ! 628 pbeta(ji,jj,jk) = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt & 629 & - 0.301985e-05_wp ) * zt & 630 & + 0.785567e-03_wp & 631 & + ( 0.515032e-08_wp * zs & 632 & + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs & 633 & + ( ( 0.121551e-17_wp * zh & 634 & - 0.602281e-15_wp * zs & 635 & - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh & 636 & + 0.408195e-10_wp * zs & 637 & + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt & 638 & - 0.121555e-07_wp ) * zh 639 ! 640 palph(ji,jj,jk) = - pbeta(ji,jj,jk) * & 641 & ((( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & 642 & - 0.203814e-03_wp ) * zt & 643 & + 0.170907e-01_wp ) * zt & 644 & + 0.665157e-01_wp & 645 & + ( - 0.678662e-05_wp * zs & 646 & - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs & 647 & + ( ( - 0.302285e-13_wp * zh & 648 & - 0.251520e-11_wp * zs & 649 & + 0.512857e-12_wp * zt * zt ) * zh & 650 & - 0.164759e-06_wp * zs & 651 & +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt & 652 & + 0.380374e-04_wp ) * zh) 626 zh = fsdept(ji,jj,jk) ! depth in meters 627 ! 628 palpbet(ji,jj,jk) = & 629 & ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & 630 & - 0.203814e-03_wp ) * zt & 631 & + 0.170907e-01_wp ) * zt & 632 & + 0.665157e-01_wp & 633 & + ( - 0.678662e-05_wp * zs & 634 & - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs & 635 & + ( ( - 0.302285e-13_wp * zh & 636 & - 0.251520e-11_wp * zs & 637 & + 0.512857e-12_wp * zt * zt ) * zh & 638 & - 0.164759e-06_wp * zs & 639 & +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt & 640 & + 0.380374e-04_wp ) * zh 653 641 END DO 654 642 END DO 655 643 END DO 656 ! 657 CASE ( 1 ) 658 palph(:,:,:) = - rn_alpha 659 pbeta(:,:,:) = 0._wp 660 ! 661 CASE ( 2 ) 662 palph(:,:,:) = - rn_alpha 663 pbeta(:,:,:) = rn_beta 644 beta0 = 1._wp 645 ! 646 CASE ( 1 ) !== Linear formulation = F( temperature ) ==! 647 palpbet(:,:,:) = rn_alpha 648 beta0 = 0._wp 649 ! 650 CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==! 651 palpbet(:,:,:) = ralpbet 652 beta0 = 1._wp 664 653 ! 665 654 CASE DEFAULT … … 747 736 748 737 !!====================================================================== 749 END MODULE eosbn2 738 END MODULE eosbn2 -
branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_grif.F90
r2715 r2840 4 4 !! Ocean tracers: horizontal component of the lateral tracer mixing trend 5 5 !!====================================================================== 6 !! History : 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) 6 !! History : 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) 7 7 !! ! Griffies operator version adapted from traldf_iso.F90 8 8 !!---------------------------------------------------------------------- … … 11 11 !! 'key_ldfslp' slope of the lateral diffusive direction 12 12 !!---------------------------------------------------------------------- 13 !! tra_ldf_iso_grif : update the tracer trend with the horizontal component 14 !! of the Griffies iso-neutral laplacian operator 13 !! tra_ldf_iso_grif : update the tracer trend with the horizontal component 14 !! of the Griffies iso-neutral laplacian operator 15 15 !!---------------------------------------------------------------------- 16 16 USE oce ! ocean dynamics and active tracers … … 53 53 !! *** ROUTINE tra_ldf_iso_grif *** 54 54 !! 55 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 56 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 55 !! ** Purpose : Compute the before horizontal tracer (t & s) diffusive 56 !! trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and 57 57 !! add it to the general trend of tracer equation. 58 58 !! 59 !! ** Method : The horizontal component of the lateral diffusive trends 59 !! ** Method : The horizontal component of the lateral diffusive trends 60 60 !! is provided by a 2nd order operator rotated along neural or geopo- 61 61 !! tential surfaces to which an eddy induced advection can be added … … 67 67 !! 68 68 !! 2nd part : horizontal fluxes of the lateral mixing operator 69 !! ======== 69 !! ======== 70 70 !! zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ] 71 71 !! - aht e2u*uslp dk[ mi(mk(tb)) ] … … 99 99 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu, pgv ! tracer gradient at pstep levels 100 100 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 101 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 102 102 REAL(wp) , INTENT(in ) :: pahtb0 ! background diffusion coef 103 103 ! … … 143 143 144 144 !!---------------------------------------------------------------------- 145 !! 0 - calculate ah_wslp2, psix_eiv, psiy_eiv 145 !! 0 - calculate ah_wslp2, psix_eiv, psiy_eiv 146 146 !!---------------------------------------------------------------------- 147 147 … … 159 159 DO jj = 1, jpjm1 160 160 DO ji = 1, fs_jpim1 161 ze1ur = 1._wp / e1u(ji,jj) 161 162 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 162 163 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 163 164 zah = fsahtu(ji,jj,jk) ! fsaht(ji+ip,jj,jk) 164 165 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 165 zslope2 = zslope_skew - ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+kp) 166 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 167 ! (do this by *adding* gradient of depth) 168 zslope2 = zslope_skew + ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * ze1ur * umask(ji,jj,jk+kp) 166 169 zslope2 = zslope2 *zslope2 167 170 ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) & … … 182 185 DO jj = 1, jpjm1 183 186 DO ji=1,fs_jpim1 187 ze2vr = 1._wp / e2v(ji,jj) 184 188 ze3wr = 1.0_wp / fse3w(ji,jj+jp,jk+kp) 185 189 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 186 zah = fsaht u(ji,jj,jk) !fsaht(ji,jj+jp,jk)190 zah = fsahtv(ji,jj,jk) !fsaht(ji,jj+jp,jk) 187 191 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 188 zslope2 = zslope_skew - ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+kp) 192 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 193 ! (do this by *adding* gradient of depth) 194 zslope2 = zslope_skew + ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) * ze2vr * vmask(ji,jj,jk+kp) 189 195 zslope2 = zslope2 * zslope2 190 196 ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) & … … 207 213 zftu(:,:,:) = 0._wp 208 214 zftv(:,:,:) = 0._wp 209 ! 215 ! 210 216 DO jk = 1, jpkm1 !== before lateral T & S gradients at T-level jk ==! 211 217 DO jj = 1, jpjm1 … … 216 222 END DO 217 223 END DO 218 IF( ln_zps ) THEN! partial steps: correction at the last level224 IF( ln_zps.and.ln_grad_zps ) THEN ! partial steps: correction at the last level 219 225 # if defined key_vectopt_loop 220 226 DO jj = 1, 1 … … 224 230 DO ji = 1, jpim1 225 231 # endif 226 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 227 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 232 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 233 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 228 234 END DO 229 235 END DO … … 244 250 ENDIF 245 251 246 IF( .NOT. l_triad_iso ) THEN 247 triadi = triadi_g 248 triadj = triadj_g 249 ENDIF 250 251 DO ip = 0, 1 !== Horizontal & vertical fluxes 252 DO kp = 0, 1 253 DO jj = 1, jpjm1 254 DO ji = 1, fs_jpim1 255 ze1ur = 1._wp / e1u(ji,jj) 256 zdxt = zdit(ji,jj,jk) * ze1ur 257 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 258 zdzt = zdkt(ji+ip,jj,kp) * ze3wr 259 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 260 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 261 262 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 263 zah = fsahtu(ji,jj,jk) !*umask(ji,jj,jk+kp) !fsaht(ji+ip,jj,jk) ===>> ???? 264 zah_slp = zah * zslope_iso 265 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew !fsaeit(ji+ip,jj,jk)*zslope_skew 266 zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 267 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 252 253 IF (ln_botmix) THEN 254 DO ip = 0, 1 !== Horizontal & vertical fluxes 255 DO kp = 0, 1 256 DO jj = 1, jpjm1 257 DO ji = 1, fs_jpim1 258 ze1ur = 1._wp / e1u(ji,jj) 259 zdxt = zdit(ji,jj,jk) * ze1ur 260 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 261 zdzt = zdkt(ji+ip,jj,kp) * ze3wr 262 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 263 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 264 265 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 266 ! ln_botmix is .T. don't mask zah for bottom half cells 267 zah = fsahtu(ji,jj,jk) !*umask(ji,jj,jk+kp) !fsaht(ji+ip,jj,jk) ===>> ???? 268 zah_slp = zah * zslope_iso 269 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew !fsaeit(ji+ip,jj,jk)*zslope_skew 270 zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 271 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 272 END DO 268 273 END DO 269 274 END DO 270 275 END DO 271 END DO 272 273 DO jp = 0, 1 274 DO kp = 0, 1 275 DO jj = 1, jpjm1 276 DO ji = 1, fs_jpim1 277 ze2vr = 1._wp / e2v(ji,jj) 278 zdyt = zdjt(ji,jj,jk) * ze2vr 279 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 280 zdzt = zdkt(ji,jj+jp,kp) * ze3wr 281 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 282 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 283 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 284 zah = fsahtv(ji,jj,jk) !*vmask(ji,jj,jk+kp) !fsaht(ji,jj+jp,jk) 285 zah_slp = zah * zslope_iso 286 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew !fsaeit(ji,jj+jp,jk)*zslope_skew 287 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 288 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 276 277 DO jp = 0, 1 278 DO kp = 0, 1 279 DO jj = 1, jpjm1 280 DO ji = 1, fs_jpim1 281 ze2vr = 1._wp / e2v(ji,jj) 282 zdyt = zdjt(ji,jj,jk) * ze2vr 283 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 284 zdzt = zdkt(ji,jj+jp,kp) * ze3wr 285 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 286 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 287 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 288 ! ln_botmix is .T. don't mask zah for bottom half cells 289 zah = fsahtv(ji,jj,jk) !*vmask(ji,jj,jk+kp) !fsaht(ji,jj+jp,jk) 290 zah_slp = zah * zslope_iso 291 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew !fsaeit(ji,jj+jp,jk)*zslope_skew 292 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 293 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 294 END DO 289 295 END DO 290 296 END DO 291 297 END DO 292 END DO 293 298 ELSE 299 DO ip = 0, 1 !== Horizontal & vertical fluxes 300 DO kp = 0, 1 301 DO jj = 1, jpjm1 302 DO ji = 1, fs_jpim1 303 ze1ur = 1._wp / e1u(ji,jj) 304 zdxt = zdit(ji,jj,jk) * ze1ur 305 ze3wr = 1._wp / fse3w(ji+ip,jj,jk+kp) 306 zdzt = zdkt(ji+ip,jj,kp) * ze3wr 307 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 308 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 309 310 zbu = 0.25_wp * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 311 ! ln_botmix is .F. mask zah for bottom half cells 312 zah = fsahtu(ji,jj,jk) * umask(ji,jj,jk+kp) !fsaht(ji+ip,jj,jk) ===>> ???? 313 zah_slp = zah * zslope_iso 314 zaei_slp = fsaeiw(ji+ip,jj,jk) * zslope_skew !fsaeit(ji+ip,jj,jk)*zslope_skew 315 zftu(ji,jj,jk) = zftu(ji,jj,jk) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 316 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 317 END DO 318 END DO 319 END DO 320 END DO 321 322 DO jp = 0, 1 323 DO kp = 0, 1 324 DO jj = 1, jpjm1 325 DO ji = 1, fs_jpim1 326 ze2vr = 1._wp / e2v(ji,jj) 327 zdyt = zdjt(ji,jj,jk) * ze2vr 328 ze3wr = 1._wp / fse3w(ji,jj+jp,jk+kp) 329 zdzt = zdkt(ji,jj+jp,kp) * ze3wr 330 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 331 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 332 zbv = 0.25_wp * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 333 ! ln_botmix is .F. mask zah for bottom half cells 334 zah = fsahtv(ji,jj,jk) * vmask(ji,jj,jk+kp) !fsaht(ji,jj+jp,jk) 335 zah_slp = zah * zslope_iso 336 zaei_slp = fsaeiw(ji,jj+jp,jk) * zslope_skew !fsaeit(ji,jj+jp,jk)*zslope_skew 337 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 338 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 339 END DO 340 END DO 341 END DO 342 END DO 343 END IF 294 344 ! !== divergence and add to the general trend ==! 295 345 DO jj = 2 , jpjm1 … … 320 370 #if defined key_diaar5 321 371 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 322 z2d(:,:) = 0._wp 323 zztmp = rau0 * rcp 372 z2d(:,:) = 0._wp 373 zztmp = rau0 * rcp 324 374 DO jk = 1, jpkm1 325 375 DO jj = 2, jpjm1 326 376 DO ji = fs_2, fs_jpim1 ! vector opt. 327 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 377 z2d(ji,jj) = z2d(ji,jj) + zftu(ji,jj,jk) 328 378 END DO 329 379 END DO … … 332 382 CALL lbc_lnk( z2d, 'U', -1. ) 333 383 CALL iom_put( "udiff_heattr", z2d ) ! heat transport in i-direction 334 z2d(:,:) = 0._wp 384 z2d(:,:) = 0._wp 335 385 DO jk = 1, jpkm1 336 386 DO jj = 2, jpjm1 337 387 DO ji = fs_2, fs_jpim1 ! vector opt. 338 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 388 z2d(ji,jj) = z2d(ji,jj) + zftv(ji,jj,jk) 339 389 END DO 340 390 END DO -
branches/2011/dev_r2782_NOCS_Griffies/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90
r2715 r2840 76 76 WRITE(charout, FMT="('ldf0 ')") ; CALL prt_ctl_trc_info(charout) 77 77 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) 78 CALL tra_ldf_iso ( kt, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 ) 78 IF( ln_traldf_grif ) THEN 79 CALL tra_ldf_iso_grif( kt, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 ) 80 ELSE 81 CALL tra_ldf_iso ( kt, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 ) 82 ENDIF 79 83 WRITE(charout, FMT="('ldf1 ')") ; CALL prt_ctl_trc_info(charout) 80 84 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
Note: See TracChangeset
for help on using the changeset viewer.
