Changeset 14086 for NEMO/trunk/src/NST/agrif_oce_update.F90
- Timestamp:
- 2020-12-04T12:37:14+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/NST/agrif_oce_update.F90
r14064 r14086 1 #undef DECAL_FEEDBACK 1 #undef DECAL_FEEDBACK /* SEPARATION of INTERFACES */ 2 2 #undef DECAL_FEEDBACK_2D /* SEPARATION of INTERFACES (Barotropic mode) */ 3 3 #undef VOL_REFLUX /* VOLUME REFLUXING*/ … … 21 21 USE zdf_oce ! vertical physics: ocean variables 22 22 USE agrif_oce 23 USE dom_oce 23 24 ! 24 25 USE in_out_manager ! I/O manager … … 34 35 35 36 PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Agrif_Update_vvl, Agrif_Update_ssh 36 PUBLIC Update_Scales 37 PUBLIC Update_Scales, Agrif_Check_parent_bat 37 38 38 39 !! * Substitutions … … 54 55 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers from grid Number',Agrif_Fixed() 55 56 56 #if defined key_vertical 57 ! Effect of this has to be carrefully checked 58 ! depending on what the nesting tools ensure for 59 ! volume conservation: 57 l_vremap = ln_vert_remap 58 Agrif_UseSpecialValueInUpdate = .NOT.l_vremap 59 Agrif_SpecialValueFineGrid = 0._wp 60 ! 61 # if ! defined DECAL_FEEDBACK 62 CALL Agrif_Update_Variable(ts_update_id, procname=updateTS) 63 ! near boundary update: 64 ! CALL Agrif_Update_Variable(ts_update_id,locupdate=(/0,2/), procname=updateTS) 65 # else 66 CALL Agrif_Update_Variable(ts_update_id, locupdate=(/1,0/),procname=updateTS) 67 ! near boundary update: 68 ! CALL Agrif_Update_Variable(ts_update_id,locupdate=(/1,2/), procname=updateTS) 69 # endif 70 ! 60 71 Agrif_UseSpecialValueInUpdate = .FALSE. 61 #else 62 Agrif_UseSpecialValueInUpdate = .TRUE. 63 #endif 64 Agrif_SpecialValueFineGrid = 0._wp 65 ! 66 # if ! defined DECAL_FEEDBACK 67 CALL Agrif_Update_Variable(tsn_id, procname=updateTS) 68 ! near boundary update: 69 ! CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS) 70 # else 71 CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS) 72 ! near boundary update: 73 ! CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS) 74 # endif 75 ! 76 Agrif_UseSpecialValueInUpdate = .FALSE. 72 l_vremap = .FALSE. 77 73 ! 78 74 ! … … 90 86 Agrif_UseSpecialValueInUpdate = .FALSE. 91 87 Agrif_SpecialValueFineGrid = 0._wp 92 93 use_sign_north = .TRUE. 94 sign_north = -1._wp 95 96 ! 88 l_vremap = ln_vert_remap 89 use_sign_north = .TRUE. 90 sign_north = -1._wp 91 ! 92 # if ! defined DECAL_FEEDBACK_2D 93 CALL Agrif_Update_Variable(unb_update_id,locupdate1=(/ nn_shift_bar,-2/),locupdate2=(/ nn_shift_bar,-2/),procname = updateU2d) 94 CALL Agrif_Update_Variable(vnb_update_id,locupdate1=(/ nn_shift_bar,-2/),locupdate2=(/ nn_shift_bar,-2/),procname = updateV2d) 95 # else 96 CALL Agrif_Update_Variable(unb_update_id,locupdate1=(/ nn_shift_bar,-2/),locupdate2=(/1+nn_shift_bar,-2/),procname = updateU2d) 97 CALL Agrif_Update_Variable(vnb_update_id,locupdate1=(/1+nn_shift_bar,-2/),locupdate2=(/ nn_shift_bar,-2/),procname = updateV2d) 98 # endif 99 ! 100 IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 101 ! Update time integrated transports 102 # if ! defined DECAL_FEEDBACK_2D 103 CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/ nn_shift_bar,-2/),locupdate2=(/ nn_shift_bar,-2/),procname = updateub2b) 104 CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ nn_shift_bar,-2/),locupdate2=(/ nn_shift_bar,-2/),procname = updatevb2b) 105 # else 106 CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/ nn_shift_bar,-2/),locupdate2=(/1+nn_shift_bar,-2/),procname = updateub2b) 107 CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1+nn_shift_bar,-2/),locupdate2=(/ nn_shift_bar,-2/),procname = updatevb2b) 108 # endif 109 END IF 110 97 111 # if ! defined DECAL_FEEDBACK 98 112 CALL Agrif_Update_Variable(un_update_id,procname = updateU) … … 108 122 ! CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV) 109 123 # endif 110 111 # if ! defined DECAL_FEEDBACK_2D112 CALL Agrif_Update_Variable(e1u_id,procname = updateU2d)113 CALL Agrif_Update_Variable(e2v_id,procname = updateV2d)114 # else115 CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d)116 CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)117 # endif118 !119 # if ! defined DECAL_FEEDBACK_2D120 ! Account for updated thicknesses at boundary edges121 IF (.NOT.ln_linssh) THEN122 ! CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy)123 ! CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy)124 ENDIF125 # endif126 !127 IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN128 ! Update time integrated transports129 # if ! defined DECAL_FEEDBACK_2D130 CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b)131 CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b)132 # else133 CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b)134 CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b)135 # endif136 END IF137 124 ! 138 125 use_sign_north = .FALSE. 126 l_vremap = .FALSE. 139 127 ! 140 128 END SUBROUTINE Agrif_Update_Dyn … … 150 138 Agrif_SpecialValueFineGrid = 0._wp 151 139 # if ! defined DECAL_FEEDBACK_2D 152 CALL Agrif_Update_Variable(sshn_id, procname = updateSSH)140 CALL Agrif_Update_Variable(sshn_id,locupdate=(/ nn_shift_bar,-2/), procname = updateSSH) 153 141 # else 154 CALL Agrif_Update_Variable(sshn_id,locupdate=(/1 ,0/),procname = updateSSH)142 CALL Agrif_Update_Variable(sshn_id,locupdate=(/1+nn_shift_bar,-2/),procname = updateSSH) 155 143 # endif 156 144 ! … … 158 146 ! 159 147 # if defined VOL_REFLUX 160 IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 148 IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 161 149 use_sign_north = .TRUE. 162 150 sign_north = -1._wp 163 151 ! Refluxing on ssh: 164 152 # if defined DECAL_FEEDBACK_2D 165 CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/ 0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu)166 CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1 , 1/),locupdate2=(/0, 0/),procname = reflux_sshv)153 CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/nn_shift_bar,nn_shift_bar/),locupdate2=(/1+nn_shift_bar,1+nn_shift_bar/),procname = reflux_sshu) 154 CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1+nn_shift_bar,1+nn_shift_bar/),locupdate2=(/nn_shift_bar,nn_shift_bar/),procname = reflux_sshv) 167 155 # else 168 CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1 ,-1/),locupdate2=(/ 0, 0/),procname = reflux_sshu)169 CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ 0, 0/),locupdate2=(/-1,-1/),procname = reflux_sshv)156 CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1+nn_shift_bar,-1+nn_shift_bar/),locupdate2=(/nn_shift_bar, nn_shift_bar/),procname = reflux_sshu) 157 CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ nn_shift_bar, nn_shift_bar/),locupdate2=(/-1+nn_shift_bar,-1+nn_shift_bar/),procname = reflux_sshv) 170 158 # endif 171 159 use_sign_north = .FALSE. … … 320 308 #endif 321 309 322 #if defined key_vertical323 310 324 311 SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) … … 340 327 ! 341 328 IF (before) THEN 342 !jc_alt 343 ! AGRIF_SpecialValue = -999._wp 344 DO jn = n1,n2-1 329 IF ( l_vremap ) THEN 330 DO jn = n1,n2-1 331 DO jk=k1,k2 332 DO jj=j1,j2 333 DO ji=i1,i2 334 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 335 END DO 336 END DO 337 END DO 338 END DO 345 339 DO jk=k1,k2 346 340 DO jj=j1,j2 347 341 DO ji=i1,i2 348 !jc_alt 349 ! tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) & 350 ! & * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp 351 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 342 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 352 343 END DO 353 344 END DO 354 345 END DO 355 END DO 356 DO jk=k1,k2 346 ELSE 347 DO jn = 1,jpts 348 DO jk=k1,k2 349 DO jj=j1,j2 350 DO ji=i1,i2 351 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) 352 END DO 353 END DO 354 END DO 355 END DO 356 357 ENDIF 358 ELSE 359 IF ( l_vremap ) THEN 360 tabres_child(:,:,:,:) = 0._wp 361 AGRIF_SpecialValue = 0._wp 357 362 DO jj=j1,j2 358 363 DO ji=i1,i2 359 !jc_alt 360 ! tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & 361 ! & + (tmask(ji,jj,jk) - 1._wp) * 999._wp 362 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 363 END DO 364 END DO 365 END DO 366 ELSE 367 tabres_child(:,:,:,:) = 0._wp 368 AGRIF_SpecialValue = 0._wp 369 DO jj=j1,j2 370 DO ji=i1,i2 371 N_in = 0 372 DO jk=k1,k2 !k2 = jpk of child grid 373 ! jc_alt 374 ! IF (tabres(ji,jj,jk,n2) < -900._wp ) EXIT 375 IF (tabres(ji,jj,jk,n2) == 0._wp ) EXIT 376 N_in = N_in + 1 377 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 378 h_in(N_in) = tabres(ji,jj,jk,n2) 364 N_in = 0 365 DO jk=k1,k2 !k2 = jpk of child grid 366 IF (tabres(ji,jj,jk,n2) <= 1.e-6_wp ) EXIT 367 N_in = N_in + 1 368 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 369 h_in(N_in) = tabres(ji,jj,jk,n2) 370 ENDDO 371 N_out = 0 372 DO jk=1,jpk ! jpk of parent grid 373 IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF 374 N_out = N_out + 1 375 h_out(N_out) = e3t(ji,jj,jk,Kmm_a) 376 ENDDO 377 IF (N_in*N_out > 0) THEN !Remove this? 378 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 379 ENDIF 379 380 ENDDO 380 N_out = 0381 DO jk=1,jpk ! jpk of parent grid382 IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF383 N_out = N_out + 1384 h_out(N_out) = e3t(ji,jj,jk,Kmm_a)385 ENDDO386 IF (N_in*N_out > 0) THEN !Remove this?387 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)388 ENDIF389 381 ENDDO 390 ENDDO 391 392 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 393 ! Add asselin part 382 383 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 384 ! Add asselin part 385 DO jn = 1,jpts 386 DO jk = 1, jpkm1 387 DO jj = j1, j2 388 DO ji = i1, i2 389 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 390 ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 391 ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a) 392 ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 393 ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) & 394 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 395 ENDIF 396 END DO 397 END DO 398 END DO 399 END DO 400 ENDIF 394 401 DO jn = 1,jpts 395 402 DO jk = 1, jpkm1 396 403 DO jj = j1, j2 397 404 DO ji = i1, i2 398 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 399 ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 400 ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a) 401 ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 402 ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) & 403 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 404 ENDIF 405 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 406 ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) 407 END IF 405 408 END DO 406 409 END DO 407 410 END DO 408 411 END DO 409 ENDIF 410 DO jn = 1,jpts 411 DO jk = 1, jpkm1 412 DO jj = j1, j2 413 DO ji = i1, i2 414 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 415 ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) 416 END IF 412 ELSE 413 DO jn = 1,jpts 414 tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 415 & * tmask(i1:i2,j1:j2,k1:k2) 416 ENDDO 417 418 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 419 ! Add asselin part 420 DO jn = 1,jpts 421 DO jk = k1, k2 422 DO jj = j1, j2 423 DO ji = i1, i2 424 IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 425 ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 426 ztnu = tabres(ji,jj,jk,jn) 427 ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 428 ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) & 429 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 430 ENDIF 431 END DO 432 END DO 417 433 END DO 418 434 END DO 419 END DO 420 END DO 421 ! 435 ENDIF 436 DO jn = 1,jpts 437 DO jk=k1,k2 438 DO jj=j1,j2 439 DO ji=i1,i2 440 IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 441 ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a) 442 END IF 443 END DO 444 END DO 445 END DO 446 END DO 447 ! 448 ENDIF 422 449 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 423 450 ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kmm_a) 424 ENDIF 451 ENDIF 425 452 ENDIF 426 453 ! 427 454 END SUBROUTINE updateTS 428 455 429 # else430 431 SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )432 !!---------------------------------------------433 !! *** ROUTINE updateT ***434 !!---------------------------------------------435 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2436 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres437 LOGICAL, INTENT(in) :: before438 !!439 INTEGER :: ji,jj,jk,jn440 REAL(wp) :: ztb, ztnu, ztno441 !!---------------------------------------------442 !443 IF (before) THEN444 DO jn = 1,jpts445 DO jk=k1,k2446 DO jj=j1,j2447 DO ji=i1,i2448 !> jc tmp449 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk)450 ! tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a)451 !< jc tmp452 END DO453 END DO454 END DO455 END DO456 ELSE457 !> jc tmp458 DO jn = 1,jpts459 tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) &460 & * tmask(i1:i2,j1:j2,k1:k2)461 ENDDO462 !< jc tmp463 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN464 ! Add asselin part465 DO jn = 1,jpts466 DO jk = k1, k2467 DO jj = j1, j2468 DO ji = i1, i2469 IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN470 ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used471 ztnu = tabres(ji,jj,jk,jn)472 ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a)473 ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) ) &474 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a)475 ENDIF476 END DO477 END DO478 END DO479 END DO480 ENDIF481 DO jn = 1,jpts482 DO jk=k1,k2483 DO jj=j1,j2484 DO ji=i1,i2485 IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN486 ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a)487 END IF488 END DO489 END DO490 END DO491 END DO492 !493 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN494 ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a)495 ENDIF496 !497 ENDIF498 !499 END SUBROUTINE updateTS500 501 # endif502 503 # if defined key_vertical504 456 505 457 SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) … … 513 465 INTEGER :: ji, jj, jk 514 466 REAL(wp):: zrhoy, zub, zunu, zuno 467 REAL(wp), DIMENSION(jpi,jpj) :: zpgu ! 2D workspace 515 468 ! VERTICAL REFINEMENT BEGIN 516 469 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 517 470 REAL(wp) :: h_in(k1:k2) 518 471 REAL(wp) :: h_out(1:jpk) 519 INTEGER :: N_in, N_out 520 REAL(wp) :: h_diff, excess, thick472 INTEGER :: N_in, N_out, N_in_save, N_out_save 473 REAL(wp) :: zhmin, zd 521 474 REAL(wp) :: tabin(k1:k2) 522 475 ! VERTICAL REFINEMENT END … … 525 478 IF( before ) THEN 526 479 zrhoy = Agrif_Rhoy() 527 !jc_alt528 ! AGRIF_SpecialValue = -999._wp529 480 DO jk=k1,k2 481 tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) & 482 & * umask(i1:i2,j1:j2,jk) * uu(i1:i2,j1:j2,jk,Kmm_a) 483 END DO 484 485 IF ( l_vremap ) THEN 486 DO jk=k1,k2 487 tabres(i1:i2,j1:j2,jk,2) = zrhoy * umask(i1:i2,j1:j2,jk) * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) 488 END DO 489 ENDIF 490 491 ELSE 492 493 tabres_child(:,:,:) = 0._wp 494 495 IF ( l_vremap ) THEN 496 530 497 DO jj=j1,j2 531 498 DO ji=i1,i2 532 !jc_alt 533 ! tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a) & 534 ! & + (umask(ji,jj,jk)-1._wp)*999._wp 535 tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a) 536 !jc_alt 537 ! tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) & 538 ! & + (umask(ji,jj,jk)-1._wp)*999._wp 539 tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 540 END DO 541 END DO 542 END DO 543 ELSE 544 tabres_child(:,:,:) = 0. 545 AGRIF_SpecialValue = 0._wp 546 DO jj=j1,j2 547 DO ji=i1,i2 548 N_in = 0 549 h_in(:) = 0._wp 550 tabin(:) = 0._wp 551 DO jk=k1,k2 !k2=jpk of child grid 552 !jc_alt 553 ! IF( tabres(ji,jj,jk,2) < -900._wp) EXIT 554 IF( tabres(ji,jj,jk,2) == 0.) EXIT 555 N_in = N_in + 1 556 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 557 h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 499 N_in = 0 500 h_in(:) = 0._wp 501 tabin(:) = 0._wp 502 DO jk=k1,k2 !k2=jpk of child grid 503 IF( tabres(ji,jj,jk,2)*r1_e2u(ji,jj) <= 1.e-6_wp ) EXIT 504 N_in = N_in + 1 505 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 506 h_in(N_in) = tabres(ji,jj,jk,2) * r1_e2u(ji,jj) 507 ENDDO 508 N_out = 0 509 DO jk=1,jpk 510 IF (umask(ji,jj,jk) == 0._wp) EXIT 511 N_out = N_out + 1 512 h_out(N_out) = e3u(ji,jj,jk,Kmm_a) 513 ENDDO 514 IF (N_in * N_out > 0) THEN 515 ! Deal with potentially different depths at velocity points: 516 N_in_save = N_in 517 N_out_save = N_out 518 IF ( ABS(sum(h_out(1:N_out))-sum(h_in(1:N_in))) > 1.e-6_wp ) THEN 519 zhmin = MIN(sum(h_out(1:N_out)), sum(h_in(1:N_in))) 520 zd = 0._wp 521 DO jk=1, N_in_save 522 IF ( (zd + h_in(jk)) > zhmin-1.e-6) THEN 523 N_in = jk 524 h_in(jk) = zhmin - zd 525 EXIT 526 ENDIF 527 zd = zd + h_in(jk) 528 END DO 529 zd = 0._wp 530 DO jk=1, N_out_save 531 IF ( (zd + h_out(jk)) > zhmin-1.e-6) THEN 532 N_out = jk 533 h_out(jk) = zhmin - zd 534 EXIT 535 ENDIF 536 zd = zd + h_out(jk) 537 END DO 538 END IF 539 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,1) 540 IF (N_out < N_out_save) tabres_child(ji,jj,N_out+1:N_out_save) = tabres_child(ji,jj,N_out) 541 ENDIF 558 542 ENDDO 559 N_out = 0560 DO jk=1,jpk561 IF (umask(ji,jj,jk) == 0) EXIT562 N_out = N_out + 1563 h_out(N_out) = e3u(ji,jj,jk,Kmm_a)564 ENDDO565 IF (N_in * N_out > 0) THEN566 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))567 excess = 0._wp568 IF (h_diff < -1.e-4) THEN569 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid.570 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.571 DO jk=N_in,1,-1572 thick = MIN(-1*h_diff, h_in(jk))573 excess = excess + tabin(jk)*thick*e2u(ji,jj)574 tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))575 h_diff = h_diff + thick576 IF ( h_diff == 0) THEN577 N_in = jk578 h_in(jk) = h_in(jk) - thick579 EXIT580 ENDIF581 ENDDO582 ENDIF583 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,1)584 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out))585 ENDIF586 543 ENDDO 587 ENDDO 544 545 ELSE 546 DO jk=1,jpk 547 DO jj=j1,j2 548 DO ji=i1,i2 549 tabres_child(ji,jj,jk) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm_a) 550 END DO 551 END DO 552 END DO 553 ENDIF 588 554 ! 589 555 DO jk=1,jpk … … 603 569 END DO 604 570 ! 571 ! Correct now and before transports: 572 DO jj=j1,j2 573 DO ji=i1,i2 574 zpgu(ji,jj) = 0._wp 575 DO jk=1,jpkm1 576 zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a) 577 END DO 578 ! 579 DO jk=1,jpkm1 580 uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + & 581 & (uu_b(ji,jj,Kmm_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kmm_a)) * umask(ji,jj,jk) 582 END DO 583 ! 584 zpgu(ji,jj) = 0._wp 585 DO jk=1,jpkm1 586 zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) 587 END DO 588 ! 589 DO jk=1,jpkm1 590 uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + & 591 & (uu_b(ji,jj,Kbb_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kbb_a)) * umask(ji,jj,jk) 592 END DO 593 ! 594 END DO 595 END DO 596 ! 605 597 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 606 598 uu(i1:i2,j1:j2,1:jpkm1,Kbb_a) = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a) … … 611 603 END SUBROUTINE updateu 612 604 613 #else 614 615 SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 616 !!--------------------------------------------- 617 !! *** ROUTINE updateu *** 618 !!--------------------------------------------- 619 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 620 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 621 LOGICAL , INTENT(in ) :: before 622 ! 623 INTEGER :: ji, jj, jk 624 REAL(wp) :: zrhoy, zub, zunu, zuno 625 !!--------------------------------------------- 626 ! 627 IF( before ) THEN 628 zrhoy = Agrif_Rhoy() 629 DO jk = k1, k2 630 tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) * uu(i1:i2,j1:j2,jk,Kmm_a) 631 END DO 632 ELSE 633 DO jk=k1,k2 634 DO jj=j1,j2 635 DO ji=i1,i2 636 tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) 637 ! 638 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part 639 zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 640 zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) 641 zunu = tabres(ji,jj,jk,1) 642 uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) & 643 & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) 644 ENDIF 645 ! 646 uu(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm_a) 647 END DO 648 END DO 649 END DO 650 ! 651 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 652 uu(i1:i2,j1:j2,k1:k2,Kbb_a) = uu(i1:i2,j1:j2,k1:k2,Kmm_a) 653 ENDIF 654 ! 655 ENDIF 656 ! 657 END SUBROUTINE updateu 658 659 # endif 660 661 SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 662 !!--------------------------------------------- 663 !! *** ROUTINE correct_u_bdy *** 605 606 SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 607 !!--------------------------------------------- 608 !! *** ROUTINE updatev *** 664 609 !!--------------------------------------------- 665 610 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 666 611 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 667 612 LOGICAL , INTENT(in ) :: before 668 INTEGER , INTENT(in) :: nb, ndir669 !!670 LOGICAL :: western_side, eastern_side671 !672 INTEGER :: jj, jk673 REAL(wp) :: zcor674 !!---------------------------------------------675 !676 IF( .NOT.before ) THEN677 !678 western_side = (nb == 1).AND.(ndir == 1)679 eastern_side = (nb == 1).AND.(ndir == 2)680 !681 IF (western_side) THEN682 DO jj=j1,j2683 zcor = uu_b(i1-1,jj,Kmm_a) * hu(i1-1,jj,Krhs_a) * r1_hu(i1-1,jj,Kmm_a) - uu_b(i1-1,jj,Kmm_a)684 uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor685 DO jk=1,jpkm1686 uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk)687 END DO688 END DO689 ENDIF690 !691 IF (eastern_side) THEN692 DO jj=j1,j2693 zcor = uu_b(i2+1,jj,Kmm_a) * hu(i2+1,jj,Krhs_a) * r1_hu(i2+1,jj,Kmm_a) - uu_b(i2+1,jj,Kmm_a)694 uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor695 DO jk=1,jpkm1696 uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk)697 END DO698 END DO699 ENDIF700 !701 ENDIF702 !703 END SUBROUTINE correct_u_bdy704 705 # if defined key_vertical706 707 SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )708 !!---------------------------------------------709 !! *** ROUTINE updatev ***710 !!---------------------------------------------711 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2712 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres713 LOGICAL , INTENT(in ) :: before714 613 ! 715 614 INTEGER :: ji, jj, jk 716 615 REAL(wp) :: zrhox, zvb, zvnu, zvno 616 REAL(wp), DIMENSION(jpi,jpj) :: zpgv ! 2D workspace 717 617 ! VERTICAL REFINEMENT BEGIN 718 618 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 719 619 REAL(wp) :: h_in(k1:k2) 720 620 REAL(wp) :: h_out(1:jpk) 721 INTEGER :: N_in, N_out722 REAL(wp) :: h_diff, excess, thick621 INTEGER :: N_in, N_out, N_in_save, N_out_save 622 REAL(wp) :: zhmin, zd 723 623 REAL(wp) :: tabin(k1:k2) 724 624 ! VERTICAL REFINEMENT END … … 727 627 IF( before ) THEN 728 628 zrhox = Agrif_Rhox() 729 !jc_alt730 ! AGRIF_SpecialValue = -999._wp731 629 DO jk=k1,k2 630 tabres(i1:i2,j1:j2,jk,1) = zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Kmm_a) & 631 & * vmask(i1:i2,j1:j2,jk) * vv(i1:i2,j1:j2,jk,Kmm_a) 632 END DO 633 634 IF ( l_vremap ) THEN 635 DO jk=k1,k2 636 tabres(i1:i2,j1:j2,jk,2) = zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk) 637 END DO 638 ENDIF 639 640 ELSE 641 642 tabres_child(:,:,:) = 0._wp 643 644 IF ( l_vremap ) THEN 645 732 646 DO jj=j1,j2 733 647 DO ji=i1,i2 734 !jc_alt 735 ! tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) & 736 ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp 737 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) 738 !jc_alt 739 ! tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) & 740 ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp 741 tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) 742 END DO 743 END DO 744 END DO 745 ELSE 746 tabres_child(:,:,:) = 0. 747 AGRIF_SpecialValue = 0._wp 748 DO jj=j1,j2 749 DO ji=i1,i2 750 N_in = 0 751 DO jk=k1,k2 752 !jc_alt 753 ! IF (tabres(ji,jj,jk,2) < -900._wp) EXIT 754 IF (tabres(ji,jj,jk,2) == 0) EXIT 755 N_in = N_in + 1 756 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 757 h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) 648 N_in = 0 649 DO jk=k1,k2 650 IF (tabres(ji,jj,jk,2)* r1_e1v(ji,jj) <= 1.e-6_wp) EXIT 651 N_in = N_in + 1 652 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 653 h_in(N_in) = tabres(ji,jj,jk,2) * r1_e1v(ji,jj) 654 ENDDO 655 N_out = 0 656 DO jk=1,jpk 657 IF (vmask(ji,jj,jk) == 0) EXIT 658 N_out = N_out + 1 659 h_out(N_out) = e3v(ji,jj,jk,Kmm_a) 660 ENDDO 661 IF (N_in * N_out > 0) THEN 662 ! Deal with potentially different depths at velocity points: 663 N_in_save = N_in 664 N_out_save = N_out 665 IF ( ABS(sum(h_out(1:N_out))-sum(h_in(1:N_in))) > 1.e-6_wp ) THEN 666 zhmin = MIN(sum(h_out(1:N_out)), sum(h_in(1:N_in))) 667 zd = 0._wp 668 DO jk=1, N_in_save 669 IF ( (zd + h_in(jk)) > zhmin-1.e-6) THEN 670 N_in = jk 671 h_in(jk) = zhmin - zd 672 EXIT 673 ENDIF 674 zd = zd + h_in(jk) 675 END DO 676 zd = 0._wp 677 DO jk=1, N_out_save 678 IF ( (zd + h_out(jk)) > zhmin-1.e-6) THEN 679 N_out = jk 680 h_out(jk) = zhmin - zd 681 EXIT 682 ENDIF 683 zd = zd + h_out(jk) 684 END DO 685 END IF 686 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,1) 687 IF (N_out < N_out_save) tabres_child(ji,jj,N_out+1:N_out_save) = tabres_child(ji,jj,N_out) 688 ENDIF 758 689 ENDDO 759 N_out = 0760 DO jk=1,jpk761 IF (vmask(ji,jj,jk) == 0) EXIT762 N_out = N_out + 1763 h_out(N_out) = e3v(ji,jj,jk,Kmm_a)764 ENDDO765 IF (N_in * N_out > 0) THEN766 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))767 excess = 0._wp768 IF (h_diff < -1.e-4) then769 !Even if bathy at T points match it's possible for the V points to be deeper in the child grid.770 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell.771 DO jk=N_in,1,-1772 thick = MIN(-1*h_diff, h_in(jk))773 excess = excess + tabin(jk)*thick*e2u(ji,jj)774 tabin(jk) = tabin(jk)*(1. - thick/h_in(jk))775 h_diff = h_diff + thick776 IF ( h_diff == 0) THEN777 N_in = jk778 h_in(jk) = h_in(jk) - thick779 EXIT780 ENDIF781 ENDDO782 ENDIF783 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,1)784 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out))785 ENDIF786 690 ENDDO 787 ENDDO 691 692 ELSE 693 DO jk=1,jpk 694 DO jj=j1,j2 695 DO ji=i1,i2 696 tabres_child(ji,jj,jk) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm_a) 697 END DO 698 END DO 699 END DO 700 ENDIF 788 701 ! 789 702 DO jk=1,jpkm1 … … 803 716 END DO 804 717 ! 718 ! Correct now and before transports: 719 DO jj=j1,j2 720 DO ji=i1,i2 721 zpgv(ji,jj) = 0._wp 722 DO jk=1,jpkm1 723 zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 724 END DO 725 ! 726 DO jk=1,jpkm1 727 vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + & 728 & (vv_b(ji,jj,Kmm_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kmm_a)) * vmask(ji,jj,jk) 729 END DO 730 ! 731 zpgv(ji,jj) = 0._wp 732 DO jk=1,jpkm1 733 zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) 734 END DO 735 ! 736 DO jk=1,jpkm1 737 vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + & 738 & (vv_b(ji,jj,Kbb_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kbb_a)) * vmask(ji,jj,jk) 739 END DO 740 ! 741 END DO 742 END DO 743 ! 805 744 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 806 745 vv(i1:i2,j1:j2,1:jpkm1,Kbb_a) = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a) … … 810 749 ! 811 750 END SUBROUTINE updatev 812 813 # else814 815 SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)816 !!---------------------------------------------817 !! *** ROUTINE updatev ***818 !!---------------------------------------------819 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2820 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres821 LOGICAL , INTENT(in ) :: before822 !823 INTEGER :: ji, jj, jk824 REAL(wp) :: zrhox, zvb, zvnu, zvno825 !!---------------------------------------------826 !827 IF (before) THEN828 zrhox = Agrif_Rhox()829 DO jk=k1,k2830 DO jj=j1,j2831 DO ji=i1,i2832 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)833 END DO834 END DO835 END DO836 ELSE837 DO jk=k1,k2838 DO jj=j1,j2839 DO ji=i1,i2840 tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj)841 !842 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part843 zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used844 zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a)845 zvnu = tabres(ji,jj,jk,1)846 vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) &847 & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)848 ENDIF849 !850 vv(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm_a)851 END DO852 END DO853 END DO854 !855 IF ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN856 vv(i1:i2,j1:j2,k1:k2,Kbb_a) = vv(i1:i2,j1:j2,k1:k2,Kmm_a)857 ENDIF858 !859 ENDIF860 !861 END SUBROUTINE updatev862 863 # endif864 865 SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )866 !!---------------------------------------------867 !! *** ROUTINE correct_v_bdy ***868 !!---------------------------------------------869 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2870 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres871 LOGICAL , INTENT(in ) :: before872 INTEGER , INTENT(in) :: nb, ndir873 !!874 LOGICAL :: southern_side, northern_side875 !876 INTEGER :: ji, jk877 REAL(wp) :: zcor878 !!---------------------------------------------879 !880 IF( .NOT.before ) THEN881 !882 southern_side = (nb == 2).AND.(ndir == 1)883 northern_side = (nb == 2).AND.(ndir == 2)884 !885 IF (southern_side) THEN886 DO ji=i1,i2887 zcor = vv_b(ji,j1-1,Kmm_a) * hv(ji,j1-1,Krhs_a) * r1_hv(ji,j1-1,Kmm_a) - vv_b(ji,j1-1,Kmm_a)888 vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor889 DO jk=1,jpkm1890 vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk)891 END DO892 END DO893 ENDIF894 !895 IF (northern_side) THEN896 DO ji=i1,i2897 zcor = vv_b(ji,j2+1,Kmm_a) * hv(ji,j2+1,Krhs_a) * r1_hv(ji,j2+1,Kmm_a) - vv_b(ji,j2+1,Kmm_a)898 vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor899 DO jk=1,jpkm1900 vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk)901 END DO902 END DO903 ENDIF904 !905 ENDIF906 !907 END SUBROUTINE correct_v_bdy908 751 909 752 … … 935 778 tabres(ji,jj) = tabres(ji,jj) * r1_e2u(ji,jj) 936 779 ! 937 ! Update "now" 3d velocities:938 zpgu(ji,jj) = 0._wp939 DO jk=1,jpkm1940 zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a)941 END DO942 !943 zcorr = (tabres(ji,jj) - zpgu(ji,jj)) * r1_hu(ji,jj,Kmm_a)944 DO jk=1,jpkm1945 uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)946 END DO947 !948 780 ! Update barotropic velocities: 949 781 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN … … 955 787 uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1) 956 788 ! 957 ! Correct "before" velocities to hold correct bt component:958 zpgu(ji,jj) = 0.e0959 DO jk=1,jpkm1960 zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a)961 END DO962 !963 zcorr = uu_b(ji,jj,Kbb_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kbb_a)964 DO jk=1,jpkm1965 uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk)966 END DO967 !968 789 END DO 969 790 END DO … … 1002 823 DO ji=i1,i2 1003 824 tabres(ji,jj) = tabres(ji,jj) * r1_e1v(ji,jj) 1004 !1005 ! Update "now" 3d velocities:1006 zpgv(ji,jj) = 0.e01007 DO jk=1,jpkm11008 zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a)1009 END DO1010 !1011 zcorr = (tabres(ji,jj) - zpgv(ji,jj)) * r1_hv(ji,jj,Kmm_a)1012 DO jk=1,jpkm11013 vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)1014 END DO1015 !1016 825 ! Update barotropic velocities: 1017 826 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN … … 1023 832 vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1) 1024 833 ! 1025 ! Correct "before" velocities to hold correct bt component:1026 zpgv(ji,jj) = 0.e01027 DO jk=1,jpkm11028 zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a)1029 END DO1030 !1031 zcorr = vv_b(ji,jj,Kbb_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kbb_a)1032 DO jk=1,jpkm11033 vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk)1034 END DO1035 !1036 834 END DO 1037 835 END DO … … 1120 918 ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor 1121 919 ! Update corrective fluxes: 1122 un_bf(ji,jj) = un_bf(ji,jj) + zcor920 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) un_bf(ji,jj) = un_bf(ji,jj) + zcor 1123 921 ! Update half step back fluxes: 1124 922 ub2_b(ji,jj) = tabres(ji,jj) … … 1208 1006 vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor 1209 1007 ! Update corrective fluxes: 1210 vn_bf(ji,jj) = vn_bf(ji,jj) + zcor1008 IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) vn_bf(ji,jj) = vn_bf(ji,jj) + zcor 1211 1009 ! Update half step back fluxes: 1212 1010 vb2_b(ji,jj) = tabres(ji,jj) … … 1479 1277 #endif 1480 1278 1279 SUBROUTINE Agrif_Check_parent_bat( ) 1280 !!---------------------------------------------------------------------- 1281 !! *** ROUTINE Agrif_Check_parent_bat *** 1282 !!---------------------------------------------------------------------- 1283 ! 1284 IF (( .NOT.ln_agrif_2way ).OR.(.NOT.ln_chk_bathy).OR.(Agrif_Root())) RETURN 1285 ! 1286 Agrif_UseSpecialValueInUpdate = .FALSE. 1287 ! 1288 IF(lwp) WRITE(numout,*) ' ' 1289 IF(lwp) WRITE(numout,*) 'AGRIF: Check parent volume at Level:', Agrif_Level() 1290 ! 1291 # if ! defined DECAL_FEEDBACK 1292 CALL Agrif_Update_Variable(batupd_id,procname = update_bat) 1293 # else 1294 CALL Agrif_Update_Variable(batupd_id,locupdate=(/1,0/),procname = update_bat) 1295 # endif 1296 ! 1297 kindic_agr = Agrif_Parent(kindic_agr) 1298 CALL mpp_sum( 'Agrif_Check_parent_bat', kindic_agr ) 1299 1300 IF( kindic_agr /= 0 ) THEN 1301 CALL ctl_stop('==> Averaged Bathymetry does not match parent volume') 1302 ELSE 1303 IF(lwp) WRITE(numout,*) '==> Averaged Bathymetry matches parent ' 1304 IF(lwp) WRITE(numout,*) '' 1305 ENDIF 1306 ! 1307 END SUBROUTINE Agrif_Check_parent_bat 1308 1309 SUBROUTINE update_bat(ptab, i1, i2, j1, j2, before ) 1310 !!--------------------------------------------- 1311 !! *** ROUTINE update_bat *** 1312 !!--------------------------------------------- 1313 REAL(wp), DIMENSION(i1:i2,j1:j2) :: ptab 1314 INTEGER, INTENT(in) :: i1, i2, j1, j2 1315 LOGICAL, INTENT(in) :: before 1316 INTEGER :: ji, jj 1317 ! 1318 !!--------------------------------------------- 1319 ! 1320 IF( before ) THEN 1321 ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 1322 ELSE 1323 kindic_agr = 0 1324 ! 1325 DO jj=j1,j2 1326 DO ji=i1,i2 1327 IF ( (ssmask(ji,jj).NE.0._wp).AND.& 1328 & (ABS(ptab(ji,jj)-ht_0(ji,jj)).GE.1.e-6) ) THEN 1329 kindic_agr = kindic_agr + 1 1330 ENDIF 1331 END DO 1332 END DO 1333 ! 1334 ENDIF 1335 ! 1336 END SUBROUTINE update_bat 1337 1481 1338 #else 1482 1339 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.