Changeset 11802


Ignore:
Timestamp:
2019-10-25T17:15:20+02:00 (10 months ago)
Author:
jchanut
Message:

#2222, 1) add linear interpolation in vremap module.
2) Switch remapping of viscosity from polynomial to linear.
3) Move to truly volume weighted averages for parent to child update.

Location:
NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_interp.F90

    r11769 r11802  
    719719               tsa(ji,jj,:,:) = 0._wp 
    720720               N_in = mbkt_parent(ji,jj) 
     721               IF ( tmask(ji,jj,1) == 0._wp) N_in = 0 
    721722               zhtot = 0._wp 
    722723               DO jk=1,N_in !k2 = jpk of parent grid 
     
    833834               N_in = mbku_parent(ji,jj) 
    834835               zhtot = 0._wp 
     836               IF ( umask(ji,jj,1) == 0._wp) N_in = 0 
    835837               DO jk=1,N_in 
    836838                  IF (jk==N_in) THEN 
     
    928930               va(ji,jj,:) = 0._wp 
    929931               N_in = mbkv_parent(ji,jj) 
     932               IF ( vmask(ji,jj,1) == 0._wp) N_in = 0 
    930933               zhtot = 0._wp 
    931934               DO jk=1,N_in 
     
    12761279      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab 
    12771280      LOGICAL                                    , INTENT(in   ) ::   before 
    1278       REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
    1279       REAL(wp), DIMENSION(1:jpk) :: h_out 
    1280       INTEGER  :: N_in, N_out, ji, jj, jk 
     1281      ! 
     1282      INTEGER  :: ji, jj, jk 
     1283      INTEGER  :: N_in, N_out 
     1284      REAL(wp), DIMENSION(k1:k2) :: tabin, z_in 
     1285      REAL(wp), DIMENSION(1:jpk) :: z_out 
    12811286      !!----------------------------------------------------------------------   
    12821287      !       
     
    12891294           END DO 
    12901295        END DO 
    1291 #ifdef key_vertical          
     1296 
     1297# if defined key_vertical 
     1298        ! Interpolate thicknesses 
     1299        ! Warning: these are masked, hence extrapolated prior interpolation. 
    12921300        DO jk=k1,k2 
    12931301           DO jj=j1,j2 
    12941302              DO ji=i1,i2 
    1295                  ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w_n(ji,jj,jk)  
     1303                  ptab(ji,jj,jk,2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
    12961304              END DO 
    12971305           END DO 
    12981306        END DO 
    1299 #endif 
     1307 
     1308        ! Extrapolate thicknesses in partial bottom cells: 
     1309        ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 
     1310        IF (ln_zps) THEN 
     1311           DO jj=j1,j2 
     1312              DO ji=i1,i2 
     1313                  jk = mbkt(ji,jj) 
     1314                  ptab(ji,jj,jk,2) = 0._wp 
     1315              END DO 
     1316           END DO            
     1317        END IF 
     1318      
     1319        ! Save ssh at last level: 
     1320        IF (.NOT.ln_linssh) THEN 
     1321           ptab(i1:i2,j1:j2,k2,2) = sshn(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1)  
     1322        ELSE 
     1323           ptab(i1:i2,j1:j2,k2,2) = 0._wp 
     1324        END IF       
     1325# endif 
    13001326      ELSE  
    13011327#ifdef key_vertical          
    1302          avm_k(i1:i2,j1:j2,1:jpk) = 0. 
    1303          DO jj=j1,j2 
    1304             DO ji=i1,i2 
    1305                N_in = 0 
    1306                DO jk=k1,k2 !k2 = jpk of parent grid 
    1307                   IF (ptab(ji,jj,jk,2) == 0) EXIT 
    1308                   N_in = N_in + 1 
    1309                   tabin(jk) = ptab(ji,jj,jk,1) 
    1310                   h_in(N_in) = ptab(ji,jj,jk,2) 
     1328         IF (ln_linssh) ptab(i1:i2,j1:j2,k2,2) = 0._wp  
     1329         avm_k(i1:i2,j1:j2,k1:k2) = 0._wp 
     1330             
     1331         DO jj = j1, j2 
     1332            DO ji =i1, i2 
     1333               N_in = mbkt_parent(ji,jj) 
     1334               IF ( tmask(ji,jj,1) == 0._wp) N_in = 0 
     1335               z_in(N_in+1) = ht0_parent(ji,jj) + ptab(ji,jj,k2,2) 
     1336               DO jk = N_in, 1, -1  ! Parent vertical grid                
     1337                     z_in(jk) = z_in(jk+1) - ptab(ji,jj,jk,2) 
     1338                    tabin(jk) = ptab(ji,jj,jk,1) 
    13111339               END DO 
    13121340               N_out = 0 
    1313                DO jk=1,jpk ! jpk of child grid 
    1314                   IF (wmask(ji,jj,jk) == 0) EXIT  
     1341          z_out(1) = 0._wp  
     1342               DO jk = 2, jpk       ! Child vertical grid 
     1343                  IF (tmask(ji,jj,jk) == 0._wp) EXIT  
    13151344                  N_out = N_out + 1 
    1316                   h_out(jk) = e3t_n(ji,jj,jk) 
     1345                  z_out(jk) = z_out(jk-1) + e3t_n(ji,jj,jk-1) 
    13171346               ENDDO 
    1318                IF (N_in > 0) THEN 
    1319                   CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out,1) 
     1347               IF (N_in*N_out > 0) THEN 
     1348                  CALL remap_linear(tabin(1:N_in),z_in(1:N_in),avm_k(ji,jj,1:N_out),z_out(1:N_out),N_in,N_out,1) 
    13201349               ENDIF 
    13211350            ENDDO 
  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90

    r11769 r11802  
    289289         spongedoneU = .TRUE. 
    290290      ENDIF 
     291 
     292#if defined key_vertical 
     293      ! Trick to vertical remove interpolation in sponge layer in case of 2DV domains: 
     294      DO jj = 2, jpjm1 
     295         DO ji = 2, jpim1 
     296            IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0 
     297            IF ((fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0 
     298! 
     299            IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0 
     300            IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0 
     301! 
     302            IF ((fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0 
     303            IF ((fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0 
     304! 
     305            IF (mbkt(ji,jj) == 0) mbkt_parent(ji,jj) = 0 
     306            IF (mbku(ji,jj) == 0) mbku_parent(ji,jj) = 0 
     307            IF (mbkv(ji,jj) == 0) mbkv_parent(ji,jj) = 0 
     308         END DO 
     309      END DO 
     310      ! 
     311      ztabramp(:,:) = REAL( mbkt_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'T', 1. ) 
     312      mbkt_parent(:,:) = MAX( NINT( ztabramp(:,:) ), 1 ) 
     313      ztabramp(:,:) = REAL( mbku_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'U', 1. ) 
     314      mbku_parent(:,:) = MAX( NINT( ztabramp(:,:) ), 1 ) 
     315      ztabramp(:,:) = REAL( mbkv_parent(:,:), wp )   ;   CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'V', 1. ) 
     316      mbkv_parent(:,:) = MAX( NINT( ztabramp(:,:) ), 1 ) 
     317#endif 
    291318      ! 
    292319#endif 
     
    366393               tabres_child(ji,jj,:,:) = 0._wp  
    367394               N_in = mbkt_parent(ji,jj) 
     395!               IF (( tmask(ji,jj,1) == 0._wp ).OR.(fspt(ji,jj)==0._wp)) N_in = 0 
    368396               zhtot = 0._wp 
    369397               DO jk=1,N_in !k2 = jpk of parent grid 
     
    530558               tabres_child(ji,jj,:) = 0._wp 
    531559               N_in = mbku_parent(ji,jj) 
     560!               IF (( umask(ji,jj,1) == 0._wp ).OR.(fspu(ji,jj)==0._wp)) N_in = 0 
    532561               zhtot = 0._wp 
    533562               DO jk=1,N_in 
     
    706735               tabres_child(ji,jj,:) = 0._wp 
    707736               N_in = mbkv_parent(ji,jj) 
     737!               IF (( vmask(ji,jj,1) == 0._wp ).OR.(fspv(ji,jj)==0._wp)) N_in = 0 
    708738               zhtot = 0._wp 
    709739               DO jk=1,N_in 
  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90

    r11791 r11802  
    4949      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed() 
    5050 
    51 ! jc_alt      Agrif_UseSpecialValueInUpdate = .FALSE. 
    52       Agrif_UseSpecialValueInUpdate = .TRUE. 
     51      Agrif_UseSpecialValueInUpdate = .FALSE. 
     52! jc_alt      Agrif_UseSpecialValueInUpdate = .TRUE. 
    5353      Agrif_SpecialValueFineGrid    = 0._wp 
    5454      !  
     
    300300               DO jj=j1,j2 
    301301                  DO ji=i1,i2 
    302                      tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) & 
    303                                          &  * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp 
    304 !jc_alt                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) 
     302!jc_alt 
     303!                     tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) & 
     304!                                         &  * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp 
     305                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) 
    305306                  END DO 
    306307               END DO 
     
    310311            DO jj=j1,j2 
    311312               DO ji=i1,i2 
    312                   tabres(ji,jj,jk,n2) =      tmask(ji,jj,jk) * e3t_n(ji,jj,jk) & 
    313                                       &   + (tmask(ji,jj,jk) - 1._wp) * 999._wp 
    314 ! jc_alt                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
     313!jc_alt 
     314!                  tabres(ji,jj,jk,n2) =      tmask(ji,jj,jk) * e3t_n(ji,jj,jk) & 
     315!                                      &   + (tmask(ji,jj,jk) - 1._wp) * 999._wp 
     316                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 
    315317               END DO 
    316318            END DO 
     
    323325               N_in = 0 
    324326               DO jk=k1,k2 !k2 = jpk of child grid 
    325                   IF (tabres(ji,jj,jk,n2) < -900._wp  ) EXIT 
    326 ! jc_alt                  IF (tabres(ji,jj,jk,n2) == 0._wp  ) EXIT 
     327! jc_alt 
     328!                  IF (tabres(ji,jj,jk,n2) < -900._wp  ) EXIT 
     329                  IF (tabres(ji,jj,jk,n2) == 0._wp  ) EXIT 
    327330                  N_in = N_in + 1 
    328331                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     
    476479      IF( before ) THEN 
    477480         zrhoy = Agrif_Rhoy() 
    478          AGRIF_SpecialValue = -999._wp 
     481!jc_alt 
     482!         AGRIF_SpecialValue = -999._wp 
    479483         DO jk=k1,k2 
    480484            DO jj=j1,j2 
    481485               DO ji=i1,i2 
    482                   tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  & 
    483                                      &  + (umask(ji,jj,jk)-1._wp)*999._wp 
    484 ! jc_alt                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)   
    485                   tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  & 
    486                                      &  + (umask(ji,jj,jk)-1._wp)*999._wp 
    487 ! jc_alt                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk) 
     486!jc_alt 
     487!                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  & 
     488!                                     &  + (umask(ji,jj,jk)-1._wp)*999._wp 
     489                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)   
     490!jc_alt 
     491!                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  & 
     492!                                     &  + (umask(ji,jj,jk)-1._wp)*999._wp 
     493                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk) 
    488494               END DO 
    489495            END DO 
     
    498504               tabin(:) = 0._wp 
    499505               DO jk=k1,k2 !k2=jpk of child grid 
    500                   IF( tabres(ji,jj,jk,2) < -900._wp) EXIT 
    501 ! jc_alt                  IF( tabres(ji,jj,jk,2) == 0.) EXIT 
     506!jc_alt 
     507!                  IF( tabres(ji,jj,jk,2) < -900._wp) EXIT 
     508                  IF( tabres(ji,jj,jk,2) == 0.) EXIT 
    502509                  N_in = N_in + 1 
    503510                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     
    512519               IF (N_in * N_out > 0) THEN 
    513520                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    514                   excess  = 0._wp 
     521                  excess = 0._wp 
    515522                  IF (h_diff < -1.e-4) THEN 
    516523!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
     
    674681      IF( before ) THEN 
    675682         zrhox = Agrif_Rhox() 
    676          AGRIF_SpecialValue = -999._wp 
     683!jc_alt 
     684!         AGRIF_SpecialValue = -999._wp 
    677685         DO jk=k1,k2 
    678686            DO jj=j1,j2 
    679687               DO ji=i1,i2 
    680                   tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) & 
    681                                      & + (vmask(ji,jj,jk)-1._wp) * 999._wp 
    682 ! jc_alt                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk)  
    683                   tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) & 
    684                                      & + (vmask(ji,jj,jk)-1._wp) * 999._wp 
    685 ! jc_alt                  tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) 
     688!jc_alt 
     689!                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) & 
     690!                                     & + (vmask(ji,jj,jk)-1._wp) * 999._wp 
     691                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk)  
     692!jc_alt 
     693!                  tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) & 
     694!                                     & + (vmask(ji,jj,jk)-1._wp) * 999._wp 
     695                  tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) 
    686696               END DO 
    687697            END DO 
     
    694704               N_in = 0 
    695705               DO jk=k1,k2 
    696                   IF (tabres(ji,jj,jk,2) < -900._wp) EXIT 
    697 ! jc_alt                  IF (tabres(ji,jj,jk,2) == 0) EXIT 
     706!jc_alt 
     707!                  IF (tabres(ji,jj,jk,2) < -900._wp) EXIT 
     708                  IF (tabres(ji,jj,jk,2) == 0) EXIT 
    698709                  N_in = N_in + 1 
    699710                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/vremap.F90

    r11741 r11802  
    2020   PRIVATE 
    2121 
    22    PUBLIC   reconstructandremap 
     22   PUBLIC   reconstructandremap, remap_linear 
    2323 
    2424   !!---------------------------------------------------------------------- 
     
    354354#endif 
    355355 
     356   SUBROUTINE remap_linear(ptin, pzin, ptout, pzout, kjpk_in, kjpk_out, kn_var) 
     357      !!---------------------------------------------------------------------- 
     358      !!                    *** ROUTINE  remap_linear *** 
     359      !! 
     360      !! ** Purpose :   Linear interpolation based on input/ouputs depths 
     361      !! 
     362      !!----------------------------------------------------------------------- 
     363      INTEGER , INTENT(in   )                      ::   kjpk_in    ! Number of input levels 
     364      INTEGER , INTENT(in   )                      ::   kjpk_out   ! Number of output levels 
     365      INTEGER , INTENT(in   )                      ::   kn_var     ! Number of variables 
     366      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in)  ::   pzin       ! Input depths 
     367      REAL(wp), INTENT(in   ), DIMENSION(kjpk_out) ::   pzout      ! Output depths 
     368      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in , kn_var) ::   ptin       ! Input data 
     369      REAL(wp), INTENT(inout), DIMENSION(kjpk_out, kn_var) ::   ptout      ! Interpolated data 
     370      ! 
     371      INTEGER  :: jkin, jkout, jn 
     372      !!-------------------------------------------------------------------- 
     373      !       
     374      DO jkout = 1, kjpk_out !  Loop over destination grid 
     375         ! 
     376         IF     ( pzout(jkout) < pzin(  1    ) ) THEN    ; ptout(jkout,1:kn_var) = ptin(    1    ,1:kn_var)          
     377         ELSEIF ( pzout(jkout) > pzin(kjpk_in) ) THEN    ; ptout(jkout,1:kn_var) = ptin( kjpk_in ,1:kn_var) 
     378         ELSEIF ( ( pzout(jkout) >= pzin(1) ).AND.( pzout(jkout) <= pzin(kjpk_in) )) THEN 
     379            DO jkin = 1, kjpk_in - 1 !  Loop over source grid 
     380               IF ( pzout(jkout) >=  pzin(jkin) ) THEN 
     381                  DO jn = 1, kn_var 
     382                     ptout(jkout,jn) =  ptin(jkin,jn) + & 
     383                                     & (pzout(jkout) - pzin(jkin)) / (pzin(jkin+1)    - pzin(jkin)) & 
     384                                     &                             * (ptin(jkin+1,jn) - ptin(jkin,jn)) 
     385                  END DO   
     386                  EXIT 
     387               ENDIF   
     388            END DO 
     389         ENDIF 
     390         ! 
     391      END DO 
     392             
     393   END SUBROUTINE remap_linear 
     394 
    356395   !!====================================================================== 
    357396!$AGRIF_END_DO_NOT_TREAT 
Note: See TracChangeset for help on using the changeset viewer.