Changeset 8835


Ignore:
Timestamp:
2017-11-28T15:16:43+01:00 (3 years ago)
Author:
timgraham
Message:

Finished sponge layer for tsn, U and V

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_sponge.F90

    r7646 r8835  
    1 #define SPONGE && define SPONGE_TOP 
     1#undef SPONGE && undef SPONGE_TOP 
    22 
    33MODULE agrif_opa_sponge 
     
    5959      REAL(wp) :: timecoeff 
    6060      !!--------------------------------------------- 
    61  
    6261#if defined SPONGE 
    6362      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot() 
     
    203202      INTEGER  ::   iku, ikv 
    204203      REAL(wp) :: ztsa, zabe1, zabe2, zbtr 
    205       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ztu, ztv 
    206       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff 
     204      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 
     205      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff 
     206      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child 
     207      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 
     208      REAL(wp) :: h_in(k1:k2) 
     209      REAL(wp) :: h_out(1:jpk) 
     210      INTEGER :: N_in, N_out 
     211      REAL(wp) :: h_diff, zrhoxy 
    207212      ! 
    208213      IF( before ) THEN 
    209          tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 
     214         DO jn = n1,n2-1 
     215            DO jk=k1,k2 
     216               DO jj=j1,j2 
     217                  DO ji=i1,i2 
     218                     tabres(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 
     219                  END DO 
     220               END DO 
     221            END DO 
     222         END DO 
     223         DO jk=k1,k2 
     224            DO jj=j1,j2 
     225               DO ji=i1,i2 
     226                   tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)  
     227               END DO 
     228            END DO 
     229         END DO 
    210230      ELSE    
    211231         ! 
    212          tsbdiff(:,:,:,:) = tsb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:)     
     232#ifdef key_vertical 
     233         tabres_child(:,:,:,:) = 0. 
     234         zrhoxy = Agrif_rhox()*Agrif_rhoy() 
     235         DO jj=j1,j2 
     236           DO ji=i1,i2 
     237             N_in = 0 
     238             DO jk=k1,k2 !k2 = jpk of parent 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)*zrhoxy) 
     243             END DO 
     244             N_out = 0 
     245             DO jk=1,jpk ! jpk of child grid 
     246               IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF. !This doesn't seem to work at the moment in GYRE but is OK in overflow model 
     247               N_out = N_out + 1 
     248               h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
     249             ENDDO 
     250             IF (N_in > 0) THEN 
     251                h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     252                tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for????? 
     253                DO jn=1,jpts 
     254                   call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
     255                ENDDO 
     256             ENDIF 
     257             DO jk=1,jpk 
     258                tsbdiff(ji,jj,jk,n1:n2-1) = tsb(ji, jj,jk,n1:n2-1) - tabres_child(ji,jj,jk,n1:n2-1) 
     259             ENDDO 
     260           ENDDO 
     261         ENDDO 
     262#else 
     263         do jk=k1,k2 
     264            do jj=j1,j2 
     265               do ji=i1,i2 
     266                 ! This will be slow - Is there a way to speed it up and avoid divide by zero? 
     267                 IF (tabres(ji,jj,jk,n2) .NE. 0) THEN 
     268                    tsbdiff(ji,jj,jk,n1:n2-1) = tsb(ji,jj,jk,n1:n2) - tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     269                 ELSE 
     270                    tsbdiff(ji,jj,jk,n1:n2-1) = 0._wp 
     271                 ENDIF 
     272               enddo 
     273            enddo 
     274         enddo 
     275#endif 
    213276         DO jn = 1, jpts             
    214277            DO jk = 1, jpkm1 
     
    258321 
    259322 
    260    SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before) 
     323   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before) 
    261324      !!--------------------------------------------- 
    262325      !!   *** ROUTINE interpun_sponge *** 
    263326      !!---------------------------------------------     
    264       INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 
    265       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 
     327      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 
     328      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres 
    266329      LOGICAL, INTENT(in) :: before 
    267330 
    268       INTEGER :: ji,jj,jk 
     331      INTEGER :: ji,jj,jk,N_in,N_out 
    269332 
    270333      ! sponge parameters  
    271       REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 
    272       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff 
    273       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff 
     334      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff, zrhoy 
     335      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff 
     336      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff 
     337      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
     338      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
     339      REAL(wp), DIMENSION(1:jpk) :: h_out 
    274340      INTEGER :: jmax 
    275341      !!---------------------------------------------     
    276342      ! 
    277343      IF( before ) THEN 
    278          tabres = un(i1:i2,j1:j2,:) 
     344#ifdef key_vertical 
     345         DO jk=1,jpk 
     346            DO jj=j1,j2 
     347               DO ji=i1,i2 
     348                  tabres(ji,jj,jk,m1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk))  
     349                  tabres(ji,jj,jk,m2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)) 
     350               END DO 
     351            END DO 
     352         END DO 
     353 
     354#else 
     355         tabres(i1:i2,j1:j2,:,m1) = un(i1:i2,j1:j2,:) 
     356#endif 
    279357      ELSE 
     358#ifdef key_vertical 
     359         tabres_child(:,:,:) = 0._wp 
     360         zrhoy = Agrif_rhoy() 
     361         DO jj=j1,j2 
     362            DO ji=i1,i2 
     363               N_in = 0 
     364               DO jk=k1,k2 
     365                  IF (tabres(ji,jj,jk,m2) == 0) EXIT 
     366                  N_in = N_in + 1 
     367                  tabin(jk) = tabres(ji,jj,jk,m1)/tabres(ji,jj,jk,m2) 
     368                  h_in(N_in) = tabres(ji,jj,jk,m2)/(e2u(ji,jj)*zrhoy)  
     369              ENDDO 
     370              ! 
     371              IF (N_in == 0) THEN 
     372                 tabres_child(ji,jj,:) = 0. 
     373                 CYCLE 
     374              ENDIF 
     375          
     376              N_out = 0 
     377              DO jk=1,jpk 
     378                 if (umask(ji,jj,jk) == 0) EXIT 
     379                 N_out = N_out + 1 
     380                 h_out(N_out) = e3u_n(ji,jj,jk) 
     381              ENDDO 
     382          
     383              IF (N_out == 0) THEN 
     384                 tabres_child(ji,jj,:) = 0. 
     385                 CYCLE 
     386              ENDIF 
     387          
     388              IF (N_in * N_out > 0) THEN 
     389                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     390                 if (h_diff < -1.e4) then 
     391                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 
     392                 endif 
     393              ENDIF 
     394              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) 
     395          
     396            ENDDO 
     397         ENDDO 
     398         DO jk = 1, jpkm1 
     399            DO jj=j1,j2 
     400               ubdiff(i1:i2,jj,jk) = (ub(i1:i2,jj,jk) - tabres_child(i1:i2,jj,jk))*umask(i1:i2,jj,jk) 
     401            END DO 
     402         END DO          
     403#else 
    280404         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:) 
     405#endif 
    281406         ! 
    282407         DO jk = 1, jpkm1                                 ! Horizontal slab 
     
    356481 
    357482 
    358    SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir) 
     483   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir) 
    359484      !!--------------------------------------------- 
    360485      !!   *** ROUTINE interpvn_sponge *** 
    361486      !!---------------------------------------------  
    362       INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 
    363       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 
     487      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 
     488      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres 
    364489      LOGICAL, INTENT(in) :: before 
    365490      INTEGER, INTENT(in) :: nb , ndir 
    366491      ! 
    367       INTEGER  ::   ji, jj, jk 
    368       REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr 
    369       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff 
    370       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff 
     492      INTEGER  ::   ji, jj, jk, N_in, N_out 
     493      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, h_diff, zrhox 
     494      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff 
     495      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff 
     496      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
     497      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
     498      REAL(wp), DIMENSION(1:jpk) :: h_out 
    371499      INTEGER :: imax 
    372500      !!---------------------------------------------  
    373  
     501       
    374502      IF( before ) THEN  
     503#ifdef key_vertical 
     504         DO jk=1,jpk 
     505            DO jj=j1,j2 
     506               DO ji=i1,i2 
     507                  tabres(ji,jj,jk,m1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk)) 
     508                  tabres(ji,jj,jk,m2) = (vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)) 
     509               END DO 
     510            END DO 
     511         END DO 
     512#else 
    375513         tabres = vn(i1:i2,j1:j2,:) 
     514#endif 
    376515      ELSE 
     516#ifdef key_vertical 
     517         zrhox = Agrif_rhox() 
     518         tabres_child(:,:,:) = 0._wp 
     519         DO jj=j1,j2 
     520            DO ji=i1,i2 
     521               N_in = 0 
     522               DO jk=k1,k2 
     523                  IF (tabres(ji,jj,jk,m2) == 0) EXIT 
     524                  N_in = N_in + 1 
     525                  tabin(jk) = tabres(ji,jj,jk,m1)/tabres(ji,jj,jk,m2) 
     526                  h_in(N_in) = tabres(ji,jj,jk,m2)/(e1v(ji,jj)*zrhox)  
     527              ENDDO 
     528          
     529              IF (N_in == 0) THEN 
     530                 tabres_child(ji,jj,:) = 0. 
     531                 CYCLE 
     532              ENDIF 
     533          
     534              N_out = 0 
     535              DO jk=1,jpk 
     536                 if (umask(ji,jj,jk) == 0) EXIT 
     537                 N_out = N_out + 1 
     538                 h_out(N_out) = e3v_n(ji,jj,jk) 
     539              ENDDO 
     540          
     541              IF (N_in * N_out > 0) THEN 
     542                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     543                 if (h_diff < -1.e4) then 
     544                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 
     545                 endif 
     546              ENDIF 
     547              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) 
     548            ENDDO 
     549         ENDDO 
     550         DO jk = 1, jpkm1 
     551            DO jj=j1,j2 
     552               vbdiff(i1:i2,jj,jk) = (vb(i1:i2,jj,jk) - tabres_child(i1:i2,jj,jk))*vmask(i1:i2,jj,jk) 
     553            END DO 
     554         END DO          
     555#else 
    377556         ! 
    378557         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:) 
     558#endif 
    379559         ! 
    380560         DO jk = 1, jpkm1                                 ! Horizontal slab 
  • branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_user.F90

    r8135 r8835  
    347347   !--------------------------------------------------------------------- 
    348348   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_id) 
    349    CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_sponge_id) 
     349   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_sponge_id) 
    350350 
    351351   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id) 
     
    353353   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_update_id) 
    354354   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_update_id) 
    355    CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_sponge_id) 
    356    CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_sponge_id) 
     355   CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_sponge_id) 
     356   CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_sponge_id) 
    357357 
    358358   CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),e3t_id) 
Note: See TracChangeset for help on using the changeset viewer.