Changeset 8902 for branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
- Timestamp:
- 2017-12-05T16:58:31+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
r8822 r8902 187 187 !! 188 188 INTEGER :: ji,jj,jk,jn 189 ! VERTICAL REFINEMENT BEGIN190 189 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 191 190 REAL(wp) :: h_in(k1:k2) … … 195 194 REAL(wp) :: zrho_xy 196 195 REAL(wp) :: tabin(k1:k2,n1:n2) 197 ! VERTICAL REFINEMENT END198 196 !!--------------------------------------------- 199 197 ! 200 198 IF (before) THEN 201 199 # if defined key_vertical 202 AGRIF_SpecialValue = -999._wp 203 zrho_xy = Agrif_rhox() * Agrif_rhoy() 204 DO jn = n1,n2-1 205 DO jk=k1,k2 206 DO jj=j1,j2 207 DO ji=i1,i2 208 tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) & 209 * tmask(ji,jj,jk) * zrho_xy + (tmask(ji,jj,jk)-1)*999._wp 210 END DO 211 END DO 212 END DO 213 END DO 200 AGRIF_SpecialValue = -999._wp 201 zrho_xy = Agrif_rhox() * Agrif_rhoy() 202 DO jn = n1,n2-1 214 203 DO jk=k1,k2 215 204 DO jj=j1,j2 216 205 DO ji=i1,i2 217 tabres(ji,jj,jk, n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) * zrho_xy&218 206 tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) & 207 * tmask(ji,jj,jk) * zrho_xy + (tmask(ji,jj,jk)-1)*999._wp 219 208 END DO 220 209 END DO 221 210 END DO 211 END DO 212 DO jk=k1,k2 213 DO jj=j1,j2 214 DO ji=i1,i2 215 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) * zrho_xy & 216 + (tmask(ji,jj,jk)-1)*999._wp 217 END DO 218 END DO 219 END DO 222 220 #else 223 DO jn = n1,n2-1 224 DO jk=k1,k2 225 DO jj=j1,j2 226 DO ji=i1,i2 227 tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 228 END DO 221 DO jn = n1,n2-1 222 DO jk=k1,k2 223 DO jj=j1,j2 224 DO ji=i1,i2 225 tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 229 226 END DO 230 227 END DO 231 END DO 228 END DO 229 END DO 232 230 #endif 233 231 ELSE 234 ! VERTICAL REFINEMENT BEGIN235 232 tabres_child(:,:,:,:) = 0. 236 233 # if defined key_vertical 237 234 AGRIF_SpecialValue = 0._wp 238 235 DO jj=j1,j2 239 DO ji=i1,i2240 N_in = 0241 DO jk=k1,k2 !k2 = jpk of child grid242 IF (tabres(ji,jj,jk,n2) == 0 ) EXIT243 N_in = N_in + 1244 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)245 h_in(N_in) = tabres(ji,jj,jk,n2)/e1e2t(ji,jj)246 ENDDO247 N_out = 0248 DO jk=1,jpk ! jpk of parent grid249 IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF250 N_out = N_out + 1251 h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above252 ENDDO253 IF (N_in > 0) THEN254 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))255 IF (h_diff < -1.e-4) THEN256 print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,e1e2t(ji,jj),sum(h_in(1:N_in)),sum(h_out(1:N_out)), N_in, N_out257 print *,h_in(1:N_in)258 print *,h_out(1:N_out)259 STOP260 ENDIF261 DO jn=n1,n2-1262 CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out)263 ENDDO264 ENDIF265 ENDDO236 DO ji=i1,i2 237 N_in = 0 238 DO jk=k1,k2 !k2 = jpk of child grid 239 IF (tabres(ji,jj,jk,n2) == 0 ) EXIT 240 N_in = N_in + 1 241 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 242 h_in(N_in) = tabres(ji,jj,jk,n2)/e1e2t(ji,jj) 243 ENDDO 244 N_out = 0 245 DO jk=1,jpk ! jpk of parent grid 246 IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 247 N_out = N_out + 1 248 h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 249 ENDDO 250 IF (N_in > 0) THEN !Remove this? 251 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 252 IF (h_diff < -1.e-4) THEN 253 print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 254 print *,h_in(1:N_in) 255 print *,h_out(1:N_out) 256 STOP 257 ENDIF 258 DO jn=n1,n2-1 259 CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 260 ENDDO 261 ENDIF 262 ENDDO 266 263 ENDDO 267 264 #else 268 265 tabres_child(:,:,:,:) = tabres(:,:,:,1:jpts) 269 266 #endif 270 ! WARNING :271 ! tabres replaced by tabres_child in the following272 ! k1 k2 replaced by 1 jpk273 ! VERTICAL REFINEMENT END274 267 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 275 268 ! Add asselin part … … 281 274 tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 282 275 & + atfp * ( tabres_child(ji,jj,jk,jn) & 283 & 276 & - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 284 277 ENDIF 285 278 ENDDO … … 325 318 ! 326 319 IF( before ) THEN 327 print *, i1,i2,j1,j2,k1,k2328 320 zrhoy = Agrif_Rhoy() 329 321 # if defined key_vertical … … 340 332 END DO 341 333 #else 342 343 344 334 DO jk = k1,k2 335 tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 336 END DO 345 337 #endif 346 338 ELSE … … 349 341 AGRIF_SpecialValue = 0._wp 350 342 DO jj=j1,j2 351 DO ji=i1,i2352 N_in = 0353 h_in(:) = 0._wp354 tabin(:) = 0._wp355 DO jk=k1,k2 !k2=jpk of child grid356 IF( tabres(ji,jj,jk,2) < -900) EXIT357 N_in = N_in + 1358 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)359 h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj)360 ENDDO361 N_out = 0362 DO jk=1,jpk363 IF (umask(ji,jj,jk) == 0) EXIT364 N_out = N_out + 1365 h_out(N_out) = e3u_n(ji,jj,jk)366 ENDDO367 IF (N_in * N_out > 0) THEN368 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))369 IF (h_diff < -1.e-4) then343 DO ji=i1,i2 344 N_in = 0 345 h_in(:) = 0._wp 346 tabin(:) = 0._wp 347 DO jk=k1,k2 !k2=jpk of child grid 348 IF( tabres(ji,jj,jk,2) < -900) EXIT 349 N_in = N_in + 1 350 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 351 h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 352 ENDDO 353 N_out = 0 354 DO jk=1,jpk 355 IF (umask(ji,jj,jk) == 0) EXIT 356 N_out = N_out + 1 357 h_out(N_out) = e3u_n(ji,jj,jk) 358 ENDDO 359 IF (N_in * N_out > 0) THEN 360 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 361 IF (h_diff < -1.e-4) THEN 370 362 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 371 363 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 372 excess = 0._wp 373 DO jk=N_in,1,-1 374 thick = MIN(-1*h_diff, h_in(jk)) 375 excess = excess + tabin(jk)*thick*e2u(ji,jj) 376 tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 377 h_diff = h_diff + thick 378 IF ( h_diff == 0) THEN 379 N_in = jk 380 h_in(jk) = h_in(jk) - thick 381 EXIT 382 ENDIF 383 ENDDO 384 ENDIF 385 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 386 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 387 ENDIF 364 excess = 0._wp 365 DO jk=N_in,1,-1 366 thick = MIN(-1*h_diff, h_in(jk)) 367 excess = excess + tabin(jk)*thick*e2u(ji,jj) 368 tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 369 h_diff = h_diff + thick 370 IF ( h_diff == 0) THEN 371 N_in = jk 372 h_in(jk) = h_in(jk) - thick 373 EXIT 374 ENDIF 375 ENDDO 376 ENDIF 377 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 378 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 379 ENDIF 380 ENDDO 388 381 ENDDO 389 ENDDO390 382 #else 391 383 DO jk=1,jpk … … 463 455 AGRIF_SpecialValue = 0._wp 464 456 DO jj=j1,j2 465 DO ji=i1,i2466 N_in = 0467 DO jk=k1,k2468 IF (tabres(ji,jj,jk,2) < -900) EXIT469 N_in = N_in + 1470 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2)471 h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj)472 ENDDO473 N_out = 0474 DO jk=1,jpk475 IF (vmask(ji,jj,jk) == 0) EXIT476 N_out = N_out + 1477 h_out(N_out) = e3v_n(ji,jj,jk)478 ENDDO479 IF (N_in * N_out > 0) THEN480 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))481 IF (h_diff < -1.e-4) then457 DO ji=i1,i2 458 N_in = 0 459 DO jk=k1,k2 460 IF (tabres(ji,jj,jk,2) < -900) EXIT 461 N_in = N_in + 1 462 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 463 h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) 464 ENDDO 465 N_out = 0 466 DO jk=1,jpk 467 IF (vmask(ji,jj,jk) == 0) EXIT 468 N_out = N_out + 1 469 h_out(N_out) = e3v_n(ji,jj,jk) 470 ENDDO 471 IF (N_in * N_out > 0) THEN 472 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 473 IF (h_diff < -1.e-4) then 482 474 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 483 475 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 484 excess = 0._wp485 DO jk=N_in,1,-1486 thick = MIN(-1*h_diff, h_in(jk))487 excess = excess + tabin(jk)*thick*e2u(ji,jj)488 tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))489 h_diff = h_diff + thick490 IF ( h_diff == 0) THEN491 N_in = jk492 h_in(jk) = h_in(jk) - thick493 EXIT494 ENDIF495 ENDDO496 ENDIF497 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)498 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))499 ENDIF500 ENDDO476 excess = 0._wp 477 DO jk=N_in,1,-1 478 thick = MIN(-1*h_diff, h_in(jk)) 479 excess = excess + tabin(jk)*thick*e2u(ji,jj) 480 tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 481 h_diff = h_diff + thick 482 IF ( h_diff == 0) THEN 483 N_in = jk 484 h_in(jk) = h_in(jk) - thick 485 EXIT 486 ENDIF 487 ENDDO 488 ENDIF 489 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 490 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 491 ENDIF 492 ENDDO 501 493 ENDDO 502 494 #else 503 504 505 506 507 508 509 495 DO jk=k1,k2 496 DO jj=j1,j2 497 DO ji=i1,i2 498 tabres(ji,jj,jk,1) = e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) 499 END DO 500 END DO 501 END DO 510 502 #endif 511 503
Note: See TracChangeset
for help on using the changeset viewer.
