- Timestamp:
- 2011-11-14T18:39:45+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_NOC_UKMO_MERGE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r2772 r3101 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 ! local integer 411 REAL(wp) :: zfacti, zfactj, zatempw,zatempu,zatempv ! local scalars 412 REAL(wp) :: zbu, zbv, zbti, zbtj ! - - 413 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_lim2, zti_g_raw, zti_g_lim 414 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_lim2, ztj_g_raw, ztj_g_lim 410 INTEGER :: iku, ikv ! local integer 411 REAL(wp) :: zfacti, zfactj ! local scalars 412 REAL(wp) :: znot_thru_surface ! local scalars 413 REAL(wp) :: zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj 414 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_g_raw, zti_g_lim 415 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim 415 416 REAL(wp) :: zdzrho_raw 417 REAL(wp) :: zbeta0 416 418 !!---------------------------------------------------------------------- 417 419 … … 424 426 !--------------------------------! 425 427 ! 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 428 CALL eos_alpbet( tsb, zalbet, zbeta0 ) !== before local thermal/haline expension ratio at T-points ==! 429 ! 430 DO jl = 0, 1 !== unmasked before density i- j-, k-gradients ==! 431 ! 432 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln) 433 DO jk = 1, jpkm1 ! done each pair of triad 434 DO jj = 1, jpjm1 ! NB: not masked ==> a minimum value is set 435 DO ji = 1, fs_jpim1 ! vector opt. 436 zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! i-gradient of T & S at u-point 437 zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 438 zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! j-gradient of T & S at v-point 439 zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 440 zdxrho_raw = ( - zalbet(ji+ip,jj ,jk) * zdit + zbeta0*zdis ) / e1u(ji,jj) 441 zdyrho_raw = ( - zalbet(ji ,jj+jp,jk) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 442 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 443 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 444 END DO 445 END DO 446 END DO 447 ! 448 IF( ln_zps.and.l_grad_zps ) THEN ! partial steps: correction of i- & j-grad on bottom 439 449 # if defined key_vectopt_loop 440 DO jj = 1, 1441 DO ji = 1, jpij-jpi! vector opt. (forced unrolling)450 DO jj = 1, 1 451 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 442 452 # else 443 DO jj = 1, jpjm1444 DO ji = 1, jpim1453 DO jj = 1, jpjm1 454 DO ji = 1, jpim1 445 455 # 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 ) 456 iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points) 457 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 458 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 459 zdxrho_raw = ( - zalbet(ji+ip,jj ,iku) * zdit + zbeta0*zdis ) / e1u(ji,jj) 460 zdyrho_raw = ( - zalbet(ji ,jj+jp,ikv) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 461 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 462 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 471 463 END DO 472 464 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 ==! 465 ENDIF 466 ! 467 END DO 468 469 DO kp = 0, 1 !== unmasked before density i- j-, k-gradients ==! 470 DO jk = 1, jpkm1 ! done each pair of triad 471 DO jj = 1, jpj ! NB: not masked ==> a minimum value is set 472 DO ji = 1, jpi ! vector opt. 473 IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp 474 zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) ) 475 zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) ) 476 ELSE 477 zdkt = 0._wp ! 1st level gradient set to zero 478 zdks = 0._wp 479 ENDIF 480 zdzrho_raw = ( - zalbet(ji ,jj ,jk) * zdkt + zbeta0*zdks ) / 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 ==! 488 488 DO ji = 1, jpi 489 489 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth … … 492 492 END DO 493 493 ! 494 ! !== intialisations to zero ==!495 ! 496 wslp2 (:,:,:) = 0._wp 497 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp 494 ! !== intialisations to zero ==! 495 ! 496 wslp2 (:,:,:) = 0._wp ! wslp2 will be cumulated 3D field set to zero 497 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 498 498 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 499 !!gm _iso set to zero missing 500 triadi (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 501 triadj (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp 502 503 503 !-------------------------------------! 504 504 ! Triads just below the Mixed Layer ! 505 505 !-------------------------------------! 506 506 ! 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 MASKED507 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 509 509 DO jj = 1, jpjm1 510 510 DO ji = 1, fs_jpim1 511 511 ip = jl ; jp = jl 512 512 jk = MIN( nmln(ji+ip,jj) , mbkt(ji+ip,jj) ) + 1 ! ML level+1 (MIN in case ML depth is the ocean depth) 513 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 513 514 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)515 & - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj) ) * umask(ji,jj,jk) 515 516 jk = MIN( nmln(ji,jj+jp) , mbkt(ji,jj+jp) ) + 1 516 517 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)518 & - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 518 519 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 519 520 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) … … 527 528 !-------------------------------------! 528 529 ! 529 DO kp = 0, 1 ! k-index of triads530 DO kp = 0, 1 ! k-index of triads 530 531 DO jl = 0, 1 531 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes)532 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes) 532 533 DO jk = 1, jpkm1 534 ! Must mask contribution to slope from dz/dx at constant s for triads jk=1,kp=0 that poke up though ocean surface 535 znot_thru_surface = REAL( 1-1/(jk+kp), wp ) !jk+kp=1,=0.; otherwise=1.0 533 536 DO jj = 1, jpjm1 534 DO ji = 1, fs_jpim1 ! vector opt.537 DO ji = 1, fs_jpim1 ! vector opt. 535 538 ! 536 539 ! 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)540 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 538 541 ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 539 542 ! masked by umask taken at the level of dz(rho) … … 543 546 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 544 547 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 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 548 549 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 550 zti_coord = znot_thru_surface * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 551 ztj_coord = znot_thru_surface * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) ! unmasked 552 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 553 ztj_g_raw = ztj_raw - ztj_coord 550 554 zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 551 555 ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 552 556 ! 553 ! Below ML use limited zti_g as is 554 ! Inside ML replace by linearly reducing sx_mlb towards surface 557 ! Below ML use limited zti_g as is & mask 558 ! Inside ML replace by linearly reducing sx_mlb towards surface & mask 555 559 ! 556 560 zfacti = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)), wp ) ! k index of uppermost point(s) of triad is jk+kp-1 557 561 zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1 558 562 ! ! otherwise zfact=0 559 zti_g_lim = zfacti * zti_g_lim &563 zti_g_lim = ( zfacti * zti_g_lim & 560 564 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 561 & * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) 562 ztj_g_lim = zfactj * ztj_g_lim &565 & * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 566 ztj_g_lim = ( zfactj * ztj_g_lim & 563 567 & + ( 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)567 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim * vmask(ji,jj,jk+kp)568 & * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 569 ! 570 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim 571 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim 568 572 ! 569 573 ! Get coefficients of isoneutral diffusion tensor … … 574 578 ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 575 579 ! 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 ! 580 zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces 581 ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 582 ! 583 IF( ln_triad_iso ) THEN 584 zti_raw = ( zti_lim*zti_lim ) / zti_raw 585 ztj_raw = ( ztj_lim*ztj_lim ) / ztj_raw 586 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 587 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 588 zti_lim = zfacti * zti_lim & 589 & + ( 1._wp - zfacti ) * zti_raw 590 ztj_lim = zfactj * ztj_lim & 591 & + ( 1._wp - zfactj ) * ztj_raw 592 ENDIF 593 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim 594 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim 595 ! 581 596 zbu = e1u(ji ,jj) * e2u(ji ,jj) * fse3u(ji ,jj,jk ) 582 597 zbv = e1v(ji ,jj) * e2v(ji ,jj) * fse3v(ji ,jj,jk ) … … 584 599 zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 585 600 ! 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 591 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_lim2 601 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked 602 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * ( zti_g_lim * zti_g_lim ) ! masked 603 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ( ztj_g_lim * ztj_g_lim ) 592 604 END DO 593 605 END DO … … 597 609 ! 598 610 wslp2(:,:,1) = 0._wp ! force the surface wslp to zero 599 611 600 612 CALL lbc_lnk( wslp2, 'W', 1. ) ! lateral boundary confition on wslp2 only ==>>> gm : necessary ? to be checked 601 613 ! … … 610 622 !! *** ROUTINE ldf_slp_mxl *** 611 623 !! 612 !! ** Purpose : Compute the slopes of iso-neutral surface just below 624 !! ** Purpose : Compute the slopes of iso-neutral surface just below 613 625 !! the mixed layer. 614 626 !! … … 619 631 !! 620 632 !! ** Action : uslpml, wslpiml : i- & j-slopes of neutral surfaces 621 !! vslpml, wslpjml just below the mixed layer 633 !! vslpml, wslpjml just below the mixed layer 622 634 !! omlmask : mixed layer mask 623 635 !!---------------------------------------------------------------------- … … 627 639 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_dzr ! z-gradient of density (T-point) 628 640 !! 629 INTEGER :: ji , jj , jk ! dummy loop indices630 INTEGER :: iku, ikv, ik, ikm1 ! local integers641 INTEGER :: ji , jj , jk ! dummy loop indices 642 INTEGER :: iku, ikv, ik, ikm1 ! local integers 631 643 REAL(wp) :: zeps, zm1_g, zm1_2g ! local scalars 632 644 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - … … 644 656 wslpjml(1,:) = 0._wp ; wslpjml(jpi,:) = 0._wp 645 657 ! 646 ! !== surface mixed layer mask !647 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise658 ! !== surface mixed layer mask ! 659 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise 648 660 # if defined key_vectopt_loop 649 661 DO jj = 1, 1 650 DO ji = 1, jpij ! vector opt. (forced unrolling)662 DO ji = 1, jpij ! vector opt. (forced unrolling) 651 663 # else 652 664 DO jj = 1, jpj … … 679 691 DO ji = 2, jpim1 680 692 # endif 681 ! !== Slope at u- & v-points just below the Mixed Layer ==!693 ! !== Slope at u- & v-points just below the Mixed Layer ==! 682 694 ! 683 ! 695 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 684 696 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 ) ! 697 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 686 698 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 687 699 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 688 ! 700 ! !- horizontal density gradient at u- & v-points 689 701 zau = p_gru(ji,jj,iku) / e1u(ji,jj) 690 702 zav = p_grv(ji,jj,ikv) / e2v(ji,jj) 691 ! 692 ! 703 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 704 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 693 705 zbu = MIN( zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) ) 694 706 zbv = MIN( zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) ) 695 ! 707 ! !- Slope at u- & v-points (uslpml, vslpml) 696 708 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 697 709 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 698 710 ! 699 ! !== i- & j-slopes at w-points just below the Mixed Layer ==!711 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! 700 712 ! 701 713 ik = MIN( nmln(ji,jj) + 1, jpk ) 702 714 ikm1 = MAX( 1, ik-1 ) 703 ! 715 ! !- vertical density gradient for w-slope (from N^2) 704 716 zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 705 ! 717 ! !- horizontal density i- & j-gradient at w-points 706 718 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) 719 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 708 720 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 709 721 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) … … 712 724 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & 713 725 & + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik) 714 ! 715 ! 726 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 727 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 716 728 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zai ) ) 717 729 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj ) ) 718 ! 730 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 719 731 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 720 732 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 721 733 END DO 722 734 END DO 723 !!gm this lbc_lnk should be useless....735 !!gm this lbc_lnk should be useless.... 724 736 CALL lbc_lnk( uslpml , 'U', -1. ) ; CALL lbc_lnk( vslpml , 'V', -1. ) ! lateral boundary cond. (sign change) 725 737 CALL lbc_lnk( wslpiml, 'W', -1. ) ; CALL lbc_lnk( wslpjml, 'W', -1. ) ! lateral boundary conditions … … 734 746 !! ** Purpose : Initialization for the isopycnal slopes computation 735 747 !! 736 !! ** Method : read the nammbf namelist and check the parameter 737 !! 748 !! ** Method : read the nammbf namelist and check the parameter 749 !! values called by tra_dmp at the first timestep (nit000) 738 750 !!---------------------------------------------------------------------- 739 751 INTEGER :: ji, jj, jk ! dummy loop indices 740 752 INTEGER :: ierr ! local integer 741 753 !!---------------------------------------------------------------------- 742 743 IF(lwp) THEN 754 755 IF(lwp) THEN 744 756 WRITE(numout,*) 745 757 WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing' 746 758 WRITE(numout,*) '~~~~~~~~~~~~' 747 759 ENDIF 748 760 749 761 IF( ln_traldf_grif ) THEN ! Griffies operator : triad of slopes 750 762 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 767 IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 756 768 ! 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 769 ELSE ! Madec operator : slopes at u-, v-, and w-points 761 770 ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) , & … … 770 779 wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp 771 780 772 !!gm I no longer understand this.....781 !!gm I no longer understand this..... 773 782 IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 774 783 IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' … … 803 812 LOGICAL, PUBLIC, PARAMETER :: lk_ldfslp = .FALSE. !: slopes flag 804 813 CONTAINS 805 SUBROUTINE ldf_slp( kt, prd, pn2 ) 806 INTEGER, INTENT(in) :: kt 814 SUBROUTINE ldf_slp( kt, prd, pn2 ) ! Dummy routine 815 INTEGER, INTENT(in) :: kt 807 816 REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2 808 817 WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1) … … 812 821 WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt 813 822 END SUBROUTINE ldf_slp_grif 814 SUBROUTINE ldf_slp_init ! Dummy routine823 SUBROUTINE ldf_slp_init ! Dummy routine 815 824 END SUBROUTINE ldf_slp_init 816 825 #endif
Note: See TracChangeset
for help on using the changeset viewer.
