Changeset 10874 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90
- Timestamp:
- 2019-04-15T15:57:37+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traldf_triad.F90
r10806 r10874 48 48 CONTAINS 49 49 50 SUBROUTINE tra_ldf_triad( kt, kit000, ktlev, kt2lev,cdtype, pahu, pahv, pgu , pgv , &50 SUBROUTINE tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu , pgv , & 51 51 & pgui, pgvi, & 52 & pt , pt_lev0, pt_rhs, kjpt, kpass )52 & ptb , ptbb, pta , kjpt, kpass ) 53 53 !!---------------------------------------------------------------------- 54 54 !! *** ROUTINE tra_ldf_triad *** … … 66 66 !! see documentation for the desciption 67 67 !! 68 !! ** Action : pt _rhsupdated with the before rotated diffusion68 !! ** Action : pta updated with the before rotated diffusion 69 69 !! ah_wslp2 .... 70 70 !! akz stabilizing vertical diffusivity coefficient (used in trazdf_imp) … … 72 72 INTEGER , INTENT(in ) :: kt ! ocean time-step index 73 73 INTEGER , INTENT(in ) :: kit000 ! first time step index 74 INTEGER , INTENT(in ) :: ktlev ! time level index for e3t75 INTEGER , INTENT(in ) :: kt2lev ! time level index for 2-time-level thicknesses76 74 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 77 75 INTEGER , INTENT(in ) :: kjpt ! number of tracers … … 80 78 REAL(wp), DIMENSION(jpi,jpj ,kjpt), INTENT(in ) :: pgu , pgv ! tracer gradient at pstep levels 81 79 REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT(in ) :: pgui, pgvi ! tracer gradient at top levels 82 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! tracer (kpass=1) or laplacian of tracer (kpass=2)83 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt _lev0! tracer (only used in kpass=2)84 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt _rhs! tracer trend80 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! tracer (kpass=1) or laplacian of tracer (kpass=2) 81 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptbb ! tracer (only used in kpass=2) 82 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 85 83 ! 86 84 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 144 142 DO jj = 1, jpjm1 145 143 DO ji = 1, fs_jpim1 146 ze3wr = 1._wp / e3w (ji+ip,jj,jk+kp,kt2lev)147 zbu = e1e2u(ji,jj) * e3u (ji,jj,jk,ktlev)144 ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 145 zbu = e1e2u(ji,jj) * e3u_n(ji,jj,jk) 148 146 zah = 0.25_wp * pahu(ji,jj,jk) 149 147 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 150 148 ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 151 zslope2 = zslope_skew + ( gdept (ji+1,jj,jk,kt2lev) - gdept(ji,jj,jk,kt2lev) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp)149 zslope2 = zslope_skew + ( gdept_n(ji+1,jj,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 152 150 zslope2 = zslope2 *zslope2 153 151 ah_wslp2(ji+ip,jj,jk+kp) = ah_wslp2(ji+ip,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji+ip,jj) * zslope2 … … 168 166 DO jj = 1, jpjm1 169 167 DO ji = 1, fs_jpim1 170 ze3wr = 1.0_wp / e3w (ji,jj+jp,jk+kp,kt2lev)171 zbv = e1e2v(ji,jj) * e3v (ji,jj,jk,ktlev)168 ze3wr = 1.0_wp / e3w_n(ji,jj+jp,jk+kp) 169 zbv = e1e2v(ji,jj) * e3v_n(ji,jj,jk) 172 170 zah = 0.25_wp * pahv(ji,jj,jk) 173 171 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 174 172 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 175 173 ! (do this by *adding* gradient of depth) 176 zslope2 = zslope_skew + ( gdept (ji,jj+1,jk,kt2lev) - gdept(ji,jj,jk,kt2lev) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp)174 zslope2 = zslope_skew + ( gdept_n(ji,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) 177 175 zslope2 = zslope2 * zslope2 178 176 ah_wslp2(ji,jj+jp,jk+kp) = ah_wslp2(ji,jj+jp,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj+jp) * zslope2 … … 195 193 DO ji = 1, fs_jpim1 196 194 akz(ji,jj,jk) = 16._wp * ah_wslp2(ji,jj,jk) & 197 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w (ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev) ) )195 & * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ( e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) ) ) 198 196 END DO 199 197 END DO … … 203 201 DO jj = 1, jpjm1 204 202 DO ji = 1, fs_jpim1 205 ze3w_2 = e3w (ji,jj,jk,kt2lev) * e3w(ji,jj,jk,kt2lev)203 ze3w_2 = e3w_n(ji,jj,jk) * e3w_n(ji,jj,jk) 206 204 zcoef0 = z2dt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) 207 205 akz(ji,jj,jk) = MAX( zcoef0 - 0.5_wp , 0._wp ) * ze3w_2 * z1_2dt … … 231 229 DO jj = 1, jpjm1 232 230 DO ji = 1, fs_jpim1 ! vector opt. 233 zdit(ji,jj,jk) = ( pt (ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk)234 zdjt(ji,jj,jk) = ( pt (ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk)231 zdit(ji,jj,jk) = ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) ) * umask(ji,jj,jk) 232 zdjt(ji,jj,jk) = ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 235 233 END DO 236 234 END DO … … 259 257 DO jk = 1, jpkm1 260 258 ! !== Vertical tracer gradient at level jk and jk+1 261 zdkt3d(:,:,1) = ( pt (:,:,jk,jn) - pt(:,:,jk+1,jn) ) * tmask(:,:,jk+1)259 zdkt3d(:,:,1) = ( ptb(:,:,jk,jn) - ptb(:,:,jk+1,jn) ) * tmask(:,:,jk+1) 262 260 ! 263 261 ! ! surface boundary condition: zdkt3d(jk=0)=zdkt3d(jk=1) 264 262 IF( jk == 1 ) THEN ; zdkt3d(:,:,0) = zdkt3d(:,:,1) 265 ELSE ; zdkt3d(:,:,0) = ( pt (:,:,jk-1,jn) - pt(:,:,jk,jn) ) * tmask(:,:,jk)263 ELSE ; zdkt3d(:,:,0) = ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) * tmask(:,:,jk) 266 264 ENDIF 267 265 ! … … 275 273 ze1ur = r1_e1u(ji,jj) 276 274 zdxt = zdit(ji,jj,jk) * ze1ur 277 ze3wr = 1._wp / e3w (ji+ip,jj,jk+kp,kt2lev)275 ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 278 276 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 279 277 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 280 278 zslope_iso = triadi (ji+ip,jj,jk,1-ip,kp) 281 279 ! 282 zbu = 0.25_wp * e1e2u(ji,jj) * e3u (ji,jj,jk,ktlev)280 zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk) 283 281 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked.... 284 282 zah = pahu(ji,jj,jk) … … 298 296 ze2vr = r1_e2v(ji,jj) 299 297 zdyt = zdjt(ji,jj,jk) * ze2vr 300 ze3wr = 1._wp / e3w (ji,jj+jp,jk+kp,kt2lev)298 ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp) 301 299 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 302 300 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 303 301 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 304 zbv = 0.25_wp * e1e2v(ji,jj) * e3v (ji,jj,jk,ktlev)302 zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk) 305 303 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahv is masked... 306 304 zah = pahv(ji,jj,jk) … … 322 320 ze1ur = r1_e1u(ji,jj) 323 321 zdxt = zdit(ji,jj,jk) * ze1ur 324 ze3wr = 1._wp / e3w (ji+ip,jj,jk+kp,kt2lev)322 ze3wr = 1._wp / e3w_n(ji+ip,jj,jk+kp) 325 323 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 326 324 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 327 325 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 328 326 ! 329 zbu = 0.25_wp * e1e2u(ji,jj) * e3u (ji,jj,jk,ktlev)327 zbu = 0.25_wp * e1e2u(ji,jj) * e3u_n(ji,jj,jk) 330 328 ! ln_botmix_triad is .F. mask zah for bottom half cells 331 329 zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp) ! pahu(ji+ip,jj,jk) ===>> ???? … … 345 343 ze2vr = r1_e2v(ji,jj) 346 344 zdyt = zdjt(ji,jj,jk) * ze2vr 347 ze3wr = 1._wp / e3w (ji,jj+jp,jk+kp,kt2lev)345 ze3wr = 1._wp / e3w_n(ji,jj+jp,jk+kp) 348 346 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 349 347 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 350 348 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 351 zbv = 0.25_wp * e1e2v(ji,jj) * e3v (ji,jj,jk,ktlev)349 zbv = 0.25_wp * e1e2v(ji,jj) * e3v_n(ji,jj,jk) 352 350 ! ln_botmix_triad is .F. mask zah for bottom half cells 353 351 zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! pahv(ji,jj+jp,jk) ???? … … 364 362 DO jj = 2 , jpjm1 365 363 DO ji = fs_2, fs_jpim1 ! vector opt. 366 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) &364 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( zftu(ji-1,jj,jk) - zftu(ji,jj,jk) & 367 365 & + zftv(ji,jj-1,jk) - zftv(ji,jj,jk) ) & 368 & / ( e1e2t(ji,jj) * e3t (ji,jj,jk,ktlev) )366 & / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 369 367 END DO 370 368 END DO … … 377 375 DO jj = 1, jpjm1 378 376 DO ji = fs_2, fs_jpim1 379 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w (ji,jj,jk,kt2lev) * tmask(ji,jj,jk) &377 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) & 380 378 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 381 & * ( pt (ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )379 & * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 382 380 END DO 383 381 END DO … … 389 387 DO jj = 1, jpjm1 390 388 DO ji = fs_2, fs_jpim1 391 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w (ji,jj,jk,kt2lev) * tmask(ji,jj,jk) &392 & * ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) )389 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) & 390 & * ah_wslp2(ji,jj,jk) * ( ptb(ji,jj,jk-1,jn) - ptb(ji,jj,jk,jn) ) 393 391 END DO 394 392 END DO 395 393 END DO 396 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt and pt_lev0gradients, resp.394 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on ptb and ptbb gradients, resp. 397 395 DO jk = 2, jpkm1 398 396 DO jj = 1, jpjm1 399 397 DO ji = fs_2, fs_jpim1 400 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w (ji,jj,jk,kt2lev) * tmask(ji,jj,jk) &401 & * ( ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) &402 & + akz (ji,jj,jk) * ( pt _lev0(ji,jj,jk-1,jn) - pt_lev0(ji,jj,jk,jn) ) )398 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w_n(ji,jj,jk) * tmask(ji,jj,jk) & 399 & * ( ah_wslp2(ji,jj,jk) * ( ptb (ji,jj,jk-1,jn) - ptb (ji,jj,jk,jn) ) & 400 & + akz (ji,jj,jk) * ( ptbb(ji,jj,jk-1,jn) - ptbb(ji,jj,jk,jn) ) ) 403 401 END DO 404 402 END DO … … 407 405 ENDIF 408 406 ! 409 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pt _rhs==!407 DO jk = 1, jpkm1 !== Divergence of vertical fluxes added to pta ==! 410 408 DO jj = 2, jpjm1 411 409 DO ji = fs_2, fs_jpim1 ! vector opt. 412 pt _rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) &413 & / ( e1e2t(ji,jj) * e3t (ji,jj,jk,ktlev) )410 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + zsign * ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) & 411 & / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) 414 412 END DO 415 413 END DO
Note: See TracChangeset
for help on using the changeset viewer.
