- Timestamp:
- 2021-07-16T20:00:12+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/TRA/traldf_triad.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO
- Property svn:externals
-
old new 9 9 10 10 # SETTE 11 ^/utils/CI/sette@14244 sette 11 ^/utils/CI/sette@HEAD sette 12
-
- Property svn:externals
-
NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/OCE/TRA/traldf_triad.F90
r14215 r15127 13 13 USE oce ! ocean dynamics and active tracers 14 14 USE dom_oce ! ocean space and time domain 15 ! TEMP: [tiling] This change not necessary if XIOS has subdomain support16 USE domtile17 15 USE domutl, ONLY : is_tile 18 16 USE phycst ! physical constants … … 109 107 REAL(wp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) :: pt_rhs ! tracer trend 110 108 ! 111 INTEGER :: ji, jj, jk, jn ! dummy loop indices 112 INTEGER :: ip,jp,kp ! dummy loop indices 113 INTEGER :: ierr ! local integer 114 REAL(wp) :: zmsku, zabe1, zcof1, zcoef3 ! local scalars 115 REAL(wp) :: zmskv, zabe2, zcof2, zcoef4 ! - - 109 INTEGER :: ji, jj, jk, jn, kp, iij ! dummy loop indices 116 110 REAL(wp) :: zcoef0, ze3w_2, zsign ! - - 117 111 ! 118 REAL(wp) :: zslope_skew, zslope_iso, zslope2, zbu, zbv 119 REAL(wp) :: ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt 120 REAL(wp) :: zah, zah_slp, zaei_slp 121 REAL(wp), DIMENSION(A2D(nn_hls),0:1) :: zdkt3d ! vertical tracer gradient at 2 levels 122 REAL(wp), DIMENSION(A2D(nn_hls) ) :: z2d ! 2D workspace 123 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk) :: zdit, zdjt, zftu, zftv, ztfw ! 3D - 124 ! TEMP: [tiling] This can be A2D(nn_hls) if XIOS has subdomain support 125 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 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 - 126 118 !!---------------------------------------------------------------------- 127 119 ! 128 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile120 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 129 121 IF( kpass == 1 .AND. kt == kit000 ) THEN 130 122 IF(lwp) WRITE(numout,*) … … 142 134 ENDIF 143 135 ! 136 ! Define pt_rhs halo points for multi-point haloes in bilaplacian case 137 IF( nldf_tra == np_blp_it .AND. kpass == 1 ) THEN ; iij = nn_hls 138 ELSE ; iij = 1 139 ENDIF 140 141 ! 144 142 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) 145 143 ELSE ; zsign = -1._wp … … 152 150 IF( kpass == 1 ) THEN !== first pass only and whatever the tracer is ==! 153 151 ! 154 DO_3D ( 0, 0, 0, 0, 1, jpk )152 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 155 153 akz (ji,jj,jk) = 0._wp 156 154 ah_wslp2(ji,jj,jk) = 0._wp 157 155 END_3D 158 156 ! 159 DO ip = 0, 1 ! i-k triads 160 DO kp = 0, 1 161 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 162 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 163 zbu = e1e2u(ji-ip,jj) * e3u(ji-ip,jj,jk,Kmm) 164 zah = 0.25_wp * pahu(ji-ip,jj,jk) 165 zslope_skew = triadi_g(ji,jj,jk,1-ip,kp) 166 ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 167 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) 168 zslope2 = zslope2 *zslope2 169 ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2 170 akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + zah * r1_e1u(ji-ip,jj) & 171 & * r1_e1u(ji-ip,jj) * umask(ji-ip,jj,jk+kp) 172 ! 173 END_3D 174 END DO 157 DO kp = 0, 1 ! i-k triads 158 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 159 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 160 zbu = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 161 zbu1 = e1e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) 162 zah = 0.25_wp * pahu(ji,jj,jk) 163 zah1 = 0.25_wp * pahu(ji-1,jj,jk) 164 ! Subtract s-coordinate slope at t-points to give slope rel to s-surfaces (do this by *adding* gradient of depth) 165 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) 166 zslope2 = zslope2 *zslope2 167 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) 168 zslope21 = zslope21 *zslope21 169 ! round brackets added to fix the order of floating point operations 170 ! needed to ensure halo 1 - halo 2 compatibility 171 ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + ( zah * zbu * ze3wr * r1_e1e2t(ji,jj) * zslope2 & 172 & + zah1 * zbu1 * ze3wr * r1_e1e2t(ji,jj) * zslope21 & 173 & ) ! bracket for halo 1 - halo 2 compatibility 174 akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + ( zah * r1_e1u(ji,jj) * r1_e1u(ji,jj) * umask(ji,jj,jk+kp) & 175 + zah1 * r1_e1u(ji-1,jj) * r1_e1u(ji-1,jj) * umask(ji-1,jj,jk+kp) & 176 & ) ! bracket for halo 1 - halo 2 compatibility 177 END_3D 175 178 END DO 176 179 ! 177 DO jp = 0, 1 ! j-k triads 178 DO kp = 0, 1 179 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 180 ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm) 181 zbv = e1e2v(ji,jj-jp) * e3v(ji,jj-jp,jk,Kmm) 182 zah = 0.25_wp * pahv(ji,jj-jp,jk) 183 zslope_skew = triadj_g(ji,jj,jk,1-jp,kp) 184 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 185 ! (do this by *adding* gradient of depth) 186 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) 187 zslope2 = zslope2 * zslope2 188 ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2 189 akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + zah * r1_e2v(ji,jj-jp) & 190 & * r1_e2v(ji,jj-jp) * vmask(ji,jj-jp,jk+kp) 191 ! 192 END_3D 193 END DO 180 DO kp = 0, 1 ! j-k triads 181 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 182 ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm) 183 zbv = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 184 zbv1 = e1e2v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) 185 zah = 0.25_wp * pahv(ji,jj,jk) 186 zah1 = 0.25_wp * pahv(ji,jj-1,jk) 187 ! Subtract s-coordinate slope at t-points to give slope rel to s surfaces 188 ! (do this by *adding* gradient of depth) 189 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) 190 zslope2 = zslope2 * zslope2 191 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) 192 zslope21 = zslope21 * zslope21 193 ! round brackets added to fix the order of floating point operations 194 ! needed to ensure halo 1 - halo 2 compatibility 195 ah_wslp2(ji,jj,jk+kp) = ah_wslp2(ji,jj,jk+kp) + ( zah * zbv * ze3wr * r1_e1e2t(ji,jj) * zslope2 & 196 & + zah1 * zbv1 * ze3wr * r1_e1e2t(ji,jj) * zslope21 & 197 & ) ! bracket for halo 1 - halo 2 compatibility 198 akz (ji,jj,jk+kp) = akz (ji,jj,jk+kp) + ( zah * r1_e2v(ji,jj) * r1_e2v(ji,jj) * vmask(ji,jj,jk+kp) & 199 & + zah1 * r1_e2v(ji,jj-1) * r1_e2v(ji,jj-1) * vmask(ji,jj-1,jk+kp) & 200 & ) ! bracket for halo 1 - halo 2 compatibility 201 END_3D 194 202 END DO 195 203 ! … … 197 205 ! 198 206 IF( ln_traldf_blp ) THEN ! bilaplacian operator 199 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )207 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 200 208 akz(ji,jj,jk) = 16._wp & 201 209 & * ah_wslp2 (ji,jj,jk) & 202 210 & * ( akz (ji,jj,jk) & 203 211 & + ah_wslp2(ji,jj,jk) & 204 & / ( e3w 212 & / ( e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) ) ) 205 213 END_3D 206 214 ELSEIF( ln_traldf_lap ) THEN ! laplacian operator 207 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )215 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 208 216 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 209 217 zcoef0 = rDt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) … … 213 221 ! 214 222 ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 215 DO_3D ( 0, 0, 0, 0, 1, jpk )223 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 216 224 akz(ji,jj,jk) = ah_wslp2(ji,jj,jk) 217 225 END_3D 218 226 ENDIF 219 227 ! 220 ! TEMP: [tiling] These changes not necessary if XIOS has subdomain support 221 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain 222 IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN 223 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) 224 225 zpsi_uw(:,:,:) = 0._wp 226 zpsi_vw(:,:,:) = 0._wp 227 228 DO jp = 0, 1 229 DO kp = 0, 1 230 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 231 zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 232 & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+jp,jj,jk,1-jp,kp) 233 zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) & 234 & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+jp,jk,1-jp,kp) 235 END_3D 236 END DO 237 END DO 238 CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 239 240 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile ) 241 ENDIF 228 IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN 229 zpsi_uw(:,:,:) = 0._wp 230 zpsi_vw(:,:,:) = 0._wp 231 232 DO kp = 0, 1 233 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 234 ! round brackets added to fix the order of floating point operations 235 ! needed to ensure halo 1 - halo 2 compatibility 236 zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 237 & + ( 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji,jj,jk,1,kp) & 238 & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+1,jj,jk,0,kp) & 239 & ) ! bracket for halo 1 - halo 2 compatibility 240 zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) & 241 & + ( 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj,jk,1,kp) & 242 & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+1,jk,0,kp) & 243 & ) ! bracket for halo 1 - halo 2 compatibility 244 END_3D 245 END DO 246 CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 242 247 ENDIF 243 248 ! … … 252 257 zftu(:,:,:) = 0._wp 253 258 zftv(:,:,:) = 0._wp 254 ! 255 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) !== before lateral T & S gradients at T-level jk ==! 259 zdit(:,:,:) = 0._wp 260 zdjt(:,:,:) = 0._wp 261 ! 262 DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 ) !== before lateral T & S gradients at T-level jk ==! 256 263 zdit(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 257 264 zdjt(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 258 265 END_3D 259 266 IF( ln_zps .AND. l_grad_zps ) THEN ! partial steps: correction at top/bottom ocean level 260 DO_2D( 1, 0, 1, 0) ! bottom level267 DO_2D( iij, iij-1, iij, iij-1 ) ! bottom level 261 268 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 262 269 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 263 270 END_2D 264 271 IF( ln_isfcav ) THEN ! top level (ocean cavities only) 265 DO_2D( 1, 0, 1, 0)272 DO_2D( iij, iij-1, iij, iij-1 ) 266 273 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn) 267 274 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn) … … 276 283 DO jk = 1, jpkm1 277 284 ! !== Vertical tracer gradient at level jk and jk+1 278 DO_2D( 1, 1, 1, 1)285 DO_2D( iij, iij, iij, iij ) 279 286 zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 280 287 END_2D … … 283 290 IF( jk == 1 ) THEN ; zdkt3d(:,:,0) = zdkt3d(:,:,1) 284 291 ELSE 285 DO_2D( 1, 1, 1, 1)292 DO_2D( iij, iij, iij, iij ) 286 293 zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 287 294 END_2D … … 289 296 ! 290 297 zaei_slp = 0._wp 298 zaei_slp_ip1 = 0._wp 299 zaei_slp_jp1 = 0._wp 300 zaei_slp1 = 0._wp 291 301 ! 292 302 IF( ln_botmix_triad ) THEN 293 DO ip = 0, 1 !== Horizontal & vertical fluxes 294 DO kp = 0, 1 295 DO_2D( 1, 0, 1, 0 ) 296 ze1ur = r1_e1u(ji,jj) 297 zdxt = zdit(ji,jj,jk) * ze1ur 298 ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm) 299 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 300 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 301 zslope_iso = triadi (ji+ip,jj,jk,1-ip,kp) 302 ! 303 zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 304 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked.... 305 zah = pahu(ji,jj,jk) 306 zah_slp = zah * zslope_iso 307 IF( ln_ldfeiv ) zaei_slp = aeiu(ji,jj,jk) * zslope_skew 308 zftu(ji ,jj,jk ) = zftu(ji ,jj,jk ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 309 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - ( zah_slp + zaei_slp) * zdxt * zbu * ze3wr 310 END_2D 311 END DO 303 DO kp = 0, 1 !== Horizontal & vertical fluxes 304 DO_2D( iij, iij-1, iij, iij-1 ) 305 ze1ur = r1_e1u(ji,jj) 306 zdxt = zdit(ji,jj,jk) * ze1ur 307 zdxt_ip1 = zdit(ji+1,jj,jk) * r1_e1u(ji+1,jj) 308 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 309 ze3wr_ip1 = 1._wp / e3w(ji+1,jj,jk+kp,Kmm) 310 zdzt = zdkt3d(ji,jj,kp) * ze3wr 311 zdzt_ip1 = zdkt3d(ji+1,jj,kp) * ze3wr_ip1 312 ! 313 zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 314 zbu_ip1 = 0.25_wp * e1e2u(ji+1,jj) * e3u(ji+1,jj,jk,Kmm) 315 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked.... 316 zah = pahu(ji,jj,jk) 317 zah_ip1 = pahu(ji+1,jj,jk) 318 zah_slp = zah * triadi(ji,jj,jk,1,kp) 319 zah_slp_ip1 = zah_ip1 * triadi(ji+1,jj,jk,1,kp) 320 zah_slp1 = zah * triadi(ji+1,jj,jk,0,kp) 321 IF( ln_ldfeiv ) THEN 322 zaei_slp = aeiu(ji,jj,jk) * triadi_g(ji,jj,jk,1,kp) 323 zaei_slp_ip1 = aeiu(ji+1,jj,jk) * triadi_g(ji+1,jj,jk,1,kp) 324 zaei_slp1 = aeiu(ji,jj,jk) * triadi_g(ji+1,jj,jk,0,kp) 325 ENDIF 326 ! round brackets added to fix the order of floating point operations 327 ! needed to ensure halo 1 - halo 2 compatibility 328 zftu(ji ,jj,jk ) = zftu(ji ,jj,jk ) & 329 & - ( ( zah * zdxt + ( zah_slp - zaei_slp ) * zdzt ) * zbu * ze1ur & 330 & + ( zah * zdxt + zah_slp1 * zdzt_ip1 - zaei_slp1 * zdzt_ip1 ) * zbu * ze1ur & 331 & ) ! bracket for halo 1 - halo 2 compatibility 332 ztfw(ji+1,jj,jk+kp) = ztfw(ji+1,jj,jk+kp) & 333 & - ( (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1 & 334 & + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1 & 335 & ) ! bracket for halo 1 - halo 2 compatibility 336 END_2D 312 337 END DO 313 338 ! 314 DO jp = 0, 1 315 DO kp = 0, 1 316 DO_2D( 1, 0, 1, 0 ) 317 ze2vr = r1_e2v(ji,jj) 318 zdyt = zdjt(ji,jj,jk) * ze2vr 319 ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm) 320 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 321 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 322 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 323 zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 324 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahv is masked... 325 zah = pahv(ji,jj,jk) 326 zah_slp = zah * zslope_iso 327 IF( ln_ldfeiv ) zaei_slp = aeiv(ji,jj,jk) * zslope_skew 328 zftv(ji,jj ,jk ) = zftv(ji,jj ,jk ) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 329 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - ( zah_slp + zaei_slp ) * zdyt * zbv * ze3wr 330 END_2D 331 END DO 339 DO kp = 0, 1 340 DO_2D( iij, iij-1, iij, iij-1 ) 341 ze2vr = r1_e2v(ji,jj) 342 zdyt = zdjt(ji,jj,jk) * ze2vr 343 zdyt_jp1 = zdjt(ji,jj+1,jk) * r1_e2v(ji,jj+1) 344 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 345 ze3wr_jp1 = 1._wp / e3w(ji,jj+1,jk+kp,Kmm) 346 zdzt = zdkt3d(ji,jj,kp) * ze3wr 347 zdzt_jp1 = zdkt3d(ji,jj+1,kp) * ze3wr_jp1 348 zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 349 zbv_jp1 = 0.25_wp * e1e2v(ji,jj+1) * e3v(ji,jj+1,jk,Kmm) 350 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked.... 351 zah = pahv(ji,jj,jk) ! pahv(ji,jj+jp,jk) ???? 352 zah_jp1 = pahv(ji,jj+1,jk) 353 zah_slp = zah * triadj(ji,jj,jk,1,kp) 354 zah_slp1 = zah * triadj(ji,jj+1,jk,0,kp) 355 zah_slp_jp1 = zah_jp1 * triadj(ji,jj+1,jk,1,kp) 356 IF( ln_ldfeiv ) THEN 357 zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp) 358 zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp) 359 zaei_slp1 = aeiv(ji,jj,jk) * triadj_g(ji,jj+1,jk,0,kp) 360 ENDIF 361 ! round brackets added to fix the order of floating point operations 362 ! needed to ensure halo 1 - halo 2 compatibility 363 zftv(ji,jj ,jk ) = zftv(ji,jj ,jk ) & 364 & - ( ( zah * zdyt + ( zah_slp - zaei_slp ) * zdzt ) * zbv * ze2vr & 365 & + ( zah * zdyt + zah_slp1 * zdzt_jp1 - zaei_slp1 * zdzt_jp1 ) * zbv * ze2vr & 366 & ) ! bracket for halo 1 - halo 2 compatibility 367 ztfw(ji,jj+1,jk+kp) = ztfw(ji,jj+1,jk+kp) & 368 & - ( ( zah_slp_jp1 + zaei_slp_jp1) * zdyt_jp1 * zbv_jp1 * ze3wr_jp1 & 369 & + ( zah_slp1 + zaei_slp1) * zdyt * zbv * ze3wr_jp1 & 370 & ) ! bracket for halo 1 - halo 2 compatibility 371 END_2D 332 372 END DO 333 373 ! 334 374 ELSE 335 375 ! 336 DO ip = 0, 1 !== Horizontal & vertical fluxes 337 DO kp = 0, 1 338 DO_2D( 1, 0, 1, 0 ) 339 ze1ur = r1_e1u(ji,jj) 340 zdxt = zdit(ji,jj,jk) * ze1ur 341 ze3wr = 1._wp / e3w(ji+ip,jj,jk+kp,Kmm) 342 zdzt = zdkt3d(ji+ip,jj,kp) * ze3wr 343 zslope_skew = triadi_g(ji+ip,jj,jk,1-ip,kp) 344 zslope_iso = triadi(ji+ip,jj,jk,1-ip,kp) 345 ! 346 zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 347 ! ln_botmix_triad is .F. mask zah for bottom half cells 348 zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp) ! pahu(ji+ip,jj,jk) ===>> ???? 349 zah_slp = zah * zslope_iso 350 IF( ln_ldfeiv ) zaei_slp = aeiu(ji,jj,jk) * zslope_skew ! aeit(ji+ip,jj,jk)*zslope_skew 351 zftu(ji ,jj,jk ) = zftu(ji ,jj,jk ) - ( zah * zdxt + (zah_slp - zaei_slp) * zdzt ) * zbu * ze1ur 352 ztfw(ji+ip,jj,jk+kp) = ztfw(ji+ip,jj,jk+kp) - (zah_slp + zaei_slp) * zdxt * zbu * ze3wr 353 END_2D 354 END DO 376 DO kp = 0, 1 !== Horizontal & vertical fluxes 377 DO_2D( iij, iij-1, iij, iij-1 ) 378 ze1ur = r1_e1u(ji,jj) 379 zdxt = zdit(ji,jj,jk) * ze1ur 380 zdxt_ip1 = zdit(ji+1,jj,jk) * r1_e1u(ji+1,jj) 381 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 382 ze3wr_ip1 = 1._wp / e3w(ji+1,jj,jk+kp,Kmm) 383 zdzt = zdkt3d(ji,jj,kp) * ze3wr 384 zdzt_ip1 = zdkt3d(ji+1,jj,kp) * ze3wr_ip1 385 ! 386 zbu = 0.25_wp * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) 387 zbu_ip1 = 0.25_wp * e1e2u(ji+1,jj) * e3u(ji+1,jj,jk,Kmm) 388 ! ln_botmix_triad is .F. mask zah for bottom half cells 389 zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp) ! pahu(ji+ip,jj,jk) ===>> ???? 390 zah_ip1 = pahu(ji+1,jj,jk) * umask(ji+1,jj,jk+kp) 391 zah_slp = zah * triadi(ji,jj,jk,1,kp) 392 zah_slp_ip1 = zah_ip1 * triadi(ji+1,jj,jk,1,kp) 393 zah_slp1 = zah * triadi(ji+1,jj,jk,0,kp) 394 IF( ln_ldfeiv ) THEN 395 zaei_slp = aeiu(ji,jj,jk) * triadi_g(ji,jj,jk,1,kp) 396 zaei_slp_ip1 = aeiu(ji+1,jj,jk) * triadi_g(ji+1,jj,jk,1,kp) 397 zaei_slp1 = aeiu(ji,jj,jk) * triadi_g(ji+1,jj,jk,0,kp) 398 ENDIF 399 ! round brackets added to fix the order of floating point operations 400 ! needed to ensure halo 1 - halo 2 compatibility 401 zftu(ji ,jj,jk ) = zftu(ji ,jj,jk ) & 402 & - ( ( zah * zdxt + ( zah_slp - zaei_slp ) * zdzt ) * zbu * ze1ur & 403 & + ( zah * zdxt + zah_slp1 * zdzt_ip1 - zaei_slp1 * zdzt_ip1 ) * zbu * ze1ur & 404 & ) ! bracket for halo 1 - halo 2 compatibility 405 ztfw(ji+1,jj,jk+kp) = ztfw(ji+1,jj,jk+kp) & 406 & - ( (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1 & 407 & + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1 & 408 & ) ! bracket for halo 1 - halo 2 compatibility 409 END_2D 355 410 END DO 356 411 ! 357 DO jp = 0, 1 358 DO kp = 0, 1 359 DO_2D( 1, 0, 1, 0 ) 360 ze2vr = r1_e2v(ji,jj) 361 zdyt = zdjt(ji,jj,jk) * ze2vr 362 ze3wr = 1._wp / e3w(ji,jj+jp,jk+kp,Kmm) 363 zdzt = zdkt3d(ji,jj+jp,kp) * ze3wr 364 zslope_skew = triadj_g(ji,jj+jp,jk,1-jp,kp) 365 zslope_iso = triadj(ji,jj+jp,jk,1-jp,kp) 366 zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 367 ! ln_botmix_triad is .F. mask zah for bottom half cells 368 zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! pahv(ji,jj+jp,jk) ???? 369 zah_slp = zah * zslope_iso 370 IF( ln_ldfeiv ) zaei_slp = aeiv(ji,jj,jk) * zslope_skew ! aeit(ji,jj+jp,jk)*zslope_skew 371 zftv(ji,jj,jk) = zftv(ji,jj,jk) - ( zah * zdyt + (zah_slp - zaei_slp) * zdzt ) * zbv * ze2vr 372 ztfw(ji,jj+jp,jk+kp) = ztfw(ji,jj+jp,jk+kp) - (zah_slp + zaei_slp) * zdyt * zbv * ze3wr 373 END_2D 374 END DO 412 DO kp = 0, 1 413 DO_2D( iij, iij-1, iij, iij-1 ) 414 ze2vr = r1_e2v(ji,jj) 415 zdyt = zdjt(ji,jj,jk) * ze2vr 416 zdyt_jp1 = zdjt(ji,jj+1,jk) * r1_e2v(ji,jj+1) 417 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 418 ze3wr_jp1 = 1._wp / e3w(ji,jj+1,jk+kp,Kmm) 419 zdzt = zdkt3d(ji,jj,kp) * ze3wr 420 zdzt_jp1 = zdkt3d(ji,jj+1,kp) * ze3wr_jp1 421 zbv = 0.25_wp * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) 422 zbv_jp1 = 0.25_wp * e1e2v(ji,jj+1) * e3v(ji,jj+1,jk,Kmm) 423 ! ln_botmix_triad is .F. mask zah for bottom half cells 424 zah = pahv(ji,jj,jk) * vmask(ji,jj,jk+kp) ! pahv(ji,jj+jp,jk) ???? 425 zah_jp1 = pahv(ji,jj+1,jk) * vmask(ji,jj+1,jk+kp) 426 zah_slp = zah * triadj(ji,jj,jk,1,kp) 427 zah_slp1 = zah * triadj(ji,jj+1,jk,0,kp) 428 zah_slp_jp1 = zah_jp1 * triadj(ji,jj+1,jk,1,kp) 429 IF( ln_ldfeiv ) THEN 430 zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp) 431 zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp) 432 zaei_slp1 = aeiv(ji,jj,jk) * triadj_g(ji,jj+1,jk,0,kp) 433 ENDIF 434 ! round brackets added to fix the order of floating point operations 435 ! needed to ensure halo 1 - halo 2 compatibility 436 zftv(ji,jj ,jk ) = zftv(ji,jj ,jk ) & 437 & - ( ( zah * zdyt + ( zah_slp - zaei_slp ) * zdzt ) * zbv * ze2vr & 438 & + ( zah * zdyt + zah_slp1 * zdzt_jp1 - zaei_slp1 * zdzt_jp1 ) * zbv * ze2vr & 439 & ) ! bracket for halo 1 - halo 2 compatibility 440 ztfw(ji,jj+1,jk+kp) = ztfw(ji,jj+1,jk+kp) & 441 & - ( ( zah_slp_jp1 + zaei_slp_jp1) * zdyt_jp1 * zbv_jp1 * ze3wr_jp1 & 442 & + ( zah_slp1 + zaei_slp1) * zdyt * zbv * ze3wr_jp1 & 443 & ) ! bracket for halo 1 - halo 2 compatibility 444 END_2D 375 445 END DO 376 446 ENDIF 377 447 ! !== horizontal divergence and add to the general trend ==! 378 DO_2D( 0, 0, 0, 0 ) 379 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) & 380 & + zsign * ( zftu(ji-1,jj ,jk) - zftu(ji,jj,jk) & 381 & + zftv(ji,jj-1,jk) - zftv(ji,jj,jk) ) & 382 & / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) 448 DO_2D( iij-1, iij-1, iij-1, iij-1 ) 449 ! round brackets added to fix the order of floating point operations 450 ! needed to ensure halo 1 - halo 2 compatibility 451 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) & 452 & + zsign * ( ( zftu(ji-1,jj ,jk) - zftu(ji,jj,jk) & 453 & ) & ! bracket for halo 1 - halo 2 compatibility 454 & + ( zftv(ji,jj-1,jk) - zftv(ji,jj,jk) & 455 & ) & ! bracket for halo 1 - halo 2 compatibility 456 & ) / ( e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) ) 383 457 END_2D 384 458 ! … … 387 461 ! !== add the vertical 33 flux ==! 388 462 IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz 389 DO_3D( 0, 0, 1, 0, 2, jpkm1 )463 DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 390 464 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) & 391 465 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & … … 395 469 SELECT CASE( kpass ) 396 470 CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2 397 DO_3D( 0, 0, 1, 0, 2, jpkm1 )471 DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 398 472 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) & 399 473 & * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 400 474 END_3D 401 475 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt and pt2 gradients, resp. 402 DO_3D( 0, 0, 1, 0, 2, jpkm1 )476 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 403 477 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) & 404 478 & * ( ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) ) & … … 408 482 ENDIF 409 483 ! 410 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Divergence of vertical fluxes added to pta ==!484 DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 ) !== Divergence of vertical fluxes added to pta ==! 411 485 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) & 412 486 & + zsign * ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) &
Note: See TracChangeset
for help on using the changeset viewer.
