Changeset 3294 for trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
- Timestamp:
- 2012-01-28T17:44:18+01:00 (12 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90 (modified) (35 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r2772 r3294 31 31 USE in_out_manager ! I/O manager 32 32 USE prtctl ! Print control 33 USE wrk_nemo ! work arrays 34 USE timing ! Timing 33 35 34 36 IMPLICIT NONE … … 46 48 ! !! Griffies operator 47 49 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 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi_g, triadj_g !: skew flux slopes relative to geopotentials 49 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi , triadj !: isoneutral slopes relative to model-coordinate 50 52 … … 56 58 57 59 REAL(wp) :: repsln = 1.e-25_wp ! tiny value used as minium of di(rho), dj(rho) and dk(rho) 58 59 ! Workspace arrays for ldf_slp_grif. These could be replaced by several 3D and 2D workspace60 ! arrays from the wrk_nemo module with a bit of code re-writing. The 4D workspace61 ! arrays can't be used here because of the zero-indexing of some of the ranks. ARPDBG.62 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: zdzrho , zdyrho, zdxrho ! Horizontal and vertical density gradients63 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: zti_mlb, ztj_mlb ! for Griffies operator only64 60 65 61 !! * Substitutions … … 75 71 CONTAINS 76 72 77 INTEGER FUNCTION ldf_slp_alloc()78 !!----------------------------------------------------------------------79 !! *** FUNCTION ldf_slp_alloc ***80 !!----------------------------------------------------------------------81 !82 ALLOCATE( zdxrho (jpi,jpj,jpk,0:1) , zti_mlb(jpi,jpj,0:1,0:1) , &83 & zdyrho (jpi,jpj,jpk,0:1) , ztj_mlb(jpi,jpj,0:1,0:1) , &84 & zdzrho (jpi,jpj,jpk,0:1) , STAT=ldf_slp_alloc )85 !86 IF( lk_mpp ) CALL mpp_sum ( ldf_slp_alloc )87 IF( ldf_slp_alloc /= 0 ) CALL ctl_warn('ldf_slp_alloc : failed to allocate arrays.')88 !89 END FUNCTION ldf_slp_alloc90 91 92 73 SUBROUTINE ldf_slp( kt, prd, pn2 ) 93 74 !!---------------------------------------------------------------------- 94 75 !! *** ROUTINE ldf_slp *** 95 !! 76 !! 96 77 !! ** Purpose : Compute the slopes of neutral surface (slope of isopycnal 97 78 !! 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 79 !! 80 !! ** Method : The slope in the i-direction is computed at U- and 81 !! W-points (uslp, wslpi) and the slope in the j-direction is 101 82 !! computed at V- and W-points (vslp, wslpj). 102 83 !! They are bounded by 1/100 over the whole ocean, and within the … … 112 93 !! bottom slope (ln_sco=T) at level jpk in inildf] 113 94 !! 114 !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and j-slopes 95 !! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and j-slopes 115 96 !! of now neutral surfaces at u-, w- and v- w-points, resp. 116 97 !!---------------------------------------------------------------------- 117 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released118 USE oce , ONLY: zgru => ua , zww => va ! (ua,va) used as workspace119 USE oce , ONLY: zgrv => ta , zwz => sa ! (ta,sa) used as workspace120 USE wrk_nemo, ONLY: zdzr => wrk_3d_1 ! 3D workspace121 !!122 98 INTEGER , INTENT(in) :: kt ! ocean time-step index 123 99 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: prd ! in situ density … … 127 103 INTEGER :: ii0, ii1, iku ! temporary integer 128 104 INTEGER :: ij0, ij1, ikv ! temporary integer 129 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16 105 REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16, zcofw ! local scalars 130 106 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - 131 107 REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - - 132 108 REAL(wp) :: zck, zfk, zbw ! - - 133 !!---------------------------------------------------------------------- 134 135 IF( wrk_in_use(3, 1) ) THEN 136 CALL ctl_stop('ldf_slp: requested workspace arrays are unavailable') ; RETURN 137 ENDIF 109 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zww 110 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdzr 111 REAL(wp), POINTER, DIMENSION(:,:,:) :: zgru, zgrv 112 !!---------------------------------------------------------------------- 113 ! 114 IF( nn_timing == 1 ) CALL timing_start('ldf_slp') 115 ! 116 CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 138 117 139 118 zeps = 1.e-20_wp !== Local constant initialization ==! … … 148 127 DO jj = 1, jpjm1 149 128 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) ) 129 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 130 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 152 131 END DO 153 132 END DO 154 133 END DO 155 134 IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level 156 # if defined key_vectopt_loop 135 # if defined key_vectopt_loop 157 136 DO jj = 1, 1 158 137 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) … … 161 140 DO ji = 1, jpim1 162 141 # endif 163 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 164 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 142 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 143 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 165 144 END DO 166 145 END DO … … 181 160 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml 182 161 183 162 184 163 ! I. slopes at u and v point | uslp = d/di( prd ) / d/dz( prd ) 185 164 ! =========================== | vslp = d/dj( prd ) / d/dz( prd ) 186 ! 165 ! 187 166 DO jk = 2, jpkm1 !* Slopes at u and v points 188 167 DO jj = 2, jpjm1 … … 225 204 DO jk = 2, jpkm1 226 205 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 227 DO ji = 2, jpim1 206 DO ji = 2, jpim1 228 207 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 229 208 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & … … 266 245 ! II. slopes at w point | wslpi = mij( d/di( prd ) / d/dz( prd ) 267 246 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 268 ! 247 ! 269 248 DO jk = 2, jpkm1 270 249 DO jj = 2, jpjm1 … … 308 287 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 309 288 DO ji = 2, jpim1 289 zcofw = tmask(ji,jj,jk) * z1_16 310 290 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)291 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 292 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 293 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 294 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 315 295 316 296 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 297 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 298 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 299 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 300 & + 4.* zww(ji ,jj ,jk) ) * zcofw 301 END DO 302 END DO 323 303 DO jj = 3, jpj-2 ! other rows 324 304 DO ji = fs_2, fs_jpim1 ! vector opt. 305 zcofw = tmask(ji,jj,jk) * z1_16 325 306 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)307 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 308 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 309 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 310 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 330 311 331 312 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)313 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 314 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 315 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 316 & + 4.* zww(ji ,jj ,jk) ) * zcofw 336 317 END DO 337 318 END DO … … 346 327 END DO 347 328 END DO 348 349 ! III. Specific grid points 350 ! =========================== 351 ! 329 330 ! III. Specific grid points 331 ! =========================== 332 ! 352 333 IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN ! ORCA_R4 configuration: horizontal diffusion in specific area 353 334 ! ! Gibraltar Strait … … 368 349 ENDIF 369 350 370 ! IV. Lateral boundary conditions 351 352 ! IV. Lateral boundary conditions 371 353 ! =============================== 372 354 CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) … … 379 361 ENDIF 380 362 ! 381 IF( wrk_not_released(3, 1) ) CALL ctl_stop('ldf_slp: failed to release workspace arrays') 363 CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv ) 364 ! 365 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp') 382 366 ! 383 367 END SUBROUTINE ldf_slp 384 368 385 369 386 370 SUBROUTINE ldf_slp_grif ( kt ) … … 390 374 !! ** Purpose : Compute the squared slopes of neutral surfaces (slope 391 375 !! 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 376 !! at W-points using the Griffies quarter-cells. 377 !! 378 !! ** Method : calculates alpha and beta at T-points 395 379 !! 396 380 !! ** Action : - triadi_g, triadj_g T-pts i- and j-slope triads relative to geopot. (used for eiv) … … 398 382 !! - wslp2 squared slope of neutral surfaces at w-points. 399 383 !!---------------------------------------------------------------------- 400 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 405 USE wrk_nemo, ONLY: z1_mlbw => wrk_2d_1 406 ! 407 INTEGER, INTENT( in ) :: kt ! ocean time-step index 408 ! 384 INTEGER, INTENT( in ) :: kt ! ocean time-step index 385 !! 409 386 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 387 INTEGER :: iku, ikv ! local integer 388 REAL(wp) :: zfacti, zfactj ! local scalars 389 REAL(wp) :: znot_thru_surface ! local scalars 390 REAL(wp) :: zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj 391 REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_g_raw, zti_g_lim 392 REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim 415 393 REAL(wp) :: zdzrho_raw 416 !!---------------------------------------------------------------------- 417 418 IF( wrk_in_use(3, 2,3,4,5) .OR. wrk_in_use(2, 1) )THEN 419 CALL ctl_stop('ldf_slp_grif: requested workspace arrays are unavailable') ; RETURN 420 ENDIF 421 394 REAL(wp) :: zbeta0 395 REAL(wp), POINTER, DIMENSION(:,:) :: z1_mlbw 396 REAL(wp), POINTER, DIMENSION(:,:,:) :: zalbet 397 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zdxrho , zdyrho, zdzrho ! Horizontal and vertical density gradients 398 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zti_mlb, ztj_mlb ! for Griffies operator only 399 !!---------------------------------------------------------------------- 400 ! 401 IF( nn_timing == 1 ) CALL timing_start('ldf_slp_grif') 402 ! 403 CALL wrk_alloc( jpi,jpj, z1_mlbw ) 404 CALL wrk_alloc( jpi,jpj,jpk, zalbet ) 405 CALL wrk_alloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho, klstart = 0 ) 406 CALL wrk_alloc( jpi,jpj, 2,2, zti_mlb, ztj_mlb, kkstart = 0, klstart = 0 ) 407 ! 422 408 !--------------------------------! 423 409 ! Some preliminary calculation ! 424 410 !--------------------------------! 425 411 ! 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 412 CALL eos_alpbet( tsb, zalbet, zbeta0 ) !== before local thermal/haline expension ratio at T-points ==! 413 ! 414 DO jl = 0, 1 !== unmasked before density i- j-, k-gradients ==! 415 ! 416 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln) 417 DO jk = 1, jpkm1 ! done each pair of triad 418 DO jj = 1, jpjm1 ! NB: not masked ==> a minimum value is set 419 DO ji = 1, fs_jpim1 ! vector opt. 420 zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! i-gradient of T & S at u-point 421 zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 422 zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! j-gradient of T & S at v-point 423 zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 424 zdxrho_raw = ( - zalbet(ji+ip,jj ,jk) * zdit + zbeta0*zdis ) / e1u(ji,jj) 425 zdyrho_raw = ( - zalbet(ji ,jj+jp,jk) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 426 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 427 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 428 END DO 429 END DO 430 END DO 431 ! 432 IF( ln_zps.and.l_grad_zps ) THEN ! partial steps: correction of i- & j-grad on bottom 439 433 # if defined key_vectopt_loop 440 DO jj = 1, 1441 DO ji = 1, jpij-jpi! vector opt. (forced unrolling)434 DO jj = 1, 1 435 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 442 436 # else 443 DO jj = 1, jpjm1444 DO ji = 1, jpim1437 DO jj = 1, jpjm1 438 DO ji = 1, jpim1 445 439 # 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 ) 440 iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points) 441 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 442 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 443 zdxrho_raw = ( - zalbet(ji+ip,jj ,iku) * zdit + zbeta0*zdis ) / e1u(ji,jj) 444 zdyrho_raw = ( - zalbet(ji ,jj+jp,ikv) * zdjt + zbeta0*zdjs ) / e2v(ji,jj) 445 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 446 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 471 447 END DO 472 448 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 ==! 449 ENDIF 450 ! 451 END DO 452 453 DO kp = 0, 1 !== unmasked before density i- j-, k-gradients ==! 454 DO jk = 1, jpkm1 ! done each pair of triad 455 DO jj = 1, jpj ! NB: not masked ==> a minimum value is set 456 DO ji = 1, jpi ! vector opt. 457 IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp 458 zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) ) 459 zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) ) 460 ELSE 461 zdkt = 0._wp ! 1st level gradient set to zero 462 zdks = 0._wp 463 ENDIF 464 zdzrho_raw = ( - zalbet(ji ,jj ,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp) 465 zdzrho(ji ,jj ,jk, kp) = - MIN( - repsln, zdzrho_raw ) ! force zdzrho >= repsln 466 END DO 467 END DO 468 END DO 469 END DO 470 ! 471 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==! 488 472 DO ji = 1, jpi 489 473 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth … … 492 476 END DO 493 477 ! 494 ! !== intialisations to zero ==!495 ! 496 wslp2 (:,:,:) = 0._wp 497 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp 478 ! !== intialisations to zero ==! 479 ! 480 wslp2 (:,:,:) = 0._wp ! wslp2 will be cumulated 3D field set to zero 481 triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 498 482 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 483 !!gm _iso set to zero missing 484 triadi (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero 485 triadj (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp 486 503 487 !-------------------------------------! 504 488 ! Triads just below the Mixed Layer ! 505 489 !-------------------------------------! 506 490 ! 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 MASKED491 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 492 DO kp = 0, 1 ! with only the slope-max limit and MASKED 509 493 DO jj = 1, jpjm1 510 494 DO ji = 1, fs_jpim1 511 495 ip = jl ; jp = jl 512 496 jk = MIN( nmln(ji+ip,jj) , mbkt(ji+ip,jj) ) + 1 ! ML level+1 (MIN in case ML depth is the ocean depth) 497 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 513 498 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)499 & - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj) ) * umask(ji,jj,jk) 515 500 jk = MIN( nmln(ji,jj+jp) , mbkt(ji,jj+jp) ) + 1 516 501 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)502 & - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 518 503 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 519 504 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) … … 527 512 !-------------------------------------! 528 513 ! 529 DO kp = 0, 1 ! k-index of triads514 DO kp = 0, 1 ! k-index of triads 530 515 DO jl = 0, 1 531 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes)516 ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes) 532 517 DO jk = 1, jpkm1 518 ! Must mask contribution to slope from dz/dx at constant s for triads jk=1,kp=0 that poke up though ocean surface 519 znot_thru_surface = REAL( 1-1/(jk+kp), wp ) !jk+kp=1,=0.; otherwise=1.0 533 520 DO jj = 1, jpjm1 534 DO ji = 1, fs_jpim1 ! vector opt.521 DO ji = 1, fs_jpim1 ! vector opt. 535 522 ! 536 523 ! 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)524 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 538 525 ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 539 526 ! masked by umask taken at the level of dz(rho) … … 543 530 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 544 531 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 532 533 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 534 zti_coord = znot_thru_surface * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj) 535 ztj_coord = znot_thru_surface * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) ! unmasked 536 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 537 ztj_g_raw = ztj_raw - ztj_coord 550 538 zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw ) 551 539 ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw ) 552 540 ! 553 ! Below ML use limited zti_g as is 554 ! Inside ML replace by linearly reducing sx_mlb towards surface 541 ! Below ML use limited zti_g as is & mask 542 ! Inside ML replace by linearly reducing sx_mlb towards surface & mask 555 543 ! 556 544 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 545 zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1 558 546 ! ! otherwise zfact=0 559 zti_g_lim = zfacti * zti_g_lim &547 zti_g_lim = ( zfacti * zti_g_lim & 560 548 & + ( 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 &549 & * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 550 ztj_g_lim = ( zfactj * ztj_g_lim & 563 551 & + ( 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)552 & * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 553 ! 554 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim 555 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim 568 556 ! 569 557 ! Get coefficients of isoneutral diffusion tensor … … 574 562 ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 575 563 ! 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 ! 564 zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces 565 ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 566 ! 567 IF( ln_triad_iso ) THEN 568 zti_raw = zti_lim**2 / zti_raw 569 ztj_raw = ztj_lim**2 / ztj_raw 570 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 571 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 572 zti_lim = zfacti * zti_lim & 573 & + ( 1._wp - zfacti ) * zti_raw 574 ztj_lim = zfactj * ztj_lim & 575 & + ( 1._wp - zfactj ) * ztj_raw 576 ENDIF 577 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim 578 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim 579 ! 581 580 zbu = e1u(ji ,jj) * e2u(ji ,jj) * fse3u(ji ,jj,jk ) 582 581 zbv = e1v(ji ,jj) * e2v(ji ,jj) * fse3v(ji ,jj,jk ) … … 584 583 zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp) 585 584 ! 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 585 !!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked 586 wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim**2 ! masked 587 wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim**2 592 588 END DO 593 589 END DO … … 597 593 ! 598 594 wslp2(:,:,1) = 0._wp ! force the surface wslp to zero 599 595 600 596 CALL lbc_lnk( wslp2, 'W', 1. ) ! lateral boundary confition on wslp2 only ==>>> gm : necessary ? to be checked 601 597 ! 602 IF( wrk_not_released(3, 2,3,4,5) .OR. & 603 wrk_not_released(2, 1) ) CALL ctl_stop('ldf_slp_grif: failed to release workspace arrays') 598 CALL wrk_dealloc( jpi,jpj, z1_mlbw ) 599 CALL wrk_dealloc( jpi,jpj,jpk, zalbet ) 600 CALL wrk_dealloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho, klstart = 0 ) 601 CALL wrk_dealloc( jpi,jpj, 2,2, zti_mlb, ztj_mlb, kkstart = 0, klstart = 0 ) 602 ! 603 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_grif') 604 604 ! 605 605 END SUBROUTINE ldf_slp_grif … … 610 610 !! *** ROUTINE ldf_slp_mxl *** 611 611 !! 612 !! ** Purpose : Compute the slopes of iso-neutral surface just below 612 !! ** Purpose : Compute the slopes of iso-neutral surface just below 613 613 !! the mixed layer. 614 614 !! … … 619 619 !! 620 620 !! ** Action : uslpml, wslpiml : i- & j-slopes of neutral surfaces 621 !! vslpml, wslpjml just below the mixed layer 621 !! vslpml, wslpjml just below the mixed layer 622 622 !! omlmask : mixed layer mask 623 623 !!---------------------------------------------------------------------- … … 627 627 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_dzr ! z-gradient of density (T-point) 628 628 !! 629 INTEGER :: ji , jj , jk ! dummy loop indices630 INTEGER :: iku, ikv, ik, ikm1 ! local integers629 INTEGER :: ji , jj , jk ! dummy loop indices 630 INTEGER :: iku, ikv, ik, ikm1 ! local integers 631 631 REAL(wp) :: zeps, zm1_g, zm1_2g ! local scalars 632 632 REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - - … … 634 634 REAL(wp) :: zck, zfk, zbw ! - - 635 635 !!---------------------------------------------------------------------- 636 636 ! 637 IF( nn_timing == 1 ) CALL timing_start('ldf_slp_mxl') 638 ! 637 639 zeps = 1.e-20_wp !== Local constant initialization ==! 638 640 zm1_g = -1.0_wp / grav … … 644 646 wslpjml(1,:) = 0._wp ; wslpjml(jpi,:) = 0._wp 645 647 ! 646 ! !== surface mixed layer mask !647 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise648 ! !== surface mixed layer mask ! 649 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise 648 650 # if defined key_vectopt_loop 649 651 DO jj = 1, 1 650 DO ji = 1, jpij ! vector opt. (forced unrolling)652 DO ji = 1, jpij ! vector opt. (forced unrolling) 651 653 # else 652 654 DO jj = 1, jpj … … 679 681 DO ji = 2, jpim1 680 682 # endif 681 ! !== Slope at u- & v-points just below the Mixed Layer ==!683 ! !== Slope at u- & v-points just below the Mixed Layer ==! 682 684 ! 683 ! 685 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 684 686 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 ) ! 687 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 686 688 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 687 689 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 688 ! 690 ! !- horizontal density gradient at u- & v-points 689 691 zau = p_gru(ji,jj,iku) / e1u(ji,jj) 690 692 zav = p_grv(ji,jj,ikv) / e2v(ji,jj) 691 ! 692 ! 693 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 694 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 693 695 zbu = MIN( zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) ) 694 696 zbv = MIN( zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) ) 695 ! 697 ! !- Slope at u- & v-points (uslpml, vslpml) 696 698 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 697 699 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 698 700 ! 699 ! !== i- & j-slopes at w-points just below the Mixed Layer ==!701 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! 700 702 ! 701 703 ik = MIN( nmln(ji,jj) + 1, jpk ) 702 704 ikm1 = MAX( 1, ik-1 ) 703 ! 705 ! !- vertical density gradient for w-slope (from N^2) 704 706 zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 705 ! 707 ! !- horizontal density i- & j-gradient at w-points 706 708 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) 709 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 708 710 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 709 711 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) … … 712 714 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & 713 715 & + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik) 714 ! 715 ! 716 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 717 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 716 718 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zai ) ) 717 719 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj ) ) 718 ! 720 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 719 721 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 720 722 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 721 723 END DO 722 724 END DO 723 !!gm this lbc_lnk should be useless....725 !!gm this lbc_lnk should be useless.... 724 726 CALL lbc_lnk( uslpml , 'U', -1. ) ; CALL lbc_lnk( vslpml , 'V', -1. ) ! lateral boundary cond. (sign change) 725 727 CALL lbc_lnk( wslpiml, 'W', -1. ) ; CALL lbc_lnk( wslpjml, 'W', -1. ) ! lateral boundary conditions 726 728 ! 729 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_mxl') 730 ! 727 731 END SUBROUTINE ldf_slp_mxl 728 732 … … 734 738 !! ** Purpose : Initialization for the isopycnal slopes computation 735 739 !! 736 !! ** Method : read the nammbf namelist and check the parameter 737 !! 740 !! ** Method : read the nammbf namelist and check the parameter 741 !! values called by tra_dmp at the first timestep (nit000) 738 742 !!---------------------------------------------------------------------- 739 743 INTEGER :: ji, jj, jk ! dummy loop indices 740 744 INTEGER :: ierr ! local integer 741 745 !!---------------------------------------------------------------------- 742 743 IF(lwp) THEN 746 ! 747 IF( nn_timing == 1 ) CALL timing_start('ldf_slp_init') 748 ! 749 IF(lwp) THEN 744 750 WRITE(numout,*) 745 751 WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing' 746 752 WRITE(numout,*) '~~~~~~~~~~~~' 747 753 ENDIF 748 754 749 755 IF( ln_traldf_grif ) THEN ! Griffies operator : triad of slopes 750 756 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 ) 751 757 ALLOCATE( triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , STAT=ierr ) 752 758 IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' ) 753 IF( ldf_slp_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate workspace arrays' )754 759 ! 755 760 IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' ) 756 !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 761 ! 760 762 ELSE ! Madec operator : slopes at u-, v-, and w-points … … 770 772 wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp 771 773 772 !!gm I no longer understand this.....774 !!gm I no longer understand this..... 773 775 IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 774 776 IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' … … 795 797 ENDIF 796 798 ! 799 IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_init') 800 ! 797 801 END SUBROUTINE ldf_slp_init 798 802 … … 803 807 LOGICAL, PUBLIC, PARAMETER :: lk_ldfslp = .FALSE. !: slopes flag 804 808 CONTAINS 805 SUBROUTINE ldf_slp( kt, prd, pn2 ) 806 INTEGER, INTENT(in) :: kt 809 SUBROUTINE ldf_slp( kt, prd, pn2 ) ! Dummy routine 810 INTEGER, INTENT(in) :: kt 807 811 REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2 808 812 WRITE(*,*) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1) … … 812 816 WRITE(*,*) 'ldf_slp_grif: You should not have seen this print! error?', kt 813 817 END SUBROUTINE ldf_slp_grif 814 SUBROUTINE ldf_slp_init ! Dummy routine818 SUBROUTINE ldf_slp_init ! Dummy routine 815 819 END SUBROUTINE ldf_slp_init 816 820 #endif
Note: See TracChangeset
for help on using the changeset viewer.
