Changeset 8902


Ignore:
Timestamp:
2017-12-05T16:58:31+01:00 (3 years ago)
Author:
timgraham
Message:

Bug fixes for non key_vertical case, reverted to interpolating tracer rather than tracer content

Location:
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO
Files:
4 edited

Legend:

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

    r8861 r8902  
    617617               DO jj=j1,j2 
    618618                 DO ji=i1,i2 
    619                        ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 
     619                       ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 
    620620                 END DO 
    621621              END DO 
     
    648648                  IF (ptab(ji,jj,jk,n2) == 0) EXIT 
    649649                  N_in = N_in + 1 
    650                   tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/(ptab(ji,jj,jk,n2)) 
     650                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 
    651651                  h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy) 
    652652               END DO 
     
    669669            DO jj=j1,j2 
    670670               DO ji=i1,i2 
    671                   ptab_child(ji,jj,jk,n1:n2-1) = ptab(ji,jj,jk,n1:n2-1)/(zrhoxy*e1e2t(ji,jj)*e3t_n(ji,jj,jk) 
     671                  ptab_child(ji,jj,jk,n1:n2-1) = ptab(ji,jj,jk,n1:n2-1) 
    672672               ENDDO 
    673673            ENDDO 
  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90

    r8893 r8902  
    216216               DO jj=j1,j2 
    217217                  DO ji=i1,i2 
    218                      tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 
     218                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 
    219219                  END DO 
    220220               END DO 
     
    239239               IF (tabres(ji,jj,jk,n2) == 0) EXIT 
    240240               N_in = N_in + 1 
    241                tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     241               tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 
    242242               h_in(N_in) = tabres(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy) 
    243243             END DO 
     
    264264            do jj=j1,j2 
    265265               do ji=i1,i2 
    266                   tsbdiff(ji,jj,jk,n1:n2-1) = tsb(ji,jj,jk,n1:n2) - tmask(ji,jj,jk)*tabres(ji,jj,jk,n1:n2-1)/e3t_n(ji,jj,jk) 
     266                  tsbdiff(ji,jj,jk,n1:n2-1) = tsb(ji,jj,jk,n1:n2) - tmask(ji,jj,jk)*tabres(ji,jj,jk,n1:n2-1) 
    267267               enddo 
    268268            enddo 
     
    397397         END DO          
    398398#else 
    399          ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:) 
     399         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:,1))*umask(i1:i2,j1:j2,:) 
    400400#endif 
    401401         ! 
     
    506506         END DO 
    507507#else 
    508          tabres = vn(i1:i2,j1:j2,:) 
     508         tabres(i1:i2,j1:j2,:,m1) = vn(i1:i2,j1:j2,:) 
    509509#endif 
    510510      ELSE 
  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90

    r8822 r8902  
    187187      !! 
    188188      INTEGER :: ji,jj,jk,jn 
    189 ! VERTICAL REFINEMENT BEGIN 
    190189      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 
    191190      REAL(wp) :: h_in(k1:k2) 
     
    195194      REAL(wp) :: zrho_xy 
    196195      REAL(wp) :: tabin(k1:k2,n1:n2) 
    197 ! VERTICAL REFINEMENT END 
    198196      !!--------------------------------------------- 
    199197      ! 
    200198      IF (before) THEN 
    201199# 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 
    214203            DO jk=k1,k2 
    215204               DO jj=j1,j2 
    216205                  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                                              + (tmask(ji,jj,jk)-1)*999._wp 
     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 
    219208                  END DO 
    220209               END DO 
    221210            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 
    222220#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) 
    229226                  END DO 
    230227               END DO 
    231             END DO         
     228            END DO 
     229         END DO         
    232230#endif 
    233231      ELSE 
    234 ! VERTICAL REFINEMENT BEGIN 
    235232         tabres_child(:,:,:,:) = 0. 
    236233# if defined key_vertical 
    237234         AGRIF_SpecialValue = 0._wp 
    238235         DO jj=j1,j2 
    239          DO ji=i1,i2 
    240            N_in = 0 
    241            DO jk=k1,k2 !k2 = jpk of child grid 
    242              IF (tabres(ji,jj,jk,n2) == 0  ) EXIT 
    243              N_in = N_in + 1 
    244              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            ENDDO 
    247            N_out = 0 
    248            DO jk=1,jpk ! jpk of parent grid 
    249              IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 
    250              N_out = N_out + 1 
    251              h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 
    252            ENDDO 
    253            IF (N_in > 0) THEN 
    254              h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    255              IF (h_diff < -1.e-4) THEN 
    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 
    260              ENDIF 
    261              DO jn=n1,n2-1 
    262                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              ENDDO 
    264           ENDIF 
    265          ENDDO 
     236            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 
    266263         ENDDO 
    267264#else 
    268             tabres_child(:,:,:,:) = tabres(:,:,:,1:jpts) 
     265         tabres_child(:,:,:,:) = tabres(:,:,:,1:jpts) 
    269266#endif 
    270 ! WARNING : 
    271 ! tabres replaced by tabres_child in the following 
    272 ! k1 k2 replaced by 1 jpk 
    273 ! VERTICAL REFINEMENT END 
    274267         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    275268            ! Add asselin part 
     
    281274                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
    282275                                 & + atfp * ( tabres_child(ji,jj,jk,jn) & 
    283                                  &             - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
     276                                 &          - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
    284277                        ENDIF 
    285278                     ENDDO 
     
    325318      !  
    326319      IF( before ) THEN 
    327          print *, i1,i2,j1,j2,k1,k2 
    328320         zrhoy = Agrif_Rhoy() 
    329321# if defined key_vertical 
     
    340332         END DO 
    341333#else 
    342             DO jk = k1,k2 
    343                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) 
    344             END DO 
     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 
    345337#endif 
    346338      ELSE 
     
    349341         AGRIF_SpecialValue = 0._wp 
    350342         DO jj=j1,j2 
    351          DO ji=i1,i2 
    352            N_in = 0 
    353            h_in(:) = 0._wp 
    354            tabin(:) = 0._wp 
    355            DO jk=k1,k2 !k2=jpk of child grid 
    356              IF( tabres(ji,jj,jk,2) < -900) EXIT 
    357              N_in = N_in + 1 
    358              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            ENDDO 
    361            N_out = 0 
    362            DO jk=1,jpk 
    363              IF (umask(ji,jj,jk) == 0) EXIT 
    364              N_out = N_out + 1 
    365              h_out(N_out) = e3u_n(ji,jj,jk) 
    366            ENDDO 
    367            IF (N_in * N_out > 0) THEN 
    368              h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    369              IF (h_diff < -1.e-4) then 
     343            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 
    370362!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
    371363!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 
    388381         ENDDO 
    389       ENDDO 
    390382#else 
    391383            DO jk=1,jpk 
     
    463455         AGRIF_SpecialValue = 0._wp 
    464456         DO jj=j1,j2 
    465          DO ji=i1,i2 
    466            N_in = 0 
    467            DO jk=k1,k2 
    468              IF (tabres(ji,jj,jk,2) < -900) EXIT 
    469              N_in = N_in + 1 
    470              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            ENDDO 
    473            N_out = 0 
    474            DO jk=1,jpk 
    475              IF (vmask(ji,jj,jk) == 0) EXIT 
    476              N_out = N_out + 1 
    477              h_out(N_out) = e3v_n(ji,jj,jk) 
    478            ENDDO 
    479            IF (N_in * N_out > 0) THEN 
    480              h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    481              IF (h_diff < -1.e-4) then 
     457            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 
    482474!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
    483475!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 
    497              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           ENDIF 
    500          ENDDO 
     476                     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 
    501493         ENDDO 
    502494#else 
    503             DO jk=k1,k2 
    504                DO jj=j1,j2 
    505                   DO ji=i1,i2 
    506                      tabres(ji,jj,jk,1) = e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
    507                   END DO 
    508                END DO 
    509             END DO 
     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 
    510502#endif 
    511503 
  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

    r7761 r8902  
    372372      ! Save maximum number of levels in jpkglo, then define all vertical grids with this number. 
    373373      ! Suppress once vertical online interpolation is ok 
    374       IF(.NOT.Agrif_Root())   jpkglo = Agrif_Parent( jpkglo ) 
     374!      IF(.NOT.Agrif_Root())   jpkglo = Agrif_Parent( jpkglo ) 
    375375#endif 
    376376      jpim1 = jpi-1                                            ! inner domain indices 
Note: See TracChangeset for help on using the changeset viewer.