Changeset 8822


Ignore:
Timestamp:
2017-11-27T14:55:00+01:00 (3 years ago)
Author:
timgraham
Message:

Bug fixes in interp and update.
Modified update to deal with none matching depths at U/V points

Location:
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90

    r8690 r8822  
    665665             IF (ptab(ji,jj,jk,n2) == 0) EXIT 
    666666             N_in = N_in + 1 
    667              tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/ptab(ji,jj,jk,n2) 
    668              h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)) 
     667             tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/(ptab(ji,jj,jk,n2)) 
     668             h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy) 
    669669           END DO 
    670670           N_out = 0 
     
    716716 
    717717         ! 
    718          western_side  = (nb == 1).AND.(ndir == 1) 
    719          eastern_side  = (nb == 1).AND.(ndir == 2) 
    720          southern_side = (nb == 2).AND.(ndir == 1) 
    721          northern_side = (nb == 2).AND.(ndir == 2) 
    722          ! 
    723718         zrhox = Agrif_Rhox() 
    724719         !  
     
    939934                  N_in = N_in + 1 
    940935                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 
    941                   h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj))  
     936                  h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)  
    942937              ENDDO 
    943938          
     
    968963              ENDIF 
    969964              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
    970           
    971965            ENDDO 
    972966         ENDDO 
     
    10501044                  N_in = N_in + 1 
    10511045                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 
    1052                   h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)) 
     1046                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 
    10531047               enddo 
    10541048               IF (N_in == 0) THEN 
  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90

    r8690 r8822  
    206206                  DO jj=j1,j2 
    207207                     DO ji=i1,i2 
    208                         tabres(ji,jj,jk,jn) = (zrho_xy * tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) & 
    209                                               * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp 
     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 
    210210                     END DO 
    211211                  END DO 
     
    215215               DO jj=j1,j2 
    216216                  DO ji=i1,i2 
    217                      tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) & 
     217                     tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) * zrho_xy & 
    218218                                              + (tmask(ji,jj,jk)-1)*999._wp 
    219219                  END DO 
     
    240240           N_in = 0 
    241241           DO jk=k1,k2 !k2 = jpk of child grid 
    242              IF (tabres(ji,jj,jk,n2) > -900) EXIT 
     242             IF (tabres(ji,jj,jk,n2) == 0  ) EXIT 
    243243             N_in = N_in + 1 
    244244             tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     
    247247           N_out = 0 
    248248           DO jk=1,jpk ! jpk of parent grid 
    249              IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF 
     249             IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 
    250250             N_out = N_out + 1 
    251251             h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 
    252252           ENDDO 
    253 !           IF(ji.EQ.i1 .AND. jj.EQ.j1) print *,'1st parent point',sum(h_in(1:N_in)), sum(h_out(1:N_out)) 
    254253           IF (N_in > 0) THEN 
    255254             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    256 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
    257 !             IF abs(h_diff > 1.e-8) THEN 
    258 !               N_in = N_in + 1 
    259 !               h_in(N_in) = h_diff 
    260 !               tabin(N_in,:) = tabin(N_in-1,:) 
    261255             IF (h_diff < -1.e-4) THEN 
    262              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_out 
    263              print *,h_in(1:N_in) 
    264              print *,h_out(1:N_out) 
    265              STOP 
    266 !               N_out = N_out + 1 
    267 !               h_out(N_out) = - h_diff 
     256                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_out 
     257                print *,h_in(1:N_in) 
     258                print *,h_out(1:N_out) 
     259                STOP 
    268260             ENDIF 
    269261             DO jn=n1,n2-1 
     
    276268            tabres_child(:,:,:,:) = tabres(:,:,:,1:jpts) 
    277269#endif 
    278           
    279270! WARNING : 
    280271! tabres replaced by tabres_child in the following 
    281272! k1 k2 replaced by 1 jpk 
    282273! VERTICAL REFINEMENT END 
    283  
    284274         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    285275            ! Add asselin part 
     
    329319      REAL(wp) :: h_out(1:jpk) 
    330320      INTEGER :: N_in, N_out 
    331       REAL(wp) :: h_diff 
     321      REAL(wp) :: h_diff, excess, thick 
    332322      REAL(wp) :: tabin(k1:k2) 
    333323! VERTICAL REFINEMENT END 
     
    343333               DO ji=i1,i2 
    344334                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  & 
    345                                        - (umask(ji,jj,jk)-1)*999._wp 
    346                   tabres(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  & 
    347                                        - (umask(ji,jj,jk)-1)*999._wp 
     335                                       + (umask(ji,jj,jk)-1)*999._wp 
     336                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  & 
     337                                       + (umask(ji,jj,jk)-1)*999._wp 
    348338               END DO 
    349339            END DO 
     
    357347         tabres_child(:,:,:) = 0. 
    358348# if defined key_vertical 
    359          AGRIF_SpecialValue = -999._wp 
     349         AGRIF_SpecialValue = 0._wp 
    360350         DO jj=j1,j2 
    361351         DO ji=i1,i2 
    362352           N_in = 0 
     353           h_in(:) = 0._wp 
     354           tabin(:) = 0._wp 
    363355           DO jk=k1,k2 !k2=jpk of child grid 
    364              IF (tabres(ji,jj,jk,2) > -900) EXIT 
     356             IF( tabres(ji,jj,jk,2) < -900) EXIT 
    365357             N_in = N_in + 1 
    366358             tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     
    375367           IF (N_in * N_out > 0) THEN 
    376368             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    377 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
    378              if (h_diff < -1.e-4) then 
    379              print *,'CHECK YOUR bathy U points ...',ji,jj,h_diff,e2u(ji,jj),sum(h_in(1:N_in)),sum(h_out(1:N_out)), N_in, N_out 
    380              stop 
    381 !             else ! Extends with 0 
    382 !             N_in = N_in + 1 
    383 !             tabin(N_in) = 0. 
    384 !             h_in(N_in) = h_diff 
    385              endif 
     369             IF (h_diff < -1.e-4) then 
     370!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
     371!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 
    386385             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) 
    387           ENDIF 
     386             tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 
     387           ENDIF 
    388388         ENDDO 
    389          ENDDO 
     389      ENDDO 
    390390#else 
    391391            DO jk=1,jpk 
     
    397397            END DO 
    398398#endif 
    399           
    400 ! WARNING : 
    401 ! tabres replaced by tabres_child in the following 
    402 ! k1 k2 replaced by 1 jpk 
    403 ! remove division by fs e3u (already done) 
    404 ! VERTICAL REFINEMENT END 
    405399 
    406400         DO jk=1,jpk 
    407401            DO jj=j1,j2 
    408402               DO ji=i1,i2 
    409 !Following line now replaced by division higher up in vertical refinement case 
    410 !                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) 
    411                   ! 
    412403                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    413404                     ub(ji,jj,jk) = ub(ji,jj,jk) &  
     
    439430      REAL(wp) :: h_out(1:jpk) 
    440431      INTEGER :: N_in, N_out 
    441       REAL(wp) :: h_diff 
     432      REAL(wp) :: h_diff, excess, thick 
    442433      REAL(wp) :: tabin(k1:k2) 
    443434! VERTICAL REFINEMENT END 
     
    447438         zrhox = Agrif_Rhox() 
    448439#if defined key_vertical 
     440         AGRIF_SpecialValue = -999._wp 
    449441         DO jk=k1,k2 
    450442            DO jj=j1,j2 
    451443               DO ji=i1,i2 
    452                   tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) 
    453                   tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk)  
     444                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) & 
     445                                       + (vmask(ji,jj,jk)-1)*999._wp 
     446                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) & 
     447                                       + (vmask(ji,jj,jk)-1)*999._wp 
    454448               END DO 
    455449            END DO 
     
    466460      ELSE 
    467461         tabres_child(:,:,:) = 0. 
    468 ! VERTICAL REFINEMENT BEGIN 
    469462#if defined key_vertical 
     463         AGRIF_SpecialValue = 0._wp 
    470464         DO jj=j1,j2 
    471465         DO ji=i1,i2 
    472466           N_in = 0 
    473467           DO jk=k1,k2 
    474              IF (tabres(ji,jj,jk,2) == 0) EXIT 
     468             IF (tabres(ji,jj,jk,2) < -900) EXIT 
    475469             N_in = N_in + 1 
    476470             tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     
    485479           IF (N_in * N_out > 0) THEN 
    486480             h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    487 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
    488              if (h_diff < -1.e-4) then 
    489              print *,'CHECK YOUR BATHY ...' 
    490              stop 
    491 !             else ! Extends with 0 
    492 !             N_in = N_in + 1 
    493 !             tabin(N_in) = 0. 
    494 !             h_in(N_in) = h_diff 
    495              endif 
     481             IF (h_diff < -1.e-4) then 
     482!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
     483!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._wp 
     485                DO jk=N_in,1,-1 
     486                   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 + thick 
     490                   IF ( h_diff == 0) THEN 
     491                      N_in = jk 
     492                      h_in(jk) = h_in(jk) - thick 
     493                      EXIT 
     494                   ENDIF 
     495                ENDDO 
     496             ENDIF 
    496497             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)) 
    497499          ENDIF 
    498500         ENDDO 
     
    507509            END DO 
    508510#endif 
    509           
    510 ! WARNING : 
    511 ! tabres replaced by tabres_child in the following 
    512 ! k1 k2 replaced by 1 jpk 
    513 ! remove division by fs e3v (already done) 
    514 ! VERTICAL REFINEMENT END 
    515511 
    516512         DO jk=k1,jpk 
    517513            DO jj=j1,j2 
    518514               DO ji=i1,i2 
    519 !Following line now replaced by division higher up in vertical refinement case 
    520 !                  tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) 
    521515                  ! 
    522516                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     
    761755         END DO 
    762756      ENDIF 
    763       ! 
    764757   END SUBROUTINE updatevb2b 
    765758 
Note: See TracChangeset for help on using the changeset viewer.