- Timestamp:
- 2011-09-15T14:38:15+02:00 (13 years ago)
- File:
-
- 1 edited
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)
Note: See TracChangeset
for help on using the changeset viewer.
