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 11590 – NEMO

Changeset 11590


Ignore:
Timestamp:
2019-09-23T18:25:29+02:00 (4 years ago)
Author:
jchanut
Message:

#2222: 1) create remapping module (vremap) and integration of D. Engwirda piecewise polynomial recontruction package (PPR_LIB cpp key). 2) Various bug corrections with key_vertical activated.

Location:
NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce.F90

    r11574 r11590  
    1717 
    1818   PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 
    19 #if defined key_vertical 
    20    PUBLIC reconstructandremap ! remapping routine 
    21 #endif    
     19   
    2220   !                                              !!* Namelist namagrif: AGRIF parameters 
    2321   LOGICAL , PUBLIC ::   ln_agrif_2way = .TRUE.    !: activate two way nesting  
    … …  
    9593   END FUNCTION agrif_oce_alloc 
    9694 
    97 #if defined key_vertical 
    98    SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)       
    99       !!---------------------------------------------------------------------- 
    100       !!                ***  FUNCTION reconstructandremap  *** 
    101       !!---------------------------------------------------------------------- 
    102       IMPLICIT NONE 
    103       INTEGER N, Nout 
    104       REAL(wp) tabin(N), tabout(Nout) 
    105       REAL(wp) hin(N), hout(Nout) 
    106       REAL(wp) coeffremap(N,3),zwork(N,3) 
    107       REAL(wp) zwork2(N+1,3) 
    108       INTEGER jk 
    109       DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8   
    110       REAL(wp) q,q01,q02,q001,q002,q0 
    111       REAL(wp) z_win(1:N+1), z_wout(1:Nout+1) 
    112       REAL(wp),PARAMETER :: dpthin = 1.D-3 
    113       INTEGER :: k1, kbox, ktop, ka, kbot 
    114       REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 
    115  
    116       z_win(1)=0.; z_wout(1)= 0. 
    117       DO jk=1,N 
    118          z_win(jk+1)=z_win(jk)+hin(jk) 
    119       ENDDO  
    120        
    121       DO jk=1,Nout 
    122          z_wout(jk+1)=z_wout(jk)+hout(jk)        
    123       ENDDO        
    124  
    125       DO jk=2,N 
    126          zwork(jk,1)=1./(hin(jk-1)+hin(jk)) 
    127       ENDDO 
    128          
    129       DO jk=2,N-1 
    130          q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1)) 
    131          zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1) 
    132          zwork(jk,3)=q0 
    133       ENDDO        
    134       
    135       DO jk= 2,N 
    136          zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1)) 
    137       ENDDO 
    138          
    139       coeffremap(:,1) = tabin(:) 
    140   
    141       DO jk=2,N-1 
    142          q001 = hin(jk)*zwork2(jk+1,1) 
    143          q002 = hin(jk)*zwork2(jk,1)         
    144          IF (q001*q002 < 0) then 
    145             q001 = 0. 
    146             q002 = 0. 
    147          ENDIF 
    148          q=zwork(jk,2) 
    149          q01=q*zwork2(jk+1,1) 
    150          q02=q*zwork2(jk,1) 
    151          IF (abs(q001) > abs(q02)) q001 = q02 
    152          IF (abs(q002) > abs(q01)) q002 = q01 
    153          
    154          q=(q001-q002)*zwork(jk,3) 
    155          q001=q001-q*hin(jk+1) 
    156          q002=q002+q*hin(jk-1) 
    157          
    158          coeffremap(jk,3)=coeffremap(jk,1)+q001 
    159          coeffremap(jk,2)=coeffremap(jk,1)-q002 
    160          
    161          zwork2(jk,1)=(2.*q001-q002)**2 
    162          zwork2(jk,2)=(2.*q002-q001)**2 
    163       ENDDO 
    164          
    165       DO jk=1,N 
    166          IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN 
    167             coeffremap(jk,3) = coeffremap(jk,1) 
    168             coeffremap(jk,2) = coeffremap(jk,1) 
    169             zwork2(jk,1) = 0. 
    170             zwork2(jk,2) = 0. 
    171          ENDIF 
    172       ENDDO 
    173          
    174       DO jk=2,N 
    175          q002=max(zwork2(jk-1,2),dsmll) 
    176          q001=max(zwork2(jk,1),dsmll) 
    177          zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002) 
    178       ENDDO 
    179          
    180       zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) 
    181       zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) 
    182   
    183       DO jk=1,N 
    184          q01=zwork2(jk+1,3)-coeffremap(jk,1) 
    185          q02=coeffremap(jk,1)-zwork2(jk,3) 
    186          q001=2.*q01 
    187          q002=2.*q02 
    188          IF (q01*q02<0) then 
    189             q01=0. 
    190             q02=0. 
    191          ELSEIF (abs(q01)>abs(q002)) then 
    192             q01=q002 
    193          ELSEIF (abs(q02)>abs(q001)) then 
    194             q02=q001 
    195          ENDIF 
    196          coeffremap(jk,2)=coeffremap(jk,1)-q02 
    197          coeffremap(jk,3)=coeffremap(jk,1)+q01 
    198       ENDDO 
    199  
    200       zbot=0.0 
    201       kbot=1 
    202       DO jk=1,Nout 
    203          ztop=zbot  !top is bottom of previous layer 
    204          ktop=kbot 
    205          IF     (ztop.GE.z_win(ktop+1)) then 
    206             ktop=ktop+1 
    207          ENDIF 
    208          
    209          zbot=z_wout(jk+1) 
    210          zthk=zbot-ztop 
    211  
    212          IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN 
    213  
    214             kbot=ktop 
    215             DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 
    216                kbot=kbot+1 
    217             ENDDO 
    218             zbox=zbot 
    219             DO k1= jk+1,Nout 
    220                IF     (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN 
    221                   exit !thick layer 
    222                ELSE 
    223                   zbox=z_wout(k1+1)  !include thin adjacent layers 
    224                   IF(zbox.EQ.z_wout(Nout+1)) THEN 
    225                      exit !at bottom 
    226                   ENDIF 
    227                ENDIF 
    228             ENDDO 
    229             zthk=zbox-ztop 
    230  
    231             kbox=ktop 
    232             DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 
    233                kbox=kbox+1 
    234             ENDDO 
    235            
    236             IF(ktop.EQ.kbox) THEN 
    237                IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN 
    238                   IF(hin(kbox).GT.dpthin) THEN 
    239                      q001 = (zbox-z_win(kbox))/hin(kbox) 
    240                      q002 = (ztop-z_win(kbox))/hin(kbox) 
    241                      q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) 
    242                      q02=q01-1.+(q001+q002) 
    243                      q0=1.-q01-q02 
    244                   ELSE 
    245                      q0 = 1.0 
    246                      q01 = 0. 
    247                      q02 = 0. 
    248                   ENDIF 
    249                   tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 
    250                ELSE 
    251                   tabout(jk) = tabin(kbox) 
    252                ENDIF  
    253             ELSE 
    254                IF(ktop.LE.jk .AND. kbox.GE.jk) THEN 
    255                   ka = jk 
    256                ELSEIF (kbox-ktop.GE.3) THEN 
    257                   ka = (kbox+ktop)/2 
    258                ELSEIF (hin(ktop).GE.hin(kbox)) THEN 
    259                   ka = ktop 
    260                ELSE 
    261                   ka = kbox 
    262                ENDIF !choose ka 
    263      
    264                offset=coeffremap(ka,1) 
    265      
    266                qtop = z_win(ktop+1)-ztop !partial layer thickness 
    267                IF(hin(ktop).GT.dpthin) THEN 
    268                   q=(ztop-z_win(ktop))/hin(ktop) 
    269                   q01=q*(q-1.) 
    270                   q02=q01+q 
    271                   q0=1-q01-q02             
    272                ELSE 
    273                   q0 = 1. 
    274                   q01 = 0. 
    275                   q02 = 0. 
    276                ENDIF 
    277                 
    278                tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop 
    279      
    280                DO k1= ktop+1,kbox-1 
    281                   tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) 
    282                ENDDO !k1 
    283                 
    284                qbot = zbox-z_win(kbox) !partial layer thickness 
    285                IF(hin(kbox).GT.dpthin) THEN 
    286                   q=qbot/hin(kbox) 
    287                   q01=(q-1.)**2 
    288                   q02=q01-1.+q 
    289                   q0=1-q01-q02                             
    290                ELSE 
    291                   q0 = 1.0 
    292                   q01 = 0. 
    293                   q02 = 0. 
    294                ENDIF 
    295                
    296                tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot 
    297                 
    298                rpsum=1.0d0/zthk 
    299                tabout(jk)=offset+tsum*rpsum 
    300                   
    301             ENDIF !single or multiple layers 
    302          ELSE 
    303             IF (jk==1) THEN 
    304                write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1) 
    305             ENDIF 
    306             tabout(jk) = tabout(jk-1) 
    307               
    308          ENDIF !normal:thin layer 
    309       ENDDO !jk 
    310              
    311       return 
    312       end subroutine reconstructandremap 
    313 #endif 
    314  
    31595#endif 
    31696   !!====================================================================== 

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

    r11574 r11590  
    3333   USE agrif_oce_sponge 
    3434   USE lib_mpp 
     35   USE vremap 
    3536  
    3637   IMPLICIT NONE 
    … …  
    661662      REAL(wp), DIMENSION(k1:k2) :: h_in 
    662663      REAL(wp), DIMENSION(1:jpk) :: h_out 
     664      !!---------------------------------------------------------------------- 
    663665 
    664666      IF( before ) THEN          
    … …  
    698700                  IF (tmask(ji,jj,jk) == 0) EXIT  
    699701                  N_out = N_out + 1 
    700                   h_out(jk) = e3t_n(ji,jj,jk) 
     702                  h_out(jk) = e3t_a(ji,jj,jk) 
    701703               ENDDO 
    702704               IF (N_in > 0) THEN 
    703705                  DO jn=1,jpts 
    704                      call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
     706                     call reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),ptab_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 
    705707                  ENDDO 
    706708               ENDIF 
    … …  
    737739   END SUBROUTINE interpsshn 
    738740 
    739    SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir ) 
     741   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 
    740742      !!---------------------------------------------------------------------- 
    741743      !!                  *** ROUTINE interpun *** 
    … …  
    745747      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 
    746748      LOGICAL, INTENT(in) :: before 
    747       INTEGER, INTENT(in) :: nb , ndir 
    748749      !! 
    749750      INTEGER :: ji,jj,jk 
    … …  
    752753      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
    753754      REAL(wp), DIMENSION(1:jpk) :: h_out 
    754       INTEGER  :: N_in, N_out, iref 
     755      INTEGER  :: N_in, N_out 
    755756      REAL(wp) :: h_diff 
    756       LOGICAL  :: western_side, eastern_side 
    757757      !!---------------------------------------------     
    758758      ! 
    … …  
    772772# if defined key_vertical 
    773773! VERTICAL REFINEMENT BEGIN 
    774          western_side  = (nb == 1).AND.(ndir == 1) 
    775          eastern_side  = (nb == 1).AND.(ndir == 2) 
    776774 
    777775         DO ji=i1,i2 
    778             iref = ji 
    779             IF (western_side) iref = MAX(2,ji) 
    780             IF (eastern_side) iref = MIN(nlci-2,ji) 
    781776            DO jj=j1,j2 
    782777               N_in = 0 
    … …  
    795790              N_out = 0 
    796791              DO jk=1,jpk 
    797                  if (umask(iref,jj,jk) == 0) EXIT 
     792                 if (umask(ji,jj,jk) == 0) EXIT 
    798793                 N_out = N_out + 1 
    799                  h_out(N_out) = e3u_a(iref,jj,jk) 
     794                 h_out(N_out) = e3u_a(ji,jj,jk) 
    800795              ENDDO 
    801796          
    … …  
    829824   END SUBROUTINE interpun 
    830825 
    831    SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir ) 
     826   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 
    832827      !!---------------------------------------------------------------------- 
    833828      !!                  *** ROUTINE interpvn *** 
    … …  
    837832      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 
    838833      LOGICAL, INTENT(in) :: before 
    839       INTEGER, INTENT(in) :: nb , ndir 
    840834      ! 
    841835      INTEGER :: ji,jj,jk 
    … …  
    844838      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
    845839      REAL(wp), DIMENSION(1:jpk) :: h_out 
    846       INTEGER  :: N_in, N_out, jref 
     840      INTEGER  :: N_in, N_out 
    847841      REAL(wp) :: h_diff 
    848       LOGICAL  :: northern_side,southern_side 
    849842      !!---------------------------------------------     
    850843      !       
    … …  
    864857# if defined key_vertical 
    865858 
    866          southern_side = (nb == 2).AND.(ndir == 1) 
    867          northern_side = (nb == 2).AND.(ndir == 2) 
    868  
    869859         DO jj=j1,j2 
    870             jref = jj 
    871             IF (southern_side) jref = MAX(2,jj) 
    872             IF (northern_side) jref = MIN(nlcj-2,jj) 
    873860            DO ji=i1,i2 
    874861               N_in = 0 
    … …  
    886873               N_out = 0 
    887874               DO jk=1,jpk 
    888                   if (vmask(ji,jref,jk) == 0) EXIT 
     875                  if (vmask(ji,jj,jk) == 0) EXIT 
    889876                  N_out = N_out + 1 
    890                   h_out(N_out) = e3v_a(ji,jref,jk) 
     877                  h_out(N_out) = e3v_a(ji,jj,jk) 
    891878               END DO 
    892879               IF (N_out == 0) THEN 

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90

    r11574 r11590  
    2323   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    2424   USE iom 
     25   USE vremap 
    2526 
    2627   IMPLICIT NONE 
    … …  
    319320            DO jj=j1,j2 
    320321               DO ji=i1,i2 
    321                   tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)  
     322                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_b(ji,jj,jk)  
    322323               END DO 
    323324            END DO 
    … …  
    342343                  IF (tmask(ji,jj,jk) == 0) EXIT  
    343344                  N_out = N_out + 1 
    344                   h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
     345                  h_out(jk) = e3t_b(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
    345346               ENDDO 
    346347               IF (N_in > 0) THEN 
    347                   h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    348                   tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for????? 
    349348                  DO jn=1,jpts 
    350                      call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
     349                     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) 
    351350                  ENDDO 
    352351               ENDIF 
    … …  
    446445                  tabres(ji,jj,jk,m1) = ub(ji,jj,jk) 
    447446# if defined key_vertical 
    448                   tabres(ji,jj,jk,m2) = e3u_n(ji,jj,jk)*umask(ji,jj,jk) 
     447                  tabres(ji,jj,jk,m2) = e3u_b(ji,jj,jk)*umask(ji,jj,jk) 
    449448# endif 
    450449               END DO 
    … …  
    475474                 if (umask(ji,jj,jk) == 0) EXIT 
    476475                 N_out = N_out + 1 
    477                  h_out(N_out) = e3u_n(ji,jj,jk) 
     476                 h_out(N_out) = e3u_b(ji,jj,jk) 
    478477              ENDDO 
    479478          
    … …  
    507506            DO jj = j1,j2 
    508507               DO ji = i1+1,i2   ! vector opt. 
    509                   zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj) 
    510                   hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u_n(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) & 
    511                                      &   -e2u(ji-1,jj)*e3u_n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr 
     508                  zbtr = r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) * fsahm_spt(ji,jj) 
     509                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u_b(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) & 
     510                                     &   -e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr 
    512511               END DO 
    513512            END DO 
    … …  
    599598                  tabres(ji,jj,jk,m1) = vb(ji,jj,jk) 
    600599# if defined key_vertical 
    601                   tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_n(ji,jj,jk) 
     600                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_b(ji,jj,jk) 
    602601# endif 
    603602               END DO 
    … …  
    627626                 if (vmask(ji,jj,jk) == 0) EXIT 
    628627                 N_out = N_out + 1 
    629                  h_out(N_out) = e3v_n(ji,jj,jk) 
     628                 h_out(N_out) = e3v_b(ji,jj,jk) 
    630629              ENDDO 
    631630          
    … …  
    653652            DO jj = j1+1,j2 
    654653               DO ji = i1,i2   ! vector opt. 
    655                   zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj) 
    656                   hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  & 
    657                                      &  -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr 
     654                  zbtr = r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) * fsahm_spt(ji,jj) 
     655                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v_b(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  & 
     656                                     &  -e1v(ji,jj-1) * e3v_b(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr 
    658657               END DO 
    659658            END DO 

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_update.F90

    r11574 r11590  
    2424   USE lib_mpp        ! MPP library 
    2525   USE domvvl         ! Need interpolation routines  
     26   USE vremap         ! Vertical remapping 
    2627 
    2728   IMPLICIT NONE 
    … …  
    319320               N_in = 0 
    320321               DO jk=k1,k2 !k2 = jpk of child grid 
    321                   IF (tabres(ji,jj,jk,n2) < 900  ) EXIT 
     322                  IF (tabres(ji,jj,jk,n2) < -900._wp  ) EXIT 
    322323                  N_in = N_in + 1 
    323324                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_interp.F90

    r11574 r11590  
    4848   END SUBROUTINE Agrif_trc 
    4949 
    50    SUBROUTINE interptrn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
     50   SUBROUTINE interptrn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    5151      !!---------------------------------------------------------------------- 
    5252      !!                  *** ROUTINE interptrn *** 
    … …  
    5555      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2 
    5656      LOGICAL                                     , INTENT(in   ) ::   before 
    57       INTEGER                                     , INTENT(in   ) ::   nb , ndir 
    5857      ! 
    59       INTEGER  ::   ji, jj, jk, jn, iref, jref, ibdy, jbdy   ! dummy loop indices 
     58      INTEGER  ::   ji, jj, jk, jn, ibdy, jbdy   ! dummy loop indices 
    6059      INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out 
    6160      REAL(wp) ::   zrho, z1, z2, z3, z4, z5, z6, z7 
    62       LOGICAL :: western_side, eastern_side,northern_side,southern_side 
     61 
    6362      ! vertical interpolation: 
    6463      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 
    … …  
    6665      REAL(wp), DIMENSION(k1:k2) :: h_in 
    6766      REAL(wp), DIMENSION(1:jpk) :: h_out 
    68       REAL(wp) :: h_diff 
     67      !!---------------------------------------------------------------------- 
    6968 
    7069      IF( before ) THEN          
    … …  
    9190 
    9291# if defined key_vertical 
    93          western_side  = (nb == 1).AND.(ndir == 1)   ;   eastern_side  = (nb == 1).AND.(ndir == 2) 
    94          southern_side = (nb == 2).AND.(ndir == 1)   ;   northern_side = (nb == 2).AND.(ndir == 2) 
    95  
    9692         DO jj=j1,j2 
    9793            DO ji=i1,i2 
    98                iref = ji 
    99                jref = jj 
    100                if(western_side) iref=MAX(2,ji) 
    101                if(eastern_side) iref=MIN(nlci-1,ji) 
    102                if(southern_side) jref=MAX(2,jj) 
    103                if(northern_side) jref=MIN(nlcj-1,jj) 
    10494               N_in = 0 
    10595               DO jk=k1,k2 !k2 = jpk of parent grid 
    … …  
    111101               N_out = 0 
    112102               DO jk=1,jpk ! jpk of child grid 
    113                   IF (tmask(iref,jref,jk) == 0) EXIT  
     103                  IF (tmask(ji,jj,jk) == 0) EXIT  
    114104                  N_out = N_out + 1 
    115                   h_out(jk) = e3t_n(iref,jref,jk) 
     105                  h_out(jk) = e3t_a(ji,jj,jk) 
    116106               ENDDO 
    117107               IF (N_in > 0) THEN 

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_top_sponge.F90

    r10068 r11590  
    9393            DO jj=j1,j2 
    9494               DO ji=i1,i2 
    95                   tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)  
     95                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_b(ji,jj,jk)  
    9696               END DO 
    9797            END DO 
    … …  
    114114                  IF (tmask(ji,jj,jk) == 0) EXIT  
    115115                  N_out = N_out + 1 
    116                   h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
     116                  h_out(jk) = e3t_b(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
    117117               ENDDO 
    118118               IF (N_in > 0) THEN 
    119                   h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    120                   tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for????? 
    121119                  DO jn=1,jptra 
    122120                     call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 

  • NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_user.F90

    r11574 r11590  
    324324      ! 3. Location of interpolation 
    325325      !----------------------------- 
    326       CALL Agrif_Set_bc(       tsn_id, (/0,ind1/) ) 
    327       CALL Agrif_Set_bc( un_interp_id, (/0,ind1/) ) 
    328       CALL Agrif_Set_bc( vn_interp_id, (/0,ind1/) ) 
     326      CALL Agrif_Set_bc(       tsn_id, (/0,ind1-1/) ) 
     327      CALL Agrif_Set_bc( un_interp_id, (/0,ind1-1/) ) 
     328      CALL Agrif_Set_bc( vn_interp_id, (/0,ind1-1/) ) 
    329329 
    330330      CALL Agrif_Set_bc( tsn_sponge_id, (/-nn_sponge_len*Agrif_irhox()-1,0/) )  ! if west and rhox=3 and sponge=2 and ghost=1: columns 2 to 9  
    … …  
    594594      ! 3. Location of interpolation 
    595595      !----------------------------- 
    596       CALL Agrif_Set_bc(trn_id,(/0,ind1/)) 
     596      CALL Agrif_Set_bc(trn_id,(/0,ind1-1/)) 
    597597      CALL Agrif_Set_bc(trn_sponge_id,(/-nn_sponge_len*Agrif_irhox()-1,0/)) 
    598598 
Note: See TracChangeset for help on using the changeset viewer.