- Timestamp:
- 2021-12-03T20:32:50+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14318_RK3_stage1
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/TRA/traldf_iso.F90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14318_RK3_stage1
- 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_r14318_RK3_stage1/src/OCE/TRA/traldf_iso.F90
r15512 r15574 132 132 INTEGER :: ji, jj, jk, jn ! dummy loop indices 133 133 INTEGER :: ikt 134 INTEGER :: ierr 134 INTEGER :: ierr, iij ! local integer 135 135 REAL(wp) :: zmsku, zahu_w, zabe1, zcof1, zcoef3 ! local scalars 136 136 REAL(wp) :: zmskv, zahv_w, zabe2, zcof2, zcoef4 ! - - … … 141 141 ! 142 142 IF( kpass == 1 .AND. kt == kit000 ) THEN 143 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile143 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 144 144 IF(lwp) WRITE(numout,*) 145 145 IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype … … 147 147 ENDIF 148 148 ! 149 DO_3D ( 0, 0, 0, 0, 1, jpk )149 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 150 150 akz (ji,jj,jk) = 0._wp 151 151 ah_wslp2(ji,jj,jk) = 0._wp … … 153 153 ENDIF 154 154 ! 155 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile155 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 156 156 l_hst = .FALSE. 157 157 l_ptr = .FALSE. … … 161 161 ENDIF 162 162 ! 163 ! Define pt_rhs halo points for multi-point haloes in bilaplacian case 164 IF( nldf_tra == np_blp_i .AND. kpass == 1 ) THEN ; iij = nn_hls 165 ELSE ; iij = 1 166 ENDIF 167 163 168 ! 164 169 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) … … 172 177 IF( kpass == 1 ) THEN !== first pass only ==! 173 178 ! 174 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )179 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 175 180 ! 176 181 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & … … 179 184 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) 180 185 ! 181 zahu_w = ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 182 & + pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) ) * zmsku 183 zahv_w = ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 184 & + pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) ) * zmskv 186 ! round brackets added to fix the order of floating point operations 187 ! needed to ensure halo 1 - halo 2 compatibility 188 zahu_w = ( ( pahu(ji ,jj,jk-1) + pahu(ji-1,jj,jk) & 189 & ) & ! bracket for halo 1 - halo 2 compatibility 190 & + ( pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) & 191 & ) & ! bracket for halo 1 - halo 2 compatibility 192 & ) * zmsku 193 zahv_w = ( ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 194 & ) & ! bracket for halo 1 - halo 2 compatibility 195 & + ( pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) & 196 & ) & ! bracket for halo 1 - halo 2 compatibility 197 & ) * zmskv 185 198 ! 186 199 ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & … … 189 202 ! 190 203 IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient 191 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 204 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 205 ! round brackets added to fix the order of floating point operations 206 ! needed to ensure halo 1 - halo 2 compatibility 192 207 akz(ji,jj,jk) = 0.25_wp * ( & 193 & 208 & ( ( pahu(ji ,jj,jk) + pahu(ji ,jj,jk-1) ) / ( e1u(ji ,jj) * e1u(ji ,jj) ) & 194 209 & + ( pahu(ji-1,jj,jk) + pahu(ji-1,jj,jk-1) ) / ( e1u(ji-1,jj) * e1u(ji-1,jj) ) & 195 & + ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & 196 & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) ) 210 & ) & ! bracket for halo 1 - halo 2 compatibility 211 & + ( ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & 212 & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) & 213 & ) & ! bracket for halo 1 - halo 2 compatibility 214 & ) 197 215 END_3D 198 216 ! 199 217 IF( ln_traldf_blp ) THEN ! bilaplacian operator 200 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )218 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 201 219 akz(ji,jj,jk) = 16._wp & 202 220 & * ah_wslp2 (ji,jj,jk) & … … 206 224 END_3D 207 225 ELSEIF( ln_traldf_lap ) THEN ! laplacian operator 208 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )226 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 209 227 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 210 228 zcoef0 = rDt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) … … 214 232 ! 215 233 ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 216 DO_3D ( 0, 0, 0, 0, 1, jpk )234 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 217 235 akz(ji,jj,jk) = ah_wslp2(ji,jj,jk) 218 236 END_3D … … 227 245 !! I - masked horizontal derivative 228 246 !!---------------------------------------------------------------------- 229 !!gm : bug.... why (x,:,:)? (1,jpj,:) and (jpi,1,:) should be sufficient.... 230 zdit (ntsi-nn_hls,:,:) = 0._wp ; zdit (ntei+nn_hls,:,:) = 0._wp 231 zdjt (ntsi-nn_hls,:,:) = 0._wp ; zdjt (ntei+nn_hls,:,:) = 0._wp 232 !!end 247 zdit(:,:,:) = 0._wp 248 zdjt(:,:,:) = 0._wp 233 249 234 250 ! Horizontal tracer gradient 235 DO_3D( 1, 0, 1, 0, 1, jpkm1 )251 DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 ) 236 252 zdit(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 237 253 zdjt(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 238 254 END_3D 239 255 IF( ln_zps ) THEN ! botton and surface ocean correction of the horizontal gradient 240 DO_2D( 1, 0, 1, 0 )! bottom correction (partial bottom cell)256 DO_2D( iij, iij-1, iij, iij-1 ) ! bottom correction (partial bottom cell) 241 257 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 242 258 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 243 259 END_2D 244 260 IF( ln_isfcav ) THEN ! first wet level beneath a cavity 245 DO_2D( 1, 0, 1, 0)261 DO_2D( iij, iij-1, iij, iij-1 ) 246 262 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 247 263 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) … … 256 272 DO jk = 1, jpkm1 ! Horizontal slab 257 273 ! 258 DO_2D( 1, 1, 1, 1)274 DO_2D( iij, iij, iij, iij ) 259 275 ! !== Vertical tracer gradient 260 276 zdk1t(ji,jj) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) ! level jk+1 … … 265 281 END_2D 266 282 ! 267 DO_2D( 1, 0, 1, 0) !== Horizontal fluxes283 DO_2D( iij, iij-1, iij, iij-1 ) !== Horizontal fluxes 268 284 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) 269 285 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) … … 278 294 zcof2 = - pahv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv 279 295 ! 280 zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & 281 & + zcof1 * ( zdkt (ji+1,jj) + zdk1t(ji,jj) & 282 & + zdk1t(ji+1,jj) + zdkt (ji,jj) ) ) * umask(ji,jj,jk) 283 zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & 284 & + zcof2 * ( zdkt (ji,jj+1) + zdk1t(ji,jj) & 285 & + zdk1t(ji,jj+1) + zdkt (ji,jj) ) ) * vmask(ji,jj,jk) 296 ! round brackets added to fix the order of floating point operations 297 ! needed to ensure halo 1 - halo 2 compatibility 298 zftu(ji,jj,jk ) = ( zabe1 * zdit(ji,jj,jk) & 299 & + zcof1 * ( ( zdkt (ji+1,jj) + zdk1t(ji,jj) & 300 & ) & ! bracket for halo 1 - halo 2 compatibility 301 & + ( zdk1t(ji+1,jj) + zdkt (ji,jj) & 302 & ) & ! bracket for halo 1 - halo 2 compatibility 303 & ) ) * umask(ji,jj,jk) 304 zftv(ji,jj,jk) = ( zabe2 * zdjt(ji,jj,jk) & 305 & + zcof2 * ( ( zdkt (ji,jj+1) + zdk1t(ji,jj) & 306 & ) & ! bracket for halo 1 - halo 2 compatibility 307 & + ( zdk1t(ji,jj+1) + zdkt (ji,jj) & 308 & ) & ! bracket for halo 1 - halo 2 compatibility 309 & ) ) * vmask(ji,jj,jk) 286 310 END_2D 287 311 ! 288 DO_2D( 0, 0, 0, 0 ) !== horizontal divergence and add to pta 289 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) & 290 & + zsign * ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) + zftv(ji,jj,jk) - zftv(ji,jj-1,jk) ) & 291 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 312 DO_2D( iij-1, iij-1, iij-1, iij-1 ) !== horizontal divergence and add to pta 313 ! round brackets added to fix the order of floating point operations 314 ! needed to ensure halo 1 - halo 2 compatibility 315 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) & 316 & + zsign * ( ( zftu(ji,jj,jk) - zftu(ji-1,jj,jk) & 317 & ) & ! bracket for halo 1 - halo 2 compatibility 318 & + ( zftv(ji,jj,jk) - zftv(ji,jj-1,jk) & 319 & ) & ! bracket for halo 1 - halo 2 compatibility 320 & ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 292 321 END_2D 293 322 END DO ! End of slab … … 302 331 ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp 303 332 304 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! interior (2=<jk=<jpk-1)333 DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) ! interior (2=<jk=<jpk-1) 305 334 ! 306 335 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & … … 317 346 zcoef4 = - zahv_w * e1t(ji,jj) * zmskv * wslpj (ji,jj,jk) 318 347 ! 319 ztfw(ji,jj,jk) = zcoef3 * ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) & 320 & + zdit(ji-1,jj ,jk-1) + zdit(ji ,jj ,jk) ) & 321 & + zcoef4 * ( zdjt(ji ,jj ,jk-1) + zdjt(ji ,jj-1,jk) & 322 & + zdjt(ji ,jj-1,jk-1) + zdjt(ji ,jj ,jk) ) 348 ! round brackets added to fix the order of floating point operations 349 ! needed to ensure halo 1 - halo 2 compatibility 350 ztfw(ji,jj,jk) = zcoef3 * ( ( zdit(ji ,jj ,jk-1) + zdit(ji-1,jj ,jk) & 351 & ) & ! bracket for halo 1 - halo 2 compatibility 352 & + ( zdit(ji-1,jj ,jk-1) + zdit(ji ,jj ,jk) & 353 & ) & ! bracket for halo 1 - halo 2 compatibility 354 & ) & 355 & + zcoef4 * ( ( zdjt(ji ,jj ,jk-1) + zdjt(ji ,jj-1,jk) & 356 & ) & ! bracket for halo 1 - halo 2 compatibility 357 & + ( zdjt(ji ,jj-1,jk-1) + zdjt(ji ,jj ,jk) & 358 & ) & ! bracket for halo 1 - halo 2 compatibility 359 & ) 323 360 END_3D 324 361 ! !== add the vertical 33 flux ==! 325 362 IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz 326 DO_3D( 0, 0, 0, 0, 2, jpkm1 )363 DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 327 364 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) & 328 365 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & … … 333 370 SELECT CASE( kpass ) 334 371 CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2 335 DO_3D( 0, 0, 0, 0, 2, jpkm1 )372 DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 336 373 ztfw(ji,jj,jk) = & 337 374 & ztfw(ji,jj,jk) + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & … … 347 384 ENDIF 348 385 ! 349 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Divergence of vertical fluxes added to pta ==!386 DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 ) !== Divergence of vertical fluxes added to pta ==! 350 387 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) & 351 388 & / e3t(ji,jj,jk,Kmm)
Note: See TracChangeset
for help on using the changeset viewer.
