- Timestamp:
- 2019-12-12T17:41:04+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_oce_update.F90
r11949 r12229 1 # define TWO_WAY /* TWO WAY NESTING*/2 #undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/3 #undef VOL_REFLUX /* VOLUME REFLUXING*/1 #undef DECAL_FEEDBACK /* SEPARATION of INTERFACES */ 2 #undef DECAL_FEEDBACK_2D /* SEPARATION of INTERFACES (Barotropic mode) */ 3 #undef VOL_REFLUX /* VOLUME REFLUXING*/ 4 4 5 5 MODULE agrif_oce_update … … 25 25 USE lib_mpp ! MPP library 26 26 USE domvvl ! Need interpolation routines 27 USE vremap ! Vertical remapping 27 28 28 29 IMPLICIT NONE … … 46 47 IF (Agrif_Root()) RETURN 47 48 ! 48 #if defined TWO_WAY49 49 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers from grid Number',Agrif_Fixed() 50 50 51 #if defined key_vertical 52 ! Effect of this has to be carrefully checked 53 ! depending on what the nesting tools ensure for 54 ! volume conservation: 55 Agrif_UseSpecialValueInUpdate = .FALSE. 56 #else 51 57 Agrif_UseSpecialValueInUpdate = .TRUE. 58 #endif 52 59 Agrif_SpecialValueFineGrid = 0._wp 53 60 ! … … 64 71 Agrif_UseSpecialValueInUpdate = .FALSE. 65 72 ! 66 #endif67 73 ! 68 74 END SUBROUTINE Agrif_Update_Tra … … 75 81 IF (Agrif_Root()) RETURN 76 82 ! 77 #if defined TWO_WAY78 83 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed() 79 84 … … 95 100 # endif 96 101 97 # if ! defined DECAL_FEEDBACK 102 # if ! defined DECAL_FEEDBACK_2D 98 103 CALL Agrif_Update_Variable(e1u_id,procname = updateU2d) 99 104 CALL Agrif_Update_Variable(e2v_id,procname = updateV2d) … … 103 108 # endif 104 109 ! 105 # if ! defined DECAL_FEEDBACK 110 # if ! defined DECAL_FEEDBACK_2D 106 111 ! Account for updated thicknesses at boundary edges 107 112 IF (.NOT.ln_linssh) THEN … … 113 118 IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 114 119 ! Update time integrated transports 115 # if ! defined DECAL_FEEDBACK 120 # if ! defined DECAL_FEEDBACK_2D 116 121 CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b) 117 122 CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b) … … 121 126 # endif 122 127 END IF 123 #endif124 128 ! 125 129 END SUBROUTINE Agrif_Update_Dyn … … 131 135 ! 132 136 IF (Agrif_Root()) RETURN 133 !134 #if defined TWO_WAY135 137 ! 136 138 Agrif_UseSpecialValueInUpdate = .TRUE. 137 139 Agrif_SpecialValueFineGrid = 0. 138 # if ! defined DECAL_FEEDBACK 140 # if ! defined DECAL_FEEDBACK_2D 139 141 CALL Agrif_Update_Variable(sshn_id,procname = updateSSH) 140 142 # else … … 147 149 IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 148 150 ! Refluxing on ssh: 149 # if defined DECAL_FEEDBACK 151 # if defined DECAL_FEEDBACK_2D 150 152 CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu) 151 153 CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv) … … 157 159 # endif 158 160 ! 159 #endif160 !161 161 END SUBROUTINE Agrif_Update_ssh 162 162 … … 170 170 IF (Agrif_Root()) RETURN 171 171 ! 172 # if defined TWO_WAY173 174 172 Agrif_UseSpecialValueInUpdate = .TRUE. 175 173 Agrif_SpecialValueFineGrid = 0. … … 180 178 181 179 Agrif_UseSpecialValueInUpdate = .FALSE. 182 183 # endif184 180 185 181 END SUBROUTINE Agrif_Update_Tke … … 192 188 ! 193 189 IF (Agrif_Root()) RETURN 194 !195 #if defined TWO_WAY196 190 ! 197 191 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step() … … 209 203 CALL dom_vvl_update_UVF 210 204 CALL Agrif_ParentGrid_To_ChildGrid() 211 !212 #endif213 205 ! 214 206 END SUBROUTINE Agrif_Update_vvl … … 300 292 !! 301 293 INTEGER :: ji,jj,jk,jn 302 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 294 INTEGER :: N_in, N_out 295 REAL(wp) :: ztb, ztnu, ztno 303 296 REAL(wp) :: h_in(k1:k2) 304 297 REAL(wp) :: h_out(1:jpk) 305 INTEGER :: N_in, N_out 306 REAL(wp) :: zrho_xy, h_diff 307 REAL(wp) :: tabin(k1:k2,n1:n2) 298 REAL(wp) :: tabin(k1:k2,1:jpts) 299 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,1:jpts) :: tabres_child 308 300 !!--------------------------------------------- 309 301 ! 310 302 IF (before) THEN 311 AGRIF_SpecialValue = -999._wp 312 zrho_xy = Agrif_rhox() * Agrif_rhoy() 303 !jc_alt 304 ! AGRIF_SpecialValue = -999._wp 313 305 DO jn = n1,n2-1 314 306 DO jk=k1,k2 315 307 DO jj=j1,j2 316 308 DO ji=i1,i2 317 tabres(ji,jj,jk,jn) = ( ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) & 318 & * tmask(ji,jj,jk) + (tmask(ji,jj,jk) - 1._wp)*999._wp 309 !jc_alt 310 ! tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) & 311 ! & * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp 312 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 319 313 END DO 320 314 END DO … … 324 318 DO jj=j1,j2 325 319 DO ji=i1,i2 326 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & 327 & + (tmask(ji,jj,jk) - 1._wp)*999._wp 320 !jc_alt 321 ! tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & 322 ! & + (tmask(ji,jj,jk) - 1._wp) * 999._wp 323 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 328 324 END DO 329 325 END DO … … 336 332 N_in = 0 337 333 DO jk=k1,k2 !k2 = jpk of child grid 338 IF (tabres(ji,jj,jk,n2) == 0 ) EXIT 334 ! jc_alt 335 ! IF (tabres(ji,jj,jk,n2) < -900._wp ) EXIT 336 IF (tabres(ji,jj,jk,n2) == 0._wp ) EXIT 339 337 N_in = N_in + 1 340 338 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) … … 343 341 N_out = 0 344 342 DO jk=1,jpk ! jpk of parent grid 345 IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF343 IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF 346 344 N_out = N_out + 1 347 345 h_out(N_out) = e3t(ji,jj,jk,Kmm_a) 348 346 ENDDO 349 IF (N_in > 0) THEN !Remove this? 350 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 351 IF (h_diff < -1.e-4) THEN 352 print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 353 print *,h_in(1:N_in) 354 print *,h_out(1:N_out) 355 STOP 356 ENDIF 357 DO jn=n1,n2-1 358 CALL reconstructandremap( tabin(1:N_in,jn), h_in(1:N_in), tabres_child(ji,jj,1:N_out,jn), & 359 & h_out(1:N_out) , N_in , N_out ) 360 ENDDO 347 IF (N_in*N_out > 0) THEN !Remove this? 348 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) 361 349 ENDIF 362 350 ENDDO … … 365 353 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 366 354 ! Add asselin part 367 DO jn = n1,n2-1 368 DO jk=1,jpk 369 DO jj=j1,j2 370 DO ji=i1,i2 371 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 372 ts(ji,jj,jk,jn,Kbb_a) = ts(ji,jj,jk,jn,Kbb_a) & 373 & + atfp * ( tabres_child(ji,jj,jk,jn) & 374 & - ts(ji,jj,jk,jn,Kmm_a) & 375 & ) * tmask(ji,jj,jk) 355 DO jn = 1,jpts 356 DO jk = 1, jpkm1 357 DO jj = j1, j2 358 DO ji = i1, i2 359 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 360 ztb = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 361 ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a) 362 ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 363 ts(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) ) & 364 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 376 365 ENDIF 377 END DO378 END DO379 END DO380 END DO381 ENDIF 382 DO jn = n1,n2-1383 DO jk =1,jpk384 DO jj =j1,j2385 DO ji =i1,i2386 IF( tabres_child(ji,jj,jk,jn) .NE. 0.) THEN387 ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk)366 END DO 367 END DO 368 END DO 369 END DO 370 ENDIF 371 DO jn = 1,jpts 372 DO jk = 1, jpkm1 373 DO jj = j1, j2 374 DO ji = i1, i2 375 IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 376 ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) 388 377 END IF 389 378 END DO … … 391 380 END DO 392 381 END DO 382 ! 383 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 384 ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kbb_a) = ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kmm_a) 385 ENDIF 393 386 ENDIF 394 387 ! … … 416 409 !> jc tmp 417 410 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) 418 ! tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a)411 ! tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 419 412 !< jc tmp 420 413 END DO … … 440 433 ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 441 434 ts(ji,jj,jk,jn,Kbb_a) = ( ztb + atfp * ( ztnu - ztno) ) & 442 & 435 & * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 443 436 ENDIF 444 437 END DO … … 480 473 ! 481 474 INTEGER :: ji, jj, jk 482 REAL(wp):: zrhoy 475 REAL(wp):: zrhoy, zub, zunu, zuno 483 476 ! VERTICAL REFINEMENT BEGIN 484 477 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child … … 493 486 IF( before ) THEN 494 487 zrhoy = Agrif_Rhoy() 495 AGRIF_SpecialValue = -999._wp 488 !jc_alt 489 ! AGRIF_SpecialValue = -999._wp 496 490 DO jk=k1,k2 497 491 DO jj=j1,j2 498 492 DO ji=i1,i2 499 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) & 500 & + (umask(ji,jj,jk)-1._wp)*999._wp 501 tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) & 502 & + (umask(ji,jj,jk)-1._wp)*999._wp 493 !jc_alt 494 ! 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) & 495 ! & + (umask(ji,jj,jk)-1._wp)*999._wp 496 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) 497 !jc_alt 498 ! tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) & 499 ! & + (umask(ji,jj,jk)-1._wp)*999._wp 500 tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 503 501 END DO 504 502 END DO … … 513 511 tabin(:) = 0._wp 514 512 DO jk=k1,k2 !k2=jpk of child grid 515 IF( tabres(ji,jj,jk,2) < -900) EXIT 513 !jc_alt 514 ! IF( tabres(ji,jj,jk,2) < -900._wp) EXIT 515 IF( tabres(ji,jj,jk,2) == 0.) EXIT 516 516 N_in = N_in + 1 517 tabin(jk) 517 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 518 518 h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 519 519 ENDDO … … 525 525 ENDDO 526 526 IF (N_in * N_out > 0) THEN 527 h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) ) 527 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 528 excess = 0._wp 528 529 IF (h_diff < -1.e-4) THEN 529 530 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 530 531 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 531 excess = 0._wp532 532 DO jk=N_in,1,-1 533 533 thick = MIN(-1*h_diff, h_in(jk)) … … 542 542 ENDDO 543 543 ENDIF 544 CALL reconstructandremap( tabin(1:N_in) , h_in(1:N_in), tabres_child(ji,jj,1:N_out), & 545 & h_out(1:N_out), N_in , N_out ) 546 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/( e2u(ji,jj)*h_out(N_out) ) 544 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) 545 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 547 546 ENDIF 548 547 ENDDO 549 548 ENDDO 550 549 ! 551 550 DO jk=1,jpk 552 551 DO jj=j1,j2 553 552 DO ji=i1,i2 554 553 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 555 uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) & 556 & + atfp * ( tabres_child(ji,jj,jk) - uu(ji,jj,jk,Kmm_a) ) * umask(ji,jj,jk) 554 zub = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 555 zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) 556 zunu = tabres_child(ji,jj,jk) * e3u(ji,jj,jk,Kmm_a) 557 uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno) ) & 558 & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) 557 559 ENDIF 558 560 ! … … 561 563 END DO 562 564 END DO 565 ! 566 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 567 uu(i1:i2,j1:j2,1:jpkm1,Kbb_a) = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a) 568 ENDIF 569 ! 563 570 ENDIF 564 571 ! … … 582 589 zrhoy = Agrif_Rhoy() 583 590 DO jk = k1, k2 584 tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) & 585 & * uu(i1:i2,j1:j2,jk,Kmm_a) 591 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) 586 592 END DO 587 593 ELSE … … 595 601 zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) 596 602 zunu = tabres(ji,jj,jk,1) 597 uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno ) )&598 &* umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a)603 uu(ji,jj,jk,Kbb_a) = ( zub + atfp * ( zunu - zuno) ) & 604 & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) 599 605 ENDIF 600 606 ! … … 669 675 ! 670 676 INTEGER :: ji, jj, jk 671 REAL(wp) :: zrhox 677 REAL(wp) :: zrhox, zvb, zvnu, zvno 672 678 ! VERTICAL REFINEMENT BEGIN 673 679 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child … … 682 688 IF( before ) THEN 683 689 zrhox = Agrif_Rhox() 684 AGRIF_SpecialValue = -999._wp 690 !jc_alt 691 ! AGRIF_SpecialValue = -999._wp 685 692 DO jk=k1,k2 686 693 DO jj=j1,j2 687 694 DO ji=i1,i2 688 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) & 689 & * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) & 690 & + (vmask(ji,jj,jk)-1)*999._wp 691 tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) & 692 + (vmask(ji,jj,jk)-1)*999._wp 695 !jc_alt 696 ! 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) & 697 ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp 698 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) 699 !jc_alt 700 ! tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) & 701 ! & + (vmask(ji,jj,jk)-1._wp) * 999._wp 702 tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) 693 703 END DO 694 704 END DO … … 701 711 N_in = 0 702 712 DO jk=k1,k2 703 IF (tabres(ji,jj,jk,2) < -900) EXIT 713 !jc_alt 714 ! IF (tabres(ji,jj,jk,2) < -900._wp) EXIT 715 IF (tabres(ji,jj,jk,2) == 0) EXIT 704 716 N_in = N_in + 1 705 717 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) … … 713 725 ENDDO 714 726 IF (N_in * N_out > 0) THEN 715 h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) ) 727 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 728 excess = 0._wp 716 729 IF (h_diff < -1.e-4) then 717 !Even if bathy at T points match it's possible for the Upoints to be deeper in the child grid.730 !Even if bathy at T points match it's possible for the V points to be deeper in the child grid. 718 731 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 719 excess = 0._wp720 732 DO jk=N_in,1,-1 721 thick = MIN( -1._wp*h_diff, h_in(jk))733 thick = MIN(-1*h_diff, h_in(jk)) 722 734 excess = excess + tabin(jk)*thick*e2u(ji,jj) 723 tabin(jk) = tabin(jk)*(1. _wp- thick/h_in(jk))735 tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 724 736 h_diff = h_diff + thick 725 737 IF ( h_diff == 0) THEN … … 730 742 ENDDO 731 743 ENDIF 732 CALL reconstructandremap( tabin(1:N_in) , h_in(1:N_in), tabres_child(ji,jj,1:N_out), & 733 & h_out(1:N_out), N_in , N_out ) 744 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) 734 745 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 735 746 ENDIF 736 747 ENDDO 737 748 ENDDO 738 739 DO jk=1,jpk 749 ! 750 DO jk=1,jpkm1 740 751 DO jj=j1,j2 741 752 DO ji=i1,i2 742 ! 743 IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part 744 vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) & 745 & + atfp * ( tabres_child(ji,jj,jk) - vv(ji,jj,jk,Kmm_a) ) & 746 & * vmask(ji,jj,jk) 753 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 754 zvb = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 755 zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) 756 zvnu = tabres_child(ji,jj,jk) * e3v(ji,jj,jk,Kmm_a) 757 vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) ) & 758 & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) 747 759 ENDIF 748 760 ! … … 751 763 END DO 752 764 END DO 765 ! 766 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 767 vv(i1:i2,j1:j2,1:jpkm1,Kbb_a) = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a) 768 ENDIF 769 ! 753 770 ENDIF 754 771 ! … … 788 805 zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) 789 806 zvnu = tabres(ji,jj,jk,1) 790 vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) ) 791 &* vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a)807 vv(ji,jj,jk,Kbb_a) = ( zvb + atfp * ( zvnu - zvno) ) & 808 & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) 792 809 ENDIF 793 810 ! … … 890 907 ! Update barotropic velocities: 891 908 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 892 IF ( .NOT.( lk_agrif_fstep .AND. (neuler==0) ) ) THEN! Add asselin part909 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 893 910 zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a) 894 911 uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + atfp * zcorr * umask(ji,jj,1) … … 955 972 ! 956 973 ! Update barotropic velocities: 957 IF ( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. ( .NOT.ln_bt_fw )) ) THEN958 IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) ) THEN! Add asselin part974 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 975 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 959 976 zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a) 960 vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + atfp * zcorr *vmask(ji,jj,1)977 vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + atfp * zcorr * vmask(ji,jj,1) 961 978 END IF 962 979 ENDIF … … 1007 1024 DO jj=j1,j2 1008 1025 DO ji=i1,i2 1009 ssh(ji,jj,Kbb_a) = ssh(ji,jj,Kbb_a) & 1010 & + atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) & 1011 & * tmask(ji,jj,1) 1026 ssh(ji,jj,Kbb_a) = ssh(ji,jj,Kbb_a) & 1027 & + atfp * ( tabres(ji,jj) - ssh(ji,jj,Kmm_a) ) * tmask(ji,jj,1) 1012 1028 END DO 1013 1029 END DO … … 1103 1119 zcor = rdt * r1_e1e2t(i1 ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj)) 1104 1120 ssh(i1 ,jj,Kmm_a) = ssh(i1 ,jj,Kmm_a) + zcor 1105 IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) ) & 1106 & ssh(i1 ,jj,Kbb_a) = ssh(i1 ,jj,Kbb_a) + atfp * zcor 1121 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i1 ,jj,Kbb_a) = ssh(i1 ,jj,Kbb_a) + atfp * zcor 1107 1122 END DO 1108 1123 ENDIF … … 1111 1126 zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj)) 1112 1127 ssh(i2+1,jj,Kmm_a) = ssh(i2+1,jj,Kmm_a) + zcor 1113 IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) ) & 1114 & ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + atfp * zcor 1128 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(i2+1,jj,Kbb_a) = ssh(i2+1,jj,Kbb_a) + atfp * zcor 1115 1129 END DO 1116 1130 ENDIF … … 1191 1205 IF (southern_side) THEN 1192 1206 DO ji=i1,i2 1193 zcor = rdt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * ( vb2_b(ji,j1)-tabres(ji,j1))1207 zcor = rdt * r1_e1e2t(ji,j1 ) * e1v(ji,j1 ) * (vb2_b(ji,j1)-tabres(ji,j1)) 1194 1208 ssh(ji,j1 ,Kmm_a) = ssh(ji,j1 ,Kmm_a) + zcor 1195 IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) ) & 1196 & ssh(ji,j1 ,Kbb_a) = ssh(ji,j1,Kbb_a) + atfp * zcor 1209 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j1 ,Kbb_a) = ssh(ji,j1,Kbb_a) + atfp * zcor 1197 1210 END DO 1198 1211 ENDIF 1199 1212 IF (northern_side) THEN 1200 1213 DO ji=i1,i2 1201 zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * ( vb2_b(ji,j2)-tabres(ji,j2))1214 zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2 ) * (vb2_b(ji,j2)-tabres(ji,j2)) 1202 1215 ssh(ji,j2+1,Kmm_a) = ssh(ji,j2+1,Kmm_a) + zcor 1203 IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) ) & 1204 & ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + atfp * zcor 1216 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) ssh(ji,j2+1,Kbb_a) = ssh(ji,j2+1,Kbb_a) + atfp * zcor 1205 1217 END DO 1206 1218 ENDIF … … 1235 1247 END DO 1236 1248 END DO 1237 tabres(:,:,:,1) =tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy()1238 tabres(:,:,:,2) =tabres(:,:,:,2)*Agrif_Rhox()1239 tabres(:,:,:,3) =tabres(:,:,:,3)*Agrif_Rhoy()1249 tabres(:,:,:,1)=tabres(:,:,:,1)*Agrif_Rhox()*Agrif_Rhoy() 1250 tabres(:,:,:,2)=tabres(:,:,:,2)*Agrif_Rhox() 1251 tabres(:,:,:,3)=tabres(:,:,:,3)*Agrif_Rhoy() 1240 1252 ELSE 1241 1253 DO jk=k1,k2 … … 1246 1258 print *,'VAL2 = ',ji,jj,jk,tabres(ji,jj,jk,2),e1t(ji,jj)*tmask(ji,jj,jk) 1247 1259 print *,'VAL3 = ',ji,jj,jk,tabres(ji,jj,jk,3),e2t(ji,jj)*tmask(ji,jj,jk) 1248 ztemp = SQRT( tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3)))1260 ztemp = sqrt(tabres(ji,jj,jk,1)/(tabres(ji,jj,jk,2)*tabres(ji,jj,jk,3))) 1249 1261 print *,'CORR = ',ztemp-1. 1250 1262 print *,'NEW VALS = ',tabres(ji,jj,jk,2)*ztemp,tabres(ji,jj,jk,3)*ztemp, & 1251 1263 tabres(ji,jj,jk,2)*ztemp*tabres(ji,jj,jk,3)*ztemp 1252 e1t(ji,jj) = tabres(ji,jj,jk,2) *ztemp1253 e2t(ji,jj) = tabres(ji,jj,jk,3) *ztemp1264 e1t(ji,jj) = tabres(ji,jj,jk,2)*ztemp 1265 e2t(ji,jj) = tabres(ji,jj,jk,3)*ztemp 1254 1266 END IF 1255 1267 END DO … … 1331 1343 DO jj=j1,j2 1332 1344 DO ji=i1,i2 1333 ptab(ji,jj,jk) = e3t_0(ji,jj,jk) & 1334 & * ( 1._wp + ssh(ji,jj,Kmm_a) & 1335 & * ssmask(ji,jj) & 1336 & / ( ht_0(ji,jj)-1._wp + ssmask(ji,jj) ) ) 1345 ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + ssh(ji,jj,Kmm_a) & 1346 & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj))) 1337 1347 END DO 1338 1348 END DO … … 1353 1363 DO jj=j1,j2 1354 1364 DO ji=i1,i2 1355 e3t(ji,jj,jk,Kbb_a) = e3t(ji,jj,jk,Kbb_a)&1356 &+ atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) )1365 e3t(ji,jj,jk,Kbb_a) = e3t(ji,jj,jk,Kbb_a) & 1366 & + atfp * ( ptab(ji,jj,jk) - e3t(ji,jj,jk,Kmm_a) ) 1357 1367 END DO 1358 1368 END DO … … 1367 1377 DO ji = i1,i2 1368 1378 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 1369 e3w(ji,jj,jk,Kbb_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) 1370 & *( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) ) &1371 & + 0.5_wp * tmask(ji,jj,jk)&1372 & *( e3t(ji,jj,jk ,Kbb_a) - e3t_0(ji,jj,jk ) )1373 gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) +e3t(ji,jj,jk-1,Kbb_a)1374 gdept(ji,jj,jk,Kbb_a) = zcoef * ( gdepw(ji,jj,jk ,Kbb_a) + 0.5 * e3w(ji,jj,jk ,Kbb_a)) &1375 & + (1._wp - zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) + e3w(ji,jj,jk ,Kbb_a))1379 e3w(ji,jj,jk,Kbb_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * & 1380 & ( e3t(ji,jj,jk-1,Kbb_a) - e3t_0(ji,jj,jk-1) ) & 1381 & + 0.5_wp * tmask(ji,jj,jk) * & 1382 & ( e3t(ji,jj,jk ,Kbb_a) - e3t_0(ji,jj,jk ) ) 1383 gdepw(ji,jj,jk,Kbb_a) = gdepw(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk-1,Kbb_a) 1384 gdept(ji,jj,jk,Kbb_a) = zcoef * ( gdepw(ji,jj,jk ,Kbb_a) + 0.5 * e3w(ji,jj,jk,Kbb_a)) & 1385 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb_a) + e3w(ji,jj,jk,Kbb_a)) 1376 1386 END DO 1377 1387 END DO … … 1393 1403 ! 1394 1404 ! Update vertical scale factor at W-points and depths: 1395 e3w(i1:i2,j1:j2,1,Kmm_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm_a) - e3t_0(i1:i2,j1:j2,1)1405 e3w (i1:i2,j1:j2,1,Kmm_a) = e3w_0(i1:i2,j1:j2,1) + e3t(i1:i2,j1:j2,1,Kmm_a) - e3t_0(i1:i2,j1:j2,1) 1396 1406 gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a) 1397 1407 gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp 1398 gde3w(i1:i2,j1:j2,1 ) = gdept(i1:i2,j1:j2,1,Kmm_a) - ( ht(i1:i2,j1:j2) - ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh1408 gde3w(i1:i2,j1:j2,1) = gdept(i1:i2,j1:j2,1,Kmm_a) - (ht(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh 1399 1409 ! 1400 1410 DO jk = 2, jpk … … 1402 1412 DO ji = i1,i2 1403 1413 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 1404 e3w(ji,jj,jk,Kmm_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) & 1405 & * ( e3t(ji,jj,jk-1,Kmm_a) - e3t_0(ji,jj,jk-1) ) & 1406 & + 0.5_wp * tmask(ji,jj,jk) & 1407 & * ( e3t(ji,jj,jk ,Kmm_a) - e3t_0(ji,jj,jk ) ) 1408 gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a) 1409 gdept(ji,jj,jk,Kmm_a) = zcoef * ( gdepw(ji,jj,jk ,Kmm_a) + 0.5 * e3w(ji,jj,jk ,Kmm_a) ) & 1410 & + ( 1._wp - zcoef ) * ( gdept(ji,jj,jk-1,Kmm_a) + e3w(ji,jj,jk ,Kmm_a) ) 1411 gde3w(ji,jj,jk ) = gdept(ji,jj,jk ,Kmm_a) - ( ht(ji,jj)-ht_0(ji,jj) ) ! Last term in the rhs is ssh 1414 e3w(ji,jj,jk,Kmm_a) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t(ji,jj,jk-1,Kmm_a) - e3t_0(ji,jj,jk-1) ) & 1415 & + 0.5_wp * tmask(ji,jj,jk) * ( e3t(ji,jj,jk ,Kmm_a) - e3t_0(ji,jj,jk ) ) 1416 gdepw(ji,jj,jk,Kmm_a) = gdepw(ji,jj,jk-1,Kmm_a) + e3t(ji,jj,jk-1,Kmm_a) 1417 gdept(ji,jj,jk,Kmm_a) = zcoef * ( gdepw(ji,jj,jk ,Kmm_a) + 0.5 * e3w(ji,jj,jk,Kmm_a)) & 1418 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm_a) + e3w(ji,jj,jk,Kmm_a)) 1419 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm_a) - (ht(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh 1412 1420 END DO 1413 1421 END DO … … 1415 1423 ! 1416 1424 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 1417 e3t(i1:i2,j1:j2,1:jpk,Kbb_a) = e3t(i1:i2,j1:j2,1:jpk,Kmm_a)1418 e3w(i1:i2,j1:j2,1:jpk,Kbb_a) = e3w(i1:i2,j1:j2,1:jpk,Kmm_a)1425 e3t (i1:i2,j1:j2,1:jpk,Kbb_a) = e3t (i1:i2,j1:j2,1:jpk,Kmm_a) 1426 e3w (i1:i2,j1:j2,1:jpk,Kbb_a) = e3w (i1:i2,j1:j2,1:jpk,Kmm_a) 1419 1427 gdepw(i1:i2,j1:j2,1:jpk,Kbb_a) = gdepw(i1:i2,j1:j2,1:jpk,Kmm_a) 1420 1428 gdept(i1:i2,j1:j2,1:jpk,Kbb_a) = gdept(i1:i2,j1:j2,1:jpk,Kmm_a)
Note: See TracChangeset
for help on using the changeset viewer.
