New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 13565 for NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST – NEMO

Ignore:
Timestamp:
2020-10-05T16:18:53+02:00 (4 years ago)
Author:
jchanut
Message:

#2222, 1) Added parent bathymetry volume consistency check 2) Added velocity extrapolation in update 3) Corrected bdy issue #2519

Location:
NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce.F90

    r13351 r13565  
    7070   INTEGER, PUBLIC :: mbkt_id, ht0_id 
    7171   INTEGER, PUBLIC :: glamt_id, gphit_id 
     72   INTEGER, PUBLIC :: batupd_id 
    7273   INTEGER, PUBLIC :: kindic_agr 
    7374 
  • NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_sponge.F90

    r13498 r13565  
    137137          
    138138         ztabramp(:,:) = 0._wp 
    139  
    140          ! Trick to remove sponge in 2DV domains: 
    141          IF ( nbcellsx <= 3 ) ispongearea = -1 
    142          IF ( nbcellsy <= 3 ) jspongearea = -1 
    143139 
    144140         IF( lk_west ) THEN                             ! --- West --- ! 
  • NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_update.F90

    r13351 r13565  
    2121   USE zdf_oce        ! vertical physics: ocean variables  
    2222   USE agrif_oce 
     23   USE dom_oce 
    2324   ! 
    2425   USE in_out_manager ! I/O manager 
     
    3233 
    3334   PUBLIC   Agrif_Update_Tra, Agrif_Update_Dyn, Agrif_Update_vvl, Agrif_Update_ssh 
    34    PUBLIC   Update_Scales 
     35   PUBLIC   Update_Scales, Agrif_Check_parent_bat 
    3536 
    3637   !!---------------------------------------------------------------------- 
     
    5051      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed() 
    5152 
    52       Agrif_UseSpecialValueInUpdate = .NOT.l_vremap 
     53      Agrif_UseSpecialValueInUpdate = .NOT.ln_vert_remap 
    5354      Agrif_SpecialValueFineGrid    = 0._wp 
    5455      l_vremap                      = ln_vert_remap 
     
    343344                  N_in = 0 
    344345                  DO jk=k1,k2 !k2 = jpk of child grid 
    345                      IF (tabres(ji,jj,jk,n2) == 0._wp  ) EXIT 
     346                     IF (tabres(ji,jj,jk,n2) <= 1.e-6_wp  ) EXIT 
    346347                     N_in = N_in + 1 
    347348                     tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     
    448449      REAL(wp) :: h_in(k1:k2) 
    449450      REAL(wp) :: h_out(1:jpk) 
    450       INTEGER  :: N_in, N_out 
    451       REAL(wp) :: h_diff, excess, thick 
     451      INTEGER  :: N_in, N_out, N_in_save, N_out_save 
     452      REAL(wp) :: zhmin, zd 
    452453      REAL(wp) :: tabin(k1:k2) 
    453454! VERTICAL REFINEMENT END 
     
    470471 
    471472         tabres_child(:,:,:) = 0._wp 
    472          AGRIF_SpecialValue = 0._wp 
    473473 
    474474         IF ( l_vremap ) THEN 
     
    480480                  tabin(:) = 0._wp 
    481481                  DO jk=k1,k2 !k2=jpk of child grid 
    482                      IF( tabres(ji,jj,jk,2) == 0.) EXIT 
     482                     IF( tabres(ji,jj,jk,2)*r1_e2u(ji,jj) <= 1.e-6_wp ) EXIT 
    483483                     N_in = N_in + 1 
    484484                     tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     
    487487                  N_out = 0 
    488488                  DO jk=1,jpk 
    489                      IF (umask(ji,jj,jk) == 0) EXIT 
     489                     IF (umask(ji,jj,jk) == 0._wp) EXIT 
    490490                     N_out = N_out + 1 
    491491                     h_out(N_out) = e3u(ji,jj,jk,Kmm_a) 
    492492                  ENDDO 
    493493                  IF (N_in * N_out > 0) THEN 
    494                      h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    495                      excess = 0._wp 
    496                      IF (h_diff < -1.e-4) THEN 
    497                         DO jk=N_in,1,-1 
    498                            thick = MIN(-1*h_diff, h_in(jk)) 
    499                            excess = excess + tabin(jk)*thick*e2u(ji,jj) 
    500                            tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 
    501                            h_diff = h_diff + thick 
    502                            IF ( h_diff == 0) THEN 
     494                     ! Deal with potentially different depths at velocity points: 
     495                     N_in_save  = N_in 
     496                     N_out_save = N_out 
     497                     IF ( ABS(sum(h_out(1:N_out))-sum(h_in(1:N_in))) > 1.e-6_wp ) THEN 
     498                        zhmin = MIN(sum(h_out(1:N_out)), sum(h_in(1:N_in))) 
     499                        zd = 0._wp 
     500                        DO jk=1, N_in_save 
     501                           IF ( (zd +  h_in(jk)) > zhmin-1.e-6) THEN 
    503502                              N_in = jk 
    504                               h_in(jk) = h_in(jk) - thick 
    505                               EXIT 
     503                              h_in(jk) = zhmin - zd 
     504                              EXIT  
    506505                           ENDIF 
    507                         ENDDO 
    508                      ENDIF 
     506                           zd = zd + h_in(jk) 
     507                        END DO 
     508                        zd = 0._wp 
     509                        DO jk=1, N_out_save 
     510                           IF ( (zd +  h_out(jk)) > zhmin-1.e-6) THEN 
     511                              N_out = jk 
     512                              h_out(jk) = zhmin - zd 
     513                              EXIT  
     514                           ENDIF 
     515                           zd = zd + h_out(jk) 
     516                        END DO 
     517                     END IF 
    509518                     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) 
    510                      tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 
     519                     IF (N_out < N_out_save) tabres_child(ji,jj,N_out+1:N_out_save) = tabres_child(ji,jj,N_out) 
    511520                  ENDIF 
    512521               ENDDO 
     
    606615      REAL(wp) :: h_in(k1:k2) 
    607616      REAL(wp) :: h_out(1:jpk) 
    608       INTEGER :: N_in, N_out 
    609       REAL(wp) :: h_diff, excess, thick 
     617      INTEGER  :: N_in, N_out, N_in_save, N_out_save 
     618      REAL(wp) :: zhmin, zd 
    610619      REAL(wp) :: tabin(k1:k2) 
    611620! VERTICAL REFINEMENT END 
     
    628637 
    629638         tabres_child(:,:,:) = 0._wp 
    630          AGRIF_SpecialValue = 0._wp 
    631639 
    632640         IF ( l_vremap ) THEN 
     
    636644                  N_in = 0 
    637645                  DO jk=k1,k2 
    638                      IF (tabres(ji,jj,jk,2) == 0) EXIT 
     646                     IF (tabres(ji,jj,jk,2)* r1_e1v(ji,jj) <= 1.e-6_wp) EXIT 
    639647                     N_in = N_in + 1 
    640648                     tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     
    648656                  ENDDO 
    649657                  IF (N_in * N_out > 0) THEN 
    650                      h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    651                      excess = 0._wp 
    652                      IF (h_diff < -1.e-4) then 
    653 !Even if bathy at T points match it's possible for the V points to be deeper in the child grid.  
    654 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 
    655                         DO jk=N_in,1,-1 
    656                            thick = MIN(-1*h_diff, h_in(jk)) 
    657                            excess = excess + tabin(jk)*thick*e2u(ji,jj) 
    658                            tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 
    659                            h_diff = h_diff + thick 
    660                            IF ( h_diff == 0) THEN 
     658                     ! Deal with potentially different depths at velocity points: 
     659                     N_in_save  = N_in 
     660                     N_out_save = N_out 
     661                     IF ( ABS(sum(h_out(1:N_out))-sum(h_in(1:N_in))) > 1.e-6_wp ) THEN 
     662                        zhmin = MIN(sum(h_out(1:N_out)), sum(h_in(1:N_in))) 
     663                        zd = 0._wp 
     664                        DO jk=1, N_in_save 
     665                           IF ( (zd +  h_in(jk)) > zhmin-1.e-6) THEN 
    661666                              N_in = jk 
    662                               h_in(jk) = h_in(jk) - thick 
    663                               EXIT 
     667                              h_in(jk) = zhmin - zd 
     668                              EXIT  
    664669                           ENDIF 
    665                         ENDDO 
    666                      ENDIF 
     670                           zd = zd + h_in(jk) 
     671                        END DO 
     672                        zd = 0._wp 
     673                        DO jk=1, N_out_save 
     674                           IF ( (zd +  h_out(jk)) > zhmin-1.e-6) THEN 
     675                              N_out = jk 
     676                              h_out(jk) = zhmin - zd 
     677                              EXIT  
     678                           ENDIF 
     679                           zd = zd + h_out(jk) 
     680                        END DO 
     681                     END IF 
    667682                     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) 
    668                      tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 
     683                     IF (N_out < N_out_save) tabres_child(ji,jj,N_out+1:N_out_save) = tabres_child(ji,jj,N_out) 
    669684                  ENDIF 
    670685               ENDDO 
     
    13161331   END SUBROUTINE updatee3t 
    13171332 
     1333   SUBROUTINE Agrif_Check_parent_bat( ) 
     1334      !!---------------------------------------------------------------------- 
     1335      !!                   *** ROUTINE Agrif_Check_parent_bat *** 
     1336      !!---------------------------------------------------------------------- 
     1337      !  
     1338      IF (( .NOT.ln_agrif_2way ).OR.(.NOT.ln_chk_bathy).OR.(Agrif_Root())) RETURN 
     1339      ! 
     1340      Agrif_UseSpecialValueInUpdate = .FALSE. 
     1341      ! 
     1342      IF(lwp) WRITE(numout,*) ' ' 
     1343      IF(lwp) WRITE(numout,*) 'AGRIF: Check parent volume at Level:', Agrif_Level() 
     1344      ! 
     1345# if ! defined DECAL_FEEDBACK && ! defined DECAL_FEEDBACK_2D 
     1346      CALL Agrif_Update_Variable(batupd_id,procname = update_bat) 
     1347# else 
     1348      CALL Agrif_Update_Variable(batupd_id,locupdate=(/1,0/),procname = update_bat) 
     1349# endif 
     1350      ! 
     1351      kindic_agr = Agrif_Parent(kindic_agr) 
     1352      CALL mpp_sum( 'Agrif_Check_parent_bat', kindic_agr ) 
     1353 
     1354      IF( kindic_agr /= 0 ) THEN 
     1355         CALL ctl_stop('==> Averaged Bathymetry does not match parent volume')  
     1356      ELSE 
     1357         IF(lwp) WRITE(numout,*) '==> Averaged Bathymetry matches parent '  
     1358         IF(lwp) WRITE(numout,*) '' 
     1359      ENDIF 
     1360      ! 
     1361   END SUBROUTINE Agrif_Check_parent_bat 
     1362 
     1363   SUBROUTINE update_bat(ptab, i1, i2, j1, j2, before ) 
     1364      !!--------------------------------------------- 
     1365      !!           *** ROUTINE update_bat *** 
     1366      !!--------------------------------------------- 
     1367      REAL(wp), DIMENSION(i1:i2,j1:j2) :: ptab 
     1368      INTEGER, INTENT(in) :: i1, i2, j1, j2 
     1369      LOGICAL, INTENT(in) :: before 
     1370      INTEGER :: ji, jj 
     1371      ! 
     1372      !!--------------------------------------------- 
     1373      ! 
     1374      IF( before ) THEN 
     1375         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 
     1376      ELSE 
     1377         kindic_agr = 0 
     1378         ! 
     1379         DO jj=j1,j2 
     1380            DO ji=i1,i2 
     1381               IF ( (ssmask(ji,jj).NE.0._wp).AND.& 
     1382               & (ABS(ptab(ji,jj)-ht_0(ji,jj)).GE.1.e-6) ) THEN  
     1383                  kindic_agr = kindic_agr + 1  
     1384               ENDIF 
     1385            END DO 
     1386         END DO 
     1387         ! 
     1388      ENDIF 
     1389      ! 
     1390   END SUBROUTINE update_bat 
     1391 
    13181392#else 
    13191393   !!---------------------------------------------------------------------- 
  • NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_user.F90

    r13371 r13565  
    9191      CALL agrif_declare_variable((/2,2    /),(/ind2  ,ind3      /),(/'x','y'        /),(/1,1    /),(/jpi,jpj        /),sshini_id) 
    9292      !  
     93      ! Update location 
     94      CALL agrif_declare_variable((/2,2/),(/ind2  ,ind3  /),(/'x','y'/),(/1,1/),(/jpi,jpj/), batupd_id) 
    9395      
    9496      ! 2. Type of interpolation 
     
    138140      CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Average       , update2=Agrif_Update_Full_Weighting) 
    139141      CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Full_Weighting, update2=Agrif_Update_Average       ) 
     142      CALL Agrif_Set_Updatetype(batupd_id, update = Agrif_Update_Full_Weighting) 
    140143#else 
    141144      CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Copy          , update2=Agrif_Update_Average       ) 
    142145      CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Average       , update2=Agrif_Update_Copy          ) 
    143 #endif 
    144        
     146      CALL Agrif_Set_Updatetype(batupd_id, update = Agrif_Update_Average) 
     147#endif       
     148 
    145149   !   CALL Agrif_Set_ExternalMapping(nemo_mapping) 
    146150      ! 
     
    199203      IF ( ln_sco.AND.Agrif_Parent(ln_sco) ) THEN  
    200204         DO_2D( 1, 0, 1, 0 ) 
    201             hu0_parent(ji,jj) = 0.5_wp * ( ht0_parent(ji,jj)+ht0_parent(ji+1,jj) ) 
    202             hv0_parent(ji,jj) = 0.5_wp * ( ht0_parent(ji,jj)+ht0_parent(ji,jj+1) ) 
     205            hu0_parent(ji,jj) = 0.5_wp * ( ht0_parent(ji,jj)+ht0_parent(ji+1,jj) ) * ssumask(ji,jj) 
     206            hv0_parent(ji,jj) = 0.5_wp * ( ht0_parent(ji,jj)+ht0_parent(ji,jj+1) ) * ssvmask(ji,jj) 
    203207         END_2D 
    204208      ELSE 
     
    432436! 
    433437! > Divergence conserving alternative: 
     438!      CALL Agrif_Set_bcinterp( ts_interp_id,interp =AGRIF_constant) 
     439!      CALL Agrif_Set_bcinterp( un_interp_id,interp1=Agrif_linear,interp2=AGRIF_constant   ) 
     440!      CALL Agrif_Set_bcinterp( vn_interp_id,interp1=AGRIF_constant   ,interp2=Agrif_linear) 
     441! 
     442!      CALL Agrif_Set_bcinterp(  ts_sponge_id,interp =AGRIF_constant) 
     443!      CALL Agrif_Set_bcinterp(  un_sponge_id,interp1=Agrif_linear,interp2=AGRIF_constant   ) 
     444!      CALL Agrif_Set_bcinterp(  vn_sponge_id,interp1=AGRIF_constant   ,interp2=Agrif_linear) 
     445! 
    434446!      CALL Agrif_Set_bcinterp(sshn_id,interp=AGRIF_constant) 
    435447!      CALL Agrif_Set_bcinterp(unb_id,interp1=Agrif_linear,interp2=AGRIF_constant) 
     
    785797      ENDIF 
    786798 
     799! JC => side effects of lines below to be checked: 
    787800      lk_west  = .NOT. ( Agrif_Ix() == 1 ) 
    788801      lk_east  = .NOT. ( Agrif_Ix() + nbcellsx/AGRIF_Irhox() == Agrif_Parent(jpiglo) -1 ) 
    789802      lk_south = .NOT. ( Agrif_Iy() == 1 ) 
    790803      lk_north = .NOT. ( Agrif_Iy() + nbcellsy/AGRIF_Irhoy() == Agrif_Parent(jpjglo) -1 ) 
    791  
    792804      ! 
    793805      ! Set the number of ghost cells according to periodicity 
     
    798810      IF(   jperio == 1  )   nbghostcells_x   = 0 
    799811      IF( .NOT. lk_south )   nbghostcells_y_s = 0 
     812      ! For 2DV domains: 
     813      IF (( nbcellsy <= 3 ).AND.(AGRIF_Irhoy()==1)) THEN 
     814         lk_north  = .FALSE. ; lk_south = .FALSE. 
     815         nbghostcells_y_s = nbghostcells 
     816      ENDIF 
     817      IF (( nbcellsx <= 3 ).AND.(AGRIF_Irhox()==1)) THEN 
     818         lk_east  = .FALSE. ; lk_north = .FALSE. 
     819      ENDIF 
    800820      ! Some checks 
    801821      IF( jpiglo /= nbcellsx + 2 + 2*nn_hls + nbghostcells_x   + nbghostcells_x   )   CALL ctl_stop( 'STOP',    & 
Note: See TracChangeset for help on using the changeset viewer.