- Timestamp:
- 2021-04-07T19:16:18+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/traldf_triad.F90
r14632 r14680 107 107 REAL(wp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) :: pt_rhs ! tracer trend 108 108 ! 109 INTEGER :: ji, jj, jk, jn ! dummy loop indices 110 INTEGER :: ip,jp,kp ! dummy loop indices 111 INTEGER :: ierr, iij ! local integer 112 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 113 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 109 INTEGER :: ji, jj, jk, jn, kp, iij ! dummy loop indices 114 110 REAL(wp) :: zcoef0, ze3w_2, zsign ! - - 115 111 ! 116 REAL(wp) :: zslope_skew, zslope_iso, zslope2, zbu, zbv 117 REAL(wp) :: ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt 118 REAL(wp) :: zah, zah_slp, zaei_slp 119 REAL(wp), DIMENSION(A2D(nn_hls),0:1) :: zdkt3d ! vertical tracer gradient at 2 levels 120 REAL(wp), DIMENSION(A2D(nn_hls) ) :: z2d ! 2D workspace 121 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdit, zdjt, zpsi_uw, zpsi_vw ! 3D - 122 ! NOTE: [halo1-halo2] 123 REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: zftu, zftv 124 REAL(wp), DIMENSION(A2D(nn_hls),jpk,2,2) :: ztfw 112 REAL(wp) :: zslope2, zbu, zbv, zbu1, zbv1, zslope21, zah, zah1, zah_ip1, zah_jp1, zbu_ip1, zbv_jp1 113 REAL(wp) :: ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt, zdyt_jp1, ze3wr_jp1, zdzt_jp1, zah_slp1, zah_slp_jp1, zaei_slp_jp1 114 REAL(wp) :: zah_slp, zaei_slp, zdxt_ip1, ze3wr_ip1, zdzt_ip1, zah_slp_ip1, zaei_slp_ip1, zaei_slp1 115 REAL(wp), DIMENSION(A2D(nn_hls),0:1) :: zdkt3d ! vertical tracer gradient at 2 levels 116 REAL(wp), DIMENSION(A2D(nn_hls) ) :: z2d ! 2D workspace 117 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk) :: zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw ! 3D - 125 118 !!---------------------------------------------------------------------- 126 119 ! … … 163 156 END_3D 164 157 ! 165 ! NOTE: [halo1-halo2] 166 zftu(:,:,:,:) = 0._wp 167 zftv(:,:,:,:) = 0._wp 168 ! 169 DO ip = 0, 1 ! i-k triads 170 DO kp = 0, 1 171 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 172 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 173 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 174 zbu = e1e2u(ji-ip,jj) * e3u(ji-ip,jj,jk,Kmm) 175 zah = 0.25_wp * pahu(ji-ip,jj,jk) 176 zslope_skew = triadi_g(ji,jj,jk,1-ip,kp) 177 ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 178 zslope2 = zslope_skew + ( gdept(ji-ip+1,jj,jk,Kmm) - gdept(ji-ip,jj,jk,Kmm) ) * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp) 179 zslope2 = zslope2 *zslope2 180 ! NOTE: [halo1-halo2] 181 zftu(ji,jj,jk+kp,ip+1) = zftu(ji,jj,jk+kp,ip+1) + & 182 & zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2 183 zftv(ji,jj,jk+kp,ip+1) = zftv(ji,jj,jk+kp,ip+1) + & 184 & zah * r1_e1u(ji-ip,jj) * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp) 185 ! 186 END_3D 187 END DO 158 DO kp = 0, 1 ! i-k triads 159 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 160 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 161 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 162 zbu = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 163 zbu1 = e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) 164 zah = 0.25_wp * pahu(ji,jj,jk) 165 zah1 = 0.25_wp * pahu(ji-1,jj,jk) 166 ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 167 zslope2 = triadi_g(ji,jj,jk,1,kp) + ( gdept(ji+1,jj,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) 168 zslope2 = zslope2 *zslope2 169 zslope21 = triadi_g(ji,jj,jk,0,kp) + ( gdept(ji,jj,jk,Kmm) - gdept(ji-1,jj,jk,Kmm) ) * r1_e1u(ji-1,jj) * umask(ji-1,jj,jk+kp) 170 zslope21 = zslope21 *zslope21 171 ! round brackets added to fix the order of floating point operations 172 ! needed to ensure halo 1 - halo 2 compatibility 173 ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + ( zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2 & 174 & + zah1 * zbu1 * ze3wr * r1_e1e2t(ji,jj) * zslope21 & 175 & ) ! bracket for halo 1 - halo 2 compatibility 176 akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + ( zah * r1_e1u(ji,jj) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) & 177 + zah1 * r1_e1u(ji-1,jj) * r1_e1u(ji-1,jj) * umask(ji-1,jj,jk+kp) & 178 & ) ! bracket for halo 1 - halo 2 compatibility 179 END_3D 188 180 END DO 189 181 ! 190 ! NOTE: [halo1-halo2] Use DO_3D instead of DO_3D_OVR 191 ! ip loop contributions added here to ensure consistent floating point arithmetic for different nn_hls 192 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 193 ah_wslp2(ji,jj,jk) = ah_wslp2(ji,jj,jk) + (zftu(ji,jj,jk,1) + zftu(ji,jj,jk,2)) 194 akz (ji,jj,jk) = akz (ji,jj,jk) + (zftv(ji,jj,jk,1) + zftv(ji,jj,jk,2)) 195 END_3D 196 ! 197 ! NOTE: [halo1-halo2] 198 zftu(:,:,:,:) = 0._wp 199 zftv(:,:,:,:) = 0._wp 200 ! 201 DO jp = 0, 1 ! j-k triads 202 DO kp = 0, 1 203 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 204 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 205 ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm) 206 zbv = e1e2v(ji,jj-jp) * e3v(ji,jj-jp,jk,Kmm) 207 zah = 0.25_wp * pahv(ji,jj-jp,jk) 208 zslope_skew = triadj_g(ji,jj,jk,1-jp,kp) 209 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 210 ! (do this by *adding* gradient of depth) 211 zslope2 = zslope_skew + ( gdept(ji,jj-jp+1,jk,Kmm) - gdept(ji,jj-jp,jk,Kmm) ) * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp) 212 zslope2 = zslope2 * zslope2 213 ! NOTE: [halo1-halo2] 214 zftu(ji,jj,jk+kp,jp+1) = zftu(ji,jj,jk+kp,jp+1) + & 215 & zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2 216 zftv(ji,jj,jk+kp,jp+1) = zftv(ji,jj,jk+kp,jp+1) + & 217 & zah * r1_e2v(ji,jj-jp) * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp) 218 ! 219 END_3D 220 END DO 182 DO kp = 0, 1 ! j-k triads 183 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 184 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 185 ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm) 186 zbv = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 187 zbv1 = e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) 188 zah = 0.25_wp * pahv(ji,jj,jk) 189 zah1 = 0.25_wp * pahv(ji,jj-1,jk) 190 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 191 ! (do this by *adding* gradient of depth) 192 zslope2 = triadj_g(ji,jj,jk,1,kp) + ( gdept(ji,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) 193 zslope2 = zslope2 * zslope2 194 zslope21 = triadj_g(ji,jj,jk,0,kp) + ( gdept(ji,jj,jk,Kmm) - gdept(ji,jj-1,jk,Kmm) ) * r1_e2v(ji,jj-1) * vmask(ji,jj-1,jk+kp) 195 zslope21 = zslope21 * zslope21 196 ! round brackets added to fix the order of floating point operations 197 ! needed to ensure halo 1 - halo 2 compatibility 198 ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + ( zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2 & 199 & + zah1 * zbv1 * ze3wr * r1_e1e2t(ji,jj) * zslope21 & 200 & ) ! bracket for halo 1 - halo 2 compatibility 201 akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + ( zah * r1_e2v(ji,jj) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) & 202 & + zah1 * r1_e2v(ji,jj-1) * r1_e2v(ji,jj-1) * vmask(ji,jj-1,jk+kp) & 203 & ) ! bracket for halo 1 - halo 2 compatibility 204 ! 205 END_3D 221 206 END DO 222 !223 ! NOTE: [halo1-halo2] Use DO_3D instead of DO_3D_OVR224 ! jp loop contributions added here to ensure consistent floating point arithmetic for different nn_hls225 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )226 ah_wslp2(ji,jj,jk) = ah_wslp2(ji,jj,jk) + (zftu(ji,jj,jk,1) + zftu(ji,jj,jk,2))227 akz (ji,jj,jk) = akz (ji,jj,jk) + (zftv(ji,jj,jk,1) + zftv(ji,jj,jk,2))228 END_3D229 207 ! 230 208 IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient … … 259 237 zpsi_vw(:,:,:) = 0._wp 260 238 261 DO jp = 0, 1 262 DO kp = 0, 1 263 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 264 zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 265 & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+jp,jj,jk,1-jp,kp) 266 zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) & 267 & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+jp,jk,1-jp,kp) 268 END_3D 269 END DO 239 DO kp = 0, 1 240 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 241 ! round brackets added to fix the order of floating point operations 242 ! needed to ensure halo 1 - halo 2 compatibility 243 zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 244 & + ( 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji,jj,jk,1,kp) & 245 & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+1,jj,jk,0,kp) & 246 & ) ! bracket for halo 1 - halo 2 compatibility 247 zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) & 248 & + ( 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj,jk,1,kp) & 249 & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+1,jk,0,kp) & 250 & ) ! bracket for halo 1 - halo 2 compatibility 251 END_3D 270 252 END DO 271 253 CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) … … 279 261 ! Zero fluxes for each tracer 280 262 !!gm this should probably be done outside the jn loop 281 ! NOTE: [halo1-halo2] 282 ztfw(:,:,:,:,:) = 0._wp 283 zftu(:,:,:,:) = 0._wp 284 zftv(:,:,:,:) = 0._wp 263 ztfw(:,:,:) = 0._wp 264 zftu(:,:,:) = 0._wp 265 zftv(:,:,:) = 0._wp 285 266 ! 286 267 ! [comm_cleanup] ! DO_3D( 1, 0, 1, 0, 1, jpkm1 ) !== before lateral T & S gradients at T-level jk ==! … … 327 308 ! 328 309 IF( ln_botmix_triad ) THEN 329 DO ip = 0, 1 !== Horizontal & vertical fluxes 330 DO kp = 0, 1 331 ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 332 DO_2D( iij, iij-1, iij, iij-1 ) 333 ze1ur = r1_e1u(ji,jj) 334 zdxt = zdit(ji,jj,jk) * ze1ur 335 ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm) 336 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 337 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 338 zslope_iso = triadi (ji+ip,jj,jk,1-ip,kp) 339 ! 340 zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 341 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked.... 342 zah = pahu(ji,jj,jk) 343 zah_slp = zah * zslope_iso 344 IF( ln_ldfeiv ) zaei_slp = aeiu(ji,jj,jk) * zslope_skew 345 ! NOTE: [halo1-halo2] 346 zftu(ji ,jj,jk,ip+1 ) = zftu(ji ,jj,jk,ip+1 ) - & 347 & ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 348 ztfw(ji,jj,jk+kp,ip+1,1) = ztfw(ji,jj,jk+kp,ip+1,1) - & 349 & (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 350 END_2D 351 END DO 310 DO kp = 0, 1 !== Horizontal & vertical fluxes 311 ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 312 DO_2D( iij, iij-1, iij, iij-1 ) 313 ze1ur = r1_e1u(ji,jj) 314 zdxt = zdit(ji,jj,jk) * ze1ur 315 zdxt_ip1 = zdit(ji+1,jj,jk) * r1_e1u(ji+1,jj) 316 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 317 ze3wr_ip1 = 1._wp / e3w(ji+1,jj,jk+kp,Kmm) 318 zdzt = zdkt3d(ji,jj,kp) * ze3wr 319 zdzt_ip1 = zdkt3d(ji+1,jj,kp) * ze3wr_ip1 320 ! 321 zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 322 zbu_ip1 = 0.25_wp * e1e2u(ji+1,jj) * e3u(ji+1,jj,jk,Kmm) 323 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked.... 324 zah = pahu(ji,jj,jk) 325 zah_ip1 = pahu(ji+1,jj,jk) 326 zah_slp = zah * triadi(ji,jj,jk,1,kp) 327 zah_slp_ip1 = zah_ip1 * triadi(ji+1,jj,jk,1,kp) 328 zah_slp1 = zah * triadi(ji+1,jj,jk,0,kp) 329 IF( ln_ldfeiv ) THEN 330 zaei_slp = aeiu(ji,jj,jk) * triadi_g(ji,jj,jk,1,kp) 331 zaei_slp_ip1 = aeiu(ji+1,jj,jk) * triadi_g(ji+1,jj,jk,1,kp) 332 zaei_slp1 = aeiu(ji,jj,jk) * triadi_g(ji+1,jj,jk,0,kp) 333 ENDIF 334 ! round brackets added to fix the order of floating point operations 335 ! needed to ensure halo 1 - halo 2 compatibility 336 zftu(ji ,jj,jk ) = zftu(ji ,jj,jk ) & 337 & - ( ( zah * zdxt + ( zah_slp - zaei_slp ) * zdzt ) * zbu * ze1ur & 338 & + ( zah * zdxt + zah_slp1 * zdzt_ip1 - zaei_slp1 * zdzt_ip1 ) * zbu * ze1ur & 339 & ) ! bracket for halo 1 - halo 2 compatibility 340 ztfw(ji+1,jj,jk+kp) = ztfw(ji+1,jj,jk+kp) & 341 & - ( (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1 & 342 & + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1 & 343 & ) ! bracket for halo 1 - halo 2 compatibility 344 END_2D 352 345 END DO 353 346 ! 354 DO jp = 0, 1 355 DO kp = 0, 1 356 ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 357 DO_2D( iij, iij-1, iij, iij-1 ) 358 ze2vr = r1_e2v(ji,jj) 359 zdyt = zdjt(ji,jj,jk) * ze2vr 360 ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm) 361 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 362 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 363 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 364 zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 365 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahv is masked... 366 zah = pahv(ji,jj,jk) 367 zah_slp = zah * zslope_iso 368 IF( ln_ldfeiv ) zaei_slp = aeiv(ji,jj,jk) * zslope_skew 369 ! NOTE: [halo1-halo2] 370 zftv(ji,jj ,jk,jp+1 ) = zftv(ji,jj ,jk,jp+1 ) - & 371 & ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 372 ztfw(ji,jj,jk+kp,jp+1,2) = ztfw(ji,jj,jk+kp,jp+1,2) - & 373 & (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 374 END_2D 375 END DO 347 DO kp = 0, 1 348 ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 349 DO_2D( iij, iij-1, iij, iij-1 ) 350 ze2vr = r1_e2v(ji,jj) 351 zdyt = zdjt(ji,jj,jk) * ze2vr 352 zdyt_jp1 = zdjt(ji,jj+1,jk) * r1_e2v(ji,jj+1) 353 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 354 ze3wr_jp1 = 1._wp / e3w(ji,jj+1,jk+kp,Kmm) 355 zdzt = zdkt3d(ji,jj,kp) * ze3wr 356 zdzt_jp1 = zdkt3d(ji,jj+1,kp) * ze3wr_jp1 357 zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 358 zbv_jp1 = 0.25_wp * e1e2v(ji,jj+1) * e3v(ji,jj+1,jk,Kmm) 359 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked.... 360 zah = pahv(ji,jj,jk) ! pahv(ji,jj+jp,jk) ???? 361 zah_jp1 = pahv(ji,jj+1,jk) 362 zah_slp = zah * triadj(ji,jj,jk,1,kp) 363 zah_slp1 = zah * triadj(ji,jj+1,jk,0,kp) 364 zah_slp_jp1 = zah_jp1 * triadj(ji,jj+1,jk,1,kp) 365 IF( ln_ldfeiv ) THEN 366 zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp) 367 zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp) 368 zaei_slp1 = aeiv(ji,jj,jk) * triadj_g(ji,jj+1,jk,0,kp) 369 ENDIF 370 ! round brackets added to fix the order of floating point operations 371 ! needed to ensure halo 1 - halo 2 compatibility 372 zftv(ji,jj ,jk ) = zftv(ji,jj ,jk ) & 373 & - ( ( zah * zdyt + ( zah_slp - zaei_slp ) * zdzt ) * zbv * ze2vr & 374 & + ( zah * zdyt + zah_slp1 * zdzt_jp1 - zaei_slp1 * zdzt_jp1 ) * zbv * ze2vr & 375 & ) ! bracket for halo 1 - halo 2 compatibility 376 ztfw(ji,jj+1,jk+kp) = ztfw(ji,jj+1,jk+kp) & 377 & - ( ( zah_slp_jp1 + zaei_slp_jp1) * zdyt_jp1 * zbv_jp1 * ze3wr_jp1 & 378 & + ( zah_slp1 + zaei_slp1) * zdyt * zbv * ze3wr_jp1 & 379 & ) ! bracket for halo 1 - halo 2 compatibility 380 END_2D 376 381 END DO 377 382 ! 378 383 ELSE 379 384 ! 380 DO ip = 0, 1 !== Horizontal & vertical fluxes 381 DO kp = 0, 1 382 ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 383 DO_2D( iij, iij-1, iij, iij-1 ) 384 ze1ur = r1_e1u(ji,jj) 385 zdxt = zdit(ji,jj,jk) * ze1ur 386 ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm) 387 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 388 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 389 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 390 ! 391 zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 392 ! ln_botmix_triad is .F. mask zah for bottom half cells 393 zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp) ! pahu(ji+ip,jj,jk) ===>> ???? 394 zah_slp = zah * zslope_iso 395 IF( ln_ldfeiv ) zaei_slp = aeiu(ji,jj,jk) * zslope_skew ! aeit(ji+ip,jj,jk)*zslope_skew 396 ! NOTE: [halo1-halo2] 397 zftu(ji ,jj,jk,ip+1 ) = zftu(ji ,jj,jk,ip+1 ) - & 398 & ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 399 ztfw(ji,jj,jk+kp,ip+1,1) = ztfw(ji,jj,jk+kp,ip+1,1) - & 400 & (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 401 END_2D 402 END DO 385 DO kp = 0, 1 !== Horizontal & vertical fluxes 386 ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 387 DO_2D( iij, iij-1, iij, iij-1 ) 388 ze1ur = r1_e1u(ji,jj) 389 zdxt = zdit(ji,jj,jk) * ze1ur 390 zdxt_ip1 = zdit(ji+1,jj,jk) * r1_e1u(ji+1,jj) 391 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 392 ze3wr_ip1 = 1._wp / e3w(ji+1,jj,jk+kp,Kmm) 393 zdzt = zdkt3d(ji,jj,kp) * ze3wr 394 zdzt_ip1 = zdkt3d(ji+1,jj,kp) * ze3wr_ip1 395 ! 396 zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 397 zbu_ip1 = 0.25_wp * e1e2u(ji+1,jj) * e3u(ji+1,jj,jk,Kmm) 398 ! ln_botmix_triad is .F. mask zah for bottom half cells 399 zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp) ! pahu(ji+ip,jj,jk) ===>> ???? 400 zah_ip1 = pahu(ji+1,jj,jk) * umask(ji+1,jj,jk+kp) 401 zah_slp = zah * triadi(ji,jj,jk,1,kp) 402 zah_slp_ip1 = zah_ip1 * triadi(ji+1,jj,jk,1,kp) 403 zah_slp1 = zah * triadi(ji+1,jj,jk,0,kp) 404 IF( ln_ldfeiv ) THEN 405 zaei_slp = aeiu(ji,jj,jk) * triadi_g(ji,jj,jk,1,kp) 406 zaei_slp_ip1 = aeiu(ji+1,jj,jk) * triadi_g(ji+1,jj,jk,1,kp) 407 zaei_slp1 = aeiu(ji,jj,jk) * triadi_g(ji+1,jj,jk,0,kp) 408 ENDIF 409 ! zftu(ji ,jj,jk ) = zftu(ji ,jj,jk ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur - ( zah * zdxt + (zah_slp1 - zaei_slp1) * zdzt_ip1 ) * zbu * ze1ur 410 ! ztfw(ji+1,jj,jk+kp) = ztfw(ji+1,jj,jk+kp) - (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1 - (zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1 411 ! round brackets added to fix the order of floating point operations 412 ! needed to ensure halo 1 - halo 2 compatibility 413 zftu(ji ,jj,jk ) = zftu(ji ,jj,jk ) & 414 & - ( ( zah * zdxt + ( zah_slp - zaei_slp ) * zdzt ) * zbu * ze1ur & 415 & + ( zah * zdxt + zah_slp1 * zdzt_ip1 - zaei_slp1 * zdzt_ip1 ) * zbu * ze1ur & 416 & ) ! bracket for halo 1 - halo 2 compatibility 417 ztfw(ji+1,jj,jk+kp) = ztfw(ji+1,jj,jk+kp) & 418 & - ( (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1 & 419 & + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1 & 420 & ) ! bracket for halo 1 - halo 2 compatibility 421 END_2D 403 422 END DO 404 423 ! 405 DO jp = 0, 1 406 DO kp = 0, 1 407 ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 408 DO_2D( iij, iij-1, iij, iij-1 ) 409 ze2vr = r1_e2v(ji,jj) 410 zdyt = zdjt(ji,jj,jk) * ze2vr 411 ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm) 412 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 413 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 414 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 415 zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 416 ! ln_botmix_triad is .F. mask zah for bottom half cells 417 zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! pahv(ji,jj+jp,jk) ???? 418 zah_slp = zah * zslope_iso 419 IF( ln_ldfeiv ) zaei_slp = aeiv(ji,jj,jk) * zslope_skew ! aeit(ji,jj+jp,jk)*zslope_skew 420 ! NOTE: [halo1-halo2] 421 zftv(ji,jj,jk, jp+1 ) = zftv(ji,jj,jk, jp+1 ) - & 422 & ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 423 ztfw(ji,jj,jk+kp,jp+1,2) = ztfw(ji,jj,jk+kp,jp+1,2) - & 424 & (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 425 END_2D 426 END DO 424 DO kp = 0, 1 425 ! [comm_cleanup] ! DO_2D( 1, 0, 1, 0 ) 426 DO_2D( iij, iij-1, iij, iij-1 ) 427 ze2vr = r1_e2v(ji,jj) 428 zdyt = zdjt(ji,jj,jk) * ze2vr 429 zdyt_jp1 = zdjt(ji,jj+1,jk) * r1_e2v(ji,jj+1) 430 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 431 ze3wr_jp1 = 1._wp / e3w(ji,jj+1,jk+kp,Kmm) 432 zdzt = zdkt3d(ji,jj,kp) * ze3wr 433 zdzt_jp1 = zdkt3d(ji,jj+1,kp) * ze3wr_jp1 434 zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 435 zbv_jp1 = 0.25_wp * e1e2v(ji,jj+1) * e3v(ji,jj+1,jk,Kmm) 436 ! ln_botmix_triad is .F. mask zah for bottom half cells 437 zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! pahv(ji,jj+jp,jk) ???? 438 zah_jp1 = pahv(ji,jj+1,jk) * vmask(ji,jj+1,jk+kp) 439 zah_slp = zah * triadj(ji,jj,jk,1,kp) 440 zah_slp1 = zah * triadj(ji,jj+1,jk,0,kp) 441 zah_slp_jp1 = zah_jp1 * triadj(ji,jj+1,jk,1,kp) 442 IF( ln_ldfeiv ) THEN 443 zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp) 444 zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp) 445 zaei_slp1 = aeiv(ji,jj,jk) * triadj_g(ji,jj+1,jk,0,kp) 446 ENDIF 447 ! zftv(ji,jj ,jk ) = zftv(ji,jj ,jk ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr - ( zah * zdyt + (zah_slp1 - zaei_slp1) * zdzt_jp1 ) * zbv * ze2vr 448 ! ztfw(ji,jj+1,jk+kp) = ztfw(ji,jj+1,jk+kp) - ( zah_slp_jp1 + zaei_slp_jp1) * zdyt_jp1 * zbv_jp1 * ze3wr_jp1 - (zah_slp1 + zaei_slp1) * zdyt * zbv * ze3wr_jp1 449 ! round brackets added to fix the order of floating point operations 450 ! needed to ensure halo 1 - halo 2 compatibility 451 zftv(ji,jj ,jk ) = zftv(ji,jj ,jk ) & 452 & - ( ( zah * zdyt + ( zah_slp - zaei_slp ) * zdzt ) * zbv * ze2vr & 453 & + ( zah * zdyt + zah_slp1 * zdzt_jp1 - zaei_slp1 * zdzt_jp1 ) * zbv * ze2vr & 454 & ) ! bracket for halo 1 - halo 2 compatibility 455 ztfw(ji,jj+1,jk+kp) = ztfw(ji,jj+1,jk+kp) & 456 & - ( ( zah_slp_jp1 + zaei_slp_jp1) * zdyt_jp1 * zbv_jp1 * ze3wr_jp1 & 457 & + ( zah_slp1 + zaei_slp1) * zdyt * zbv * ze3wr_jp1 & 458 & ) ! bracket for halo 1 - halo 2 compatibility 459 END_2D 427 460 END DO 428 461 ENDIF 429 ! NOTE: [halo1-halo2]430 ! ip and jp loop contributions added here to ensure consistent floating point arithmetic for different nn_hls431 DO_2D( iij, iij-1, iij, iij-1 )432 zftu(ji,jj,jk,1) = zftu(ji,jj,jk,1) + zftu(ji,jj,jk,2)433 zftv(ji,jj,jk,1) = zftv(ji,jj,jk,1) + zftv(ji,jj,jk,2)434 END_2D435 DO_2D( iij-1, iij-1, iij-1, iij-1 )436 ztfw(ji,jj,jk,1,1) = (ztfw(ji,jj,jk,1,1) + ztfw(ji-1,jj,jk,2,1)) + &437 & (ztfw(ji,jj,jk,1,2) + ztfw(ji,jj-1,jk,2,2))438 END_2D439 462 ! !== horizontal divergence and add to the general trend ==! 440 463 ! [comm_cleanup] ! DO_2D( 0, 0, 0, 0 ) 441 464 DO_2D( iij-1, iij-1, iij-1, iij-1 ) 442 ! NOTE: [halo1-halo2] 443 ! Extra brackets required to ensure consistent floating point arithmetic for different nn_hls for bilaplacian 444 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) & 445 & + zsign * ( (zftu(ji-1,jj,jk,1) - zftu(ji,jj,jk,1)) & 446 & + (zftv(ji,jj-1,jk,1) - zftv(ji,jj,jk,1)) ) & 447 & / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) 465 ! round brackets added to fix the order of floating point operations 466 ! needed to ensure halo 1 - halo 2 compatibility 467 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) & 468 & + zsign * ( ( zftu(ji-1,jj ,jk) - zftu(ji,jj,jk) & 469 & ) & ! bracket for halo 1 - halo 2 compatibility 470 & + ( zftv(ji,jj-1,jk) - zftv(ji,jj,jk) & 471 & ) & ! bracket for halo 1 - halo 2 compatibility 472 & ) / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) 448 473 END_2D 449 474 ! … … 454 479 ! [comm_cleanup] ! DO_3D( 0, 0, 1, 0, 2, jpkm1 ) 455 480 DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 456 ztfw(ji,jj,jk ,1,1) = ztfw(ji,jj,jk,1,1) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) &481 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) & 457 482 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & 458 483 & * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) … … 463 488 ! [comm_cleanup] ! DO_3D( 0, 0, 1, 0, 2, jpkm1 ) 464 489 DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 465 ztfw(ji,jj,jk ,1,1) = ztfw(ji,jj,jk,1,1) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) &490 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) & 466 491 & * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 467 492 END_3D 468 493 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt and pt2 gradients, resp. 494 ! [comm_cleanup] ! DO_3D( 0, 0, 1, 0, 2, jpkm1 ) 469 495 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 470 ztfw(ji,jj,jk ,1,1) = ztfw(ji,jj,jk,1,1) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) &496 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) & 471 497 & * ( ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) ) & 472 498 & + akz (ji,jj,jk) * ( pt2(ji,jj,jk-1,jn) - pt2(ji,jj,jk,jn) ) ) … … 476 502 ! 477 503 ! [comm_cleanup] ! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Divergence of vertical fluxes added to pta ==! 478 DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 ) 479 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( ztfw(ji,jj,jk+1,1,1) - ztfw(ji,jj,jk,1,1) ) & 504 DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 ) !== Divergence of vertical fluxes added to pta ==! 505 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) & 506 & + zsign * ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) & 480 507 & / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) 481 508 END_3D … … 485 512 ! 486 513 ! ! "Poleward" diffusive heat or salt transports (T-S case only) 487 ! NOTE: [halo1-halo2] 488 IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:,1) ) 514 IF( l_ptr ) CALL dia_ptr_hst( jn, 'ldf', zftv(:,:,:) ) 489 515 ! ! Diffusive heat transports 490 ! NOTE: [halo1-halo2] 491 IF( l_hst ) CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:,1), zftv(:,:,:,1) ) 516 IF( l_hst ) CALL dia_ar5_hst( jn, 'ldf', zftu(:,:,:), zftv(:,:,:) ) 492 517 ! 493 518 ENDIF !== end pass selection ==!
Note: See TracChangeset
for help on using the changeset viewer.
