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 9031 for branches/2017 – NEMO

Changeset 9031 for branches/2017


Ignore:
Timestamp:
2017-12-14T11:10:02+01:00 (6 years ago)
Author:
timgraham
Message:

Resolved AGRIF conflicts

Location:
branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_lim3_update.F90

    r9019 r9031  
    5959      Agrif_UseSpecialValueInUpdate = .TRUE. 
    6060# if defined TWO_WAY 
    61       IF( MOD(nbcline,nbclineupdate) == 0) THEN ! update the whole basin at each nbclineupdate (=nn_cln_update) baroclinic parent time steps 
    62                                                 ! nbcline is incremented (+1) at the end of each parent time step from 0 (1st time step) 
    63          CALL Agrif_Update_Variable( tra_ice_id , procname = update_tra_ice  ) 
    64          CALL Agrif_Update_Variable( u_ice_id   , procname = update_u_ice    ) 
    65          CALL Agrif_Update_Variable( v_ice_id   , procname = update_v_ice    ) 
    66       ELSE                                      ! update only the boundaries defined par locupdate 
    67          CALL Agrif_Update_Variable( tra_ice_id , locupdate=(/0,2/), procname = update_tra_ice  ) 
    68          CALL Agrif_Update_Variable( u_ice_id   , locupdate=(/0,1/), procname = update_u_ice    ) 
    69          CALL Agrif_Update_Variable( v_ice_id   , locupdate=(/0,1/), procname = update_v_ice    ) 
    70       ENDIF 
     61      CALL Agrif_Update_Variable( tra_ice_id , procname = update_tra_ice  ) 
     62      CALL Agrif_Update_Variable( u_ice_id   , procname = update_u_ice    ) 
     63      CALL Agrif_Update_Variable( v_ice_id   , procname = update_v_ice    ) 
     64 
     65!      CALL Agrif_Update_Variable( tra_ice_id , locupdate=(/0,2/), procname = update_tra_ice  ) 
     66!      CALL Agrif_Update_Variable( u_ice_id   , locupdate=(/0,1/), procname = update_u_ice    ) 
     67!      CALL Agrif_Update_Variable( v_ice_id   , locupdate=(/0,1/), procname = update_v_ice    ) 
    7168# endif 
    7269      Agrif_SpecialValueFineGrid    = 0. 

  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90

    r9019 r9031  
    1717 
    1818   PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 
    19  
     19#if defined key_vertical 
     20   PUBLIC reconstructandremap ! remapping routine 
     21#endif    
    2022   !                                              !!* Namelist namagrif: AGRIF parameters 
    2123   LOGICAL , PUBLIC ::   ln_spc_dyn    = .FALSE.   !: 
    22    INTEGER , PUBLIC ::   nn_cln_update = 3         !: update frequency  
    2324   INTEGER , PUBLIC, PARAMETER ::   nn_sponge_len = 2  !: Sponge width (in number of parent grid points) 
    2425   REAL(wp), PUBLIC ::   rn_sponge_tra = 2800.     !: sponge coeff. for tracers 
    2526   REAL(wp), PUBLIC ::   rn_sponge_dyn = 2800.     !: sponge coeff. for dynamics 
    2627   LOGICAL , PUBLIC ::   ln_chk_bathy  = .FALSE.   !: check of parent bathymetry  
    27  
     28   LOGICAL , PUBLIC ::   lk_agrif_clp  = .TRUE.    !: Force clamped bcs 
    2829   !                                              !!! OLD namelist names 
    29    INTEGER , PUBLIC ::   nbcline = 0               !: update counter 
    30    INTEGER , PUBLIC ::   nbclineupdate             !: update frequency  
    3130   REAL(wp), PUBLIC ::   visc_tra                  !: sponge coeff. for tracers 
    3231   REAL(wp), PUBLIC ::   visc_dyn                  !: sponge coeff. for dynamics 
    … …  
    3534   LOGICAL , PUBLIC :: spongedoneU = .FALSE.       !: dynamics sponge layer indicator 
    3635   LOGICAL , PUBLIC :: lk_agrif_fstep = .TRUE.     !: if true: first step 
    37    LOGICAL , PUBLIC :: lk_agrif_doupd = .TRUE.     !: if true: send update from current grid 
    3836   LOGICAL , PUBLIC :: lk_agrif_debug = .FALSE.    !: if true: print debugging info 
    3937 
    … …  
    102100   END FUNCTION agrif_oce_alloc 
    103101 
     102#if defined key_vertical 
     103   SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)       
     104      !!---------------------------------------------------------------------- 
     105      !!                ***  FUNCTION reconstructandremap  *** 
     106      !!---------------------------------------------------------------------- 
     107      IMPLICIT NONE 
     108      INTEGER N, Nout 
     109      REAL(wp) tabin(N), tabout(Nout) 
     110      REAL(wp) hin(N), hout(Nout) 
     111      REAL(wp) coeffremap(N,3),zwork(N,3) 
     112      REAL(wp) zwork2(N+1,3) 
     113      INTEGER jk 
     114      DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8   
     115      REAL(wp) q,q01,q02,q001,q002,q0 
     116      REAL(wp) z_win(1:N+1), z_wout(1:Nout+1) 
     117      REAL(wp),PARAMETER :: dpthin = 1.D-3 
     118      INTEGER :: k1, kbox, ktop, ka, kbot 
     119      REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 
     120 
     121      z_win(1)=0.; z_wout(1)= 0. 
     122      DO jk=1,N 
     123         z_win(jk+1)=z_win(jk)+hin(jk) 
     124      ENDDO  
     125       
     126      DO jk=1,Nout 
     127         z_wout(jk+1)=z_wout(jk)+hout(jk)        
     128      ENDDO        
     129 
     130      DO jk=2,N 
     131         zwork(jk,1)=1./(hin(jk-1)+hin(jk)) 
     132      ENDDO 
     133         
     134      DO jk=2,N-1 
     135         q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1)) 
     136         zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1) 
     137         zwork(jk,3)=q0 
     138      ENDDO        
     139      
     140      DO jk= 2,N 
     141         zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1)) 
     142      ENDDO 
     143         
     144      coeffremap(:,1) = tabin(:) 
     145  
     146      DO jk=2,N-1 
     147         q001 = hin(jk)*zwork2(jk+1,1) 
     148         q002 = hin(jk)*zwork2(jk,1)         
     149         IF (q001*q002 < 0) then 
     150            q001 = 0. 
     151            q002 = 0. 
     152         ENDIF 
     153         q=zwork(jk,2) 
     154         q01=q*zwork2(jk+1,1) 
     155         q02=q*zwork2(jk,1) 
     156         IF (abs(q001) > abs(q02)) q001 = q02 
     157         IF (abs(q002) > abs(q01)) q002 = q01 
     158         
     159         q=(q001-q002)*zwork(jk,3) 
     160         q001=q001-q*hin(jk+1) 
     161         q002=q002+q*hin(jk-1) 
     162         
     163         coeffremap(jk,3)=coeffremap(jk,1)+q001 
     164         coeffremap(jk,2)=coeffremap(jk,1)-q002 
     165         
     166         zwork2(jk,1)=(2.*q001-q002)**2 
     167         zwork2(jk,2)=(2.*q002-q001)**2 
     168      ENDDO 
     169         
     170      DO jk=1,N 
     171         IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN 
     172            coeffremap(jk,3) = coeffremap(jk,1) 
     173            coeffremap(jk,2) = coeffremap(jk,1) 
     174            zwork2(jk,1) = 0. 
     175            zwork2(jk,2) = 0. 
     176         ENDIF 
     177      ENDDO 
     178         
     179      DO jk=2,N 
     180         q002=max(zwork2(jk-1,2),dsmll) 
     181         q001=max(zwork2(jk,1),dsmll) 
     182         zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002) 
     183      ENDDO 
     184         
     185      zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) 
     186      zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) 
     187  
     188      DO jk=1,N 
     189         q01=zwork2(jk+1,3)-coeffremap(jk,1) 
     190         q02=coeffremap(jk,1)-zwork2(jk,3) 
     191         q001=2.*q01 
     192         q002=2.*q02 
     193         IF (q01*q02<0) then 
     194            q01=0. 
     195            q02=0. 
     196         ELSEIF (abs(q01)>abs(q002)) then 
     197            q01=q002 
     198         ELSEIF (abs(q02)>abs(q001)) then 
     199            q02=q001 
     200         ENDIF 
     201         coeffremap(jk,2)=coeffremap(jk,1)-q02 
     202         coeffremap(jk,3)=coeffremap(jk,1)+q01 
     203      ENDDO 
     204 
     205      zbot=0.0 
     206      kbot=1 
     207      DO jk=1,Nout 
     208         ztop=zbot  !top is bottom of previous layer 
     209         ktop=kbot 
     210         IF     (ztop.GE.z_win(ktop+1)) then 
     211            ktop=ktop+1 
     212         ENDIF 
     213         
     214         zbot=z_wout(jk+1) 
     215         zthk=zbot-ztop 
     216 
     217         IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN 
     218 
     219            kbot=ktop 
     220            DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 
     221               kbot=kbot+1 
     222            ENDDO 
     223            zbox=zbot 
     224            DO k1= jk+1,Nout 
     225               IF     (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN 
     226                  exit !thick layer 
     227               ELSE 
     228                  zbox=z_wout(k1+1)  !include thin adjacent layers 
     229                  IF(zbox.EQ.z_wout(Nout+1)) THEN 
     230                     exit !at bottom 
     231                  ENDIF 
     232               ENDIF 
     233            ENDDO 
     234            zthk=zbox-ztop 
     235 
     236            kbox=ktop 
     237            DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 
     238               kbox=kbox+1 
     239            ENDDO 
     240           
     241            IF(ktop.EQ.kbox) THEN 
     242               IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN 
     243                  IF(hin(kbox).GT.dpthin) THEN 
     244                     q001 = (zbox-z_win(kbox))/hin(kbox) 
     245                     q002 = (ztop-z_win(kbox))/hin(kbox) 
     246                     q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) 
     247                     q02=q01-1.+(q001+q002) 
     248                     q0=1.-q01-q02 
     249                  ELSE 
     250                     q0 = 1.0 
     251                     q01 = 0. 
     252                     q02 = 0. 
     253                  ENDIF 
     254                  tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 
     255               ELSE 
     256                  tabout(jk) = tabin(kbox) 
     257               ENDIF  
     258            ELSE 
     259               IF(ktop.LE.jk .AND. kbox.GE.jk) THEN 
     260                  ka = jk 
     261               ELSEIF (kbox-ktop.GE.3) THEN 
     262                  ka = (kbox+ktop)/2 
     263               ELSEIF (hin(ktop).GE.hin(kbox)) THEN 
     264                  ka = ktop 
     265               ELSE 
     266                  ka = kbox 
     267               ENDIF !choose ka 
     268     
     269               offset=coeffremap(ka,1) 
     270     
     271               qtop = z_win(ktop+1)-ztop !partial layer thickness 
     272               IF(hin(ktop).GT.dpthin) THEN 
     273                  q=(ztop-z_win(ktop))/hin(ktop) 
     274                  q01=q*(q-1.) 
     275                  q02=q01+q 
     276                  q0=1-q01-q02             
     277               ELSE 
     278                  q0 = 1. 
     279                  q01 = 0. 
     280                  q02 = 0. 
     281               ENDIF 
     282                
     283               tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop 
     284     
     285               DO k1= ktop+1,kbox-1 
     286                  tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) 
     287               ENDDO !k1 
     288                
     289               qbot = zbox-z_win(kbox) !partial layer thickness 
     290               IF(hin(kbox).GT.dpthin) THEN 
     291                  q=qbot/hin(kbox) 
     292                  q01=(q-1.)**2 
     293                  q02=q01-1.+q 
     294                  q0=1-q01-q02                             
     295               ELSE 
     296                  q0 = 1.0 
     297                  q01 = 0. 
     298                  q02 = 0. 
     299               ENDIF 
     300               
     301               tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot 
     302                
     303               rpsum=1.0d0/zthk 
     304               tabout(jk)=offset+tsum*rpsum 
     305                  
     306            ENDIF !single or multiple layers 
     307         ELSE 
     308            IF (jk==1) THEN 
     309               write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1) 
     310            ENDIF 
     311            tabout(jk) = tabout(jk-1) 
     312              
     313         ENDIF !normal:thin layer 
     314      ENDDO !jk 
     315             
     316      return 
     317      end subroutine reconstructandremap 
     318#endif 
     319 
    104320#endif 
    105321   !!====================================================================== 

  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90

    r9019 r9031  
    2828   USE agrif_oce 
    2929   USE phycst 
     30   USE dynspg_ts, ONLY: un_adv, vn_adv 
    3031   ! 
    3132   USE in_out_manager 
    … …  
    4243   PUBLIC   interpe3t, interpumsk, interpvmsk 
    4344   PUBLIC   Agrif_avm, interpavm 
     45>>>>>>> .merge-right.r9019 
    4446 
    4547   INTEGER ::   bdy_tinterp = 0 
    … …  
    196198               END DO 
    197199            END DO 
    198              
    199             zub(nlci-2,:) = 0._wp        ! Correct transport 
     200         ENDIF 
     201         zub(nlci-2,:) = 0._wp        ! Correct transport 
     202         DO jk = 1, jpkm1 
     203            DO jj = 1, jpj 
     204               zub(nlci-2,jj) = zub(nlci-2,jj) + e3u_a(nlci-2,jj,jk) * ua(nlci-2,jj,jk) 
     205            END DO 
     206         END DO 
     207         DO jj = 1, jpj 
     208            zub(nlci-2,jj) = zub(nlci-2,jj) * r1_hu_a(nlci-2,jj) 
     209         END DO 
     210 
     211         DO jk = 1, jpkm1 
     212            DO jj = 1, jpj 
     213               ua(nlci-2,jj,jk) = ( ua(nlci-2,jj,jk) + ua_b(nlci-2,jj) - zub(nlci-2,jj) ) * umask(nlci-2,jj,jk) 
     214            END DO 
     215         END DO 
     216         ! 
     217         ! Set tangential velocities to time splitting estimate 
     218         !----------------------------------------------------- 
     219         IF( ln_dynspg_ts ) THEN 
     220            zvb(nlci-1,:) = 0._wp 
    200221            DO jk = 1, jpkm1 
    201222               DO jj = 1, jpj 
    … …  
    258279         ENDIF 
    259280         ! 
    260          ! Smoothing if only 1 ghostcell 
    261          ! ----------------------------- 
    262          IF( nbghostcells == 1 ) THEN 
     281         IF (.NOT.lk_agrif_clp) THEN 
    263282            DO jk = 1, jpkm1              ! Smooth 
    264283               DO ji = i1, i2 
    … …  
    267286               END DO 
    268287            END DO 
    269             ! 
    270             zvb(:,2) = 0._wp              ! Correct transport 
    271             DO jk=1,jpkm1 
    272                DO ji=1,jpi 
    273                   zvb(ji,2) = zvb(ji,2) + e3v_a(ji,2,jk) * va(ji,2,jk) * vmask(ji,2,jk) 
    274                END DO 
     288         ENDIF 
     289         ! 
     290         zvb(:,2) = 0._wp              ! Correct transport 
     291         DO jk=1,jpkm1 
     292            DO ji=1,jpi 
     293               zvb(ji,2) = zvb(ji,2) + e3v_a(ji,2,jk) * va(ji,2,jk) * vmask(ji,2,jk) 
    275294            END DO 
    276295            DO ji = 1, jpi 
    … …  
    330349         ENDIF 
    331350         ! 
    332          ! Smoothing if only 1 ghostcell 
    333          ! ----------------------------- 
    334          IF( nbghostcells == 1 ) THEN 
     351         IF (.NOT.lk_agrif_clp) THEN 
    335352            DO jk = 1, jpkm1              ! Smooth 
    336353               DO ji = i1, i2 
    … …  
    339356               END DO 
    340357            END DO 
    341             ! 
    342             zvb(:,nlcj-2) = 0._wp         ! Correct transport 
    343             DO jk = 1, jpkm1 
    344                DO ji = 1, jpi 
    345                   zvb(ji,nlcj-2) = zvb(ji,nlcj-2) + e3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk) 
    346                END DO 
     358         END IF 
     359         ! 
     360         zvb(:,nlcj-2) = 0._wp         ! Correct transport 
     361         DO jk = 1, jpkm1 
     362            DO ji = 1, jpi 
     363               zvb(ji,nlcj-2) = zvb(ji,nlcj-2) + e3v_a(ji,nlcj-2,jk) * va(ji,nlcj-2,jk) * vmask(ji,nlcj-2,jk) 
    347364            END DO 
    348365            DO ji = 1, jpi 
    … …  
    455472      INTEGER :: ji, jj 
    456473      LOGICAL :: ll_int_cons 
    457       REAL(wp) :: zrhot, zt 
    458474      !!----------------------------------------------------------------------   
    459475      ! 
    … …  
    462478      ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only 
    463479      ! 
    464       zrhot = Agrif_rhot() 
    465       ! 
    466       ! "Central" time index for interpolation: 
    467       IF( ln_bt_fw ) THEN 
    468          zt = REAL( Agrif_NbStepint()+0.5_wp, wp ) / zrhot 
    469       ELSE 
    470          zt = REAL( Agrif_NbStepint()       , wp ) / zrhot 
    471       ENDIF 
    472       ! 
    473       ! Linear interpolation of sea level 
    474       Agrif_SpecialValue    = 0._wp 
    475       Agrif_UseSpecialValue = .TRUE. 
    476       CALL Agrif_Bc_variable( sshn_id, calledweight=zt, procname=interpsshn ) 
    477       Agrif_UseSpecialValue = .FALSE. 
     480      ! Enforce volume conservation if no time refinement:   
     481      IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE.   
    478482      ! 
    479483      ! Interpolate barotropic fluxes 
    480       Agrif_SpecialValue=0. 
     484      Agrif_SpecialValue=0._wp 
    481485      Agrif_UseSpecialValue = ln_spc_dyn 
    482486      ! 
    … …  
    497501         ubdy_n(:) = 0._wp   ;   vbdy_n(:) = 0._wp  
    498502         ubdy_s(:) = 0._wp   ;   vbdy_s(:) = 0._wp 
    499          CALL Agrif_Bc_variable( unb_id, calledweight=zt, procname=interpunb ) 
    500          CALL Agrif_Bc_variable( vnb_id, calledweight=zt, procname=interpvnb ) 
     503         CALL Agrif_Bc_variable( unb_id, procname=interpunb ) 
     504         CALL Agrif_Bc_variable( vnb_id, procname=interpvnb ) 
    501505      ENDIF 
    502506      Agrif_UseSpecialValue = .FALSE. 
    … …  
    507511   SUBROUTINE Agrif_ssh( kt ) 
    508512      !!---------------------------------------------------------------------- 
    509       !!                  ***  ROUTINE Agrif_DYN  *** 
     513      !!                  ***  ROUTINE Agrif_ssh  *** 
    510514      !!----------------------------------------------------------------------   
    511515      INTEGER, INTENT(in) ::   kt 
    512516      ! 
    513517      INTEGER  :: ji, jj, indx 
     518      INTEGER :: ji, jj 
    514519      !!----------------------------------------------------------------------   
    515520      ! 
    516521      IF( Agrif_Root() )   RETURN 
    517       !! clem ghost 
    518       ! --- West --- ! 
     522      !       
     523      ! Linear interpolation in time of sea level 
     524      ! 
     525      Agrif_SpecialValue    = 0._wp 
     526      Agrif_UseSpecialValue = .TRUE. 
     527      CALL Agrif_Bc_variable(sshn_id, procname=interpsshn ) 
     528      Agrif_UseSpecialValue = .FALSE. 
     529      ! 
    519530      IF((nbondi == -1).OR.(nbondi == 2)) THEN 
    520531         indx = 1+nbghostcells 
    521532         DO jj = 1, jpj 
    522533            DO ji = 2, indx 
    523                ssha(ji,jj)=ssha(indx+1,jj) 
    524                sshn(ji,jj)=sshn(indx+1,jj) 
     534               ssha(ji,jj) = hbdy_w(jj) 
    525535            ENDDO 
    526536         ENDDO 
    … …  
    532542         DO jj = 1, jpj 
    533543            DO ji = indx, nlci-1 
    534                ssha(ji,jj)=ssha(indx-1,jj) 
    535                sshn(ji,jj)=sshn(indx-1,jj) 
     544               ssha(indx,jj) = hbdy_e(jj) 
    536545            ENDDO 
    537546         ENDDO 
    … …  
    543552         DO jj = 2, indx 
    544553            DO ji = 1, jpi 
    545                ssha(ji,jj)=ssha(ji,indx+1) 
    546                sshn(ji,jj)=sshn(ji,indx+1) 
     554               ssha(ji,indx) = hbdy_s(ji) 
    547555            ENDDO 
    548556         ENDDO 
    … …  
    554562         DO jj = indx, nlcj-1 
    555563            DO ji = 1, jpi 
    556                ssha(ji,jj)=ssha(ji,indx-1) 
    557                sshn(ji,jj)=sshn(ji,indx-1) 
     564               ssha(ji,indx) = hbdy_n(ji) 
    558565            ENDDO 
    559566         ENDDO 
    … …  
    572579      !!----------------------------------------------------------------------   
    573580      !! clem ghost (starting at i,j=1 is important I think otherwise you introduce a grad(ssh)/=0 at point 2) 
     581      ! 
     582      IF( Agrif_Root() )   RETURN 
     583      ! 
    574584      IF((nbondi == -1).OR.(nbondi == 2)) THEN 
    575585         DO jj = 1, jpj 
    … …  
    598608   END SUBROUTINE Agrif_ssh_ts 
    599609 
    600  
    601610   SUBROUTINE Agrif_avm 
    602611      !!---------------------------------------------------------------------- 
    … …  
    606615      !!----------------------------------------------------------------------   
    607616      ! 
    608       zalpha = 1._wp   ! proper time interpolation impossible  ==> use last available value from parent  
    609       ! 
    610       Agrif_SpecialValue    = 0._wp 
     617      IF( Agrif_Root() )   RETURN 
     618      ! 
     619      zalpha = 1._wp ! JC: proper time interpolation impossible   
     620                     ! => use last available value from parent  
     621      ! 
     622      Agrif_SpecialValue    = 0.e0 
    611623      Agrif_UseSpecialValue = .TRUE. 
    612624      ! 
    … …  
    627639      INTEGER                                     , INTENT(in   ) ::   nb , ndir 
    628640      ! 
    629       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    630       INTEGER  ::   imin, imax, jmin, jmax 
     641      INTEGER  ::   ji, jj, jk, jn, iref, jref   ! dummy loop indices 
     642      INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out 
    631643      REAL(wp) ::   zrhox, z1, z2, z3, z4, z5, z6, z7 
    632       LOGICAL  ::   western_side, eastern_side,northern_side,southern_side 
    633       !!---------------------------------------------------------------------- 
    634       ! 
     644      LOGICAL :: western_side, eastern_side,northern_side,southern_side 
     645      ! vertical interpolation: 
     646      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 
     647      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 
     648      REAL(wp), DIMENSION(k1:k2) :: h_in 
     649      REAL(wp), DIMENSION(1:jpk) :: h_out(1:jpk) 
     650      REAL(wp) :: h_diff, zrhoxy 
     651 
     652      zrhoxy = Agrif_rhox()*Agrif_rhoy() 
    635653      IF( before ) THEN          
    636          ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 
    637       ELSE 
    638          ! 
     654         DO jn = 1,jpts 
     655            DO jk=k1,k2 
     656               DO jj=j1,j2 
     657                 DO ji=i1,i2 
     658                       ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 
     659                 END DO 
     660              END DO 
     661           END DO 
     662        END DO 
     663 
     664# if defined key_vertical 
     665        DO jk=k1,k2 
     666           DO jj=j1,j2 
     667              DO ji=i1,i2 
     668                 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)  
     669              END DO 
     670           END DO 
     671        END DO 
     672# endif 
     673      ELSE  
     674 
    639675         western_side  = (nb == 1).AND.(ndir == 1)   ;   eastern_side  = (nb == 1).AND.(ndir == 2) 
    640676         southern_side = (nb == 2).AND.(ndir == 1)   ;   northern_side = (nb == 2).AND.(ndir == 2) 
     677 
     678# if defined key_vertical               
     679         DO jj=j1,j2 
     680            DO ji=i1,i2 
     681               iref = ji 
     682               jref = jj 
     683               if(western_side) iref=MAX(2,ji) 
     684               if(eastern_side) iref=MIN(nlci-1,ji) 
     685               if(southern_side) jref=MAX(2,jj) 
     686               if(northern_side) jref=MIN(nlcj-1,jj) 
     687               N_in = 0 
     688               DO jk=k1,k2 !k2 = jpk of parent grid 
     689                  IF (ptab(ji,jj,jk,n2) == 0) EXIT 
     690                  N_in = N_in + 1 
     691                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 
     692                  h_in(N_in) = ptab(ji,jj,jk,n2) 
     693               END DO 
     694               N_out = 0 
     695               DO jk=1,jpk ! jpk of child grid 
     696                  IF (tmask(iref,jref,jk) == 0) EXIT  
     697                  N_out = N_out + 1 
     698                  h_out(jk) = e3t_n(iref,jref,jk) 
     699               ENDDO 
     700               IF (N_in > 0) THEN 
     701                  DO jn=1,jpts 
     702                     call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
     703                  ENDDO 
     704               ENDIF 
     705            ENDDO 
     706         ENDDO 
     707# else 
     708         ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts) 
     709# endif 
    641710         ! 
    642711         IF( nbghostcells > 1 ) THEN  ! no smoothing 
    643             tsa(i1:i2,j1:j2,k1:k2,n1:n2) = ptab(i1:i2,j1:j2,k1:k2,n1:n2) 
     712            tsa(i1:i2,j1:j2,k1:k2,n1:n2) = ptab_child(i1:i2,j1:j2,k1:k2,n1:n2) 
    644713         ELSE                         ! smoothing 
    645714            ! 
    … …  
    665734            IF( eastern_side ) THEN 
    666735               DO jn = 1, jpts 
    667                   tsa(nlci,j1:j2,k1:k2,jn) = z1 * ptab(nlci,j1:j2,k1:k2,jn) + z2 * ptab(nlci-1,j1:j2,k1:k2,jn) 
     736                  tsa(nlci,j1:j2,k1:k2,jn) = z1 * ptab_child(nlci,j1:j2,k1:k2,jn) + z2 * ptab_child(nlci-1,j1:j2,k1:k2,jn) 
    668737                  DO jk = 1, jpkm1 
    669738                     DO jj = jmin,jmax 
    … …  
    685754            IF( northern_side ) THEN             
    686755               DO jn = 1, jpts 
    687                   tsa(i1:i2,nlcj,k1:k2,jn) = z1 * ptab(i1:i2,nlcj,k1:k2,jn) + z2 * ptab(i1:i2,nlcj-1,k1:k2,jn) 
     756                  tsa(i1:i2,nlcj,k1:k2,jn) = z1 * ptab_child(i1:i2,nlcj,k1:k2,jn) + z2 * ptab_child(i1:i2,nlcj-1,k1:k2,jn) 
    688757                  DO jk = 1, jpkm1 
    689758                     DO ji = imin,imax 
    … …  
    705774            IF( western_side ) THEN             
    706775               DO jn = 1, jpts 
    707                   tsa(1,j1:j2,k1:k2,jn) = z1 * ptab(1,j1:j2,k1:k2,jn) + z2 * ptab(2,j1:j2,k1:k2,jn) 
     776                  tsa(1,j1:j2,k1:k2,jn) = z1 * ptab_child(1,j1:j2,k1:k2,jn) + z2 * ptab_child(2,j1:j2,k1:k2,jn) 
    708777                  DO jk = 1, jpkm1 
    709778                     DO jj = jmin,jmax 
    … …  
    724793            IF( southern_side ) THEN            
    725794               DO jn = 1, jpts 
    726                   tsa(i1:i2,1,k1:k2,jn) = z1 * ptab(i1:i2,1,k1:k2,jn) + z2 * ptab(i1:i2,2,k1:k2,jn) 
     795                  tsa(i1:i2,1,k1:k2,jn) = z1 * ptab_child(i1:i2,1,k1:k2,jn) + z2 * ptab_child(i1:i2,2,k1:k2,jn) 
    727796                  DO jk = 1, jpk       
    728797                     DO ji=imin,imax 
    … …  
    741810            ENDIF 
    742811            ! 
     812            ! 
    743813            ! Treatment of corners 
    744             IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2)))   tsa(nlci-1,2,:,:) = ptab(nlci-1,2,:,:)            ! East south 
    745             IF ((eastern_side).AND.((nbondj ==  1).OR.(nbondj == 2)))   tsa(nlci-1,nlcj-1,:,:) = ptab(nlci-1,nlcj-1,:,:)  ! East north 
    746             IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2)))   tsa(2,2,:,:) = ptab(2,2,:,:)                      ! West south 
    747             IF ((western_side).AND.((nbondj ==  1).OR.(nbondj == 2)))   tsa(2,nlcj-1,:,:) = ptab(2,nlcj-1,:,:)            ! West north 
     814            !  
     815            ! East south 
     816            IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 
     817               tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,1:jpts) 
     818            ENDIF 
     819            ! East north 
     820            IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 
     821               tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,1:jpts) 
     822            ENDIF 
     823            ! West south 
     824            IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 
     825               tsa(2,2,:,:) = ptab_child(2,2,:,1:jpts) 
     826            ENDIF 
     827            ! West north 
     828            IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 
     829               tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,1:jpts) 
     830            ENDIF 
    748831            ! 
    749832         ENDIF 
    … …  
    751834      ! 
    752835   END SUBROUTINE interptsn 
    753  
    754836 
    755837   SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir ) 
    … …  
    781863   END SUBROUTINE interpsshn 
    782864 
    783  
    784    SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, before ) 
     865   SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir ) 
    785866      !!---------------------------------------------------------------------- 
    786867      !!                  *** ROUTINE interpun *** 
    787       !!---------------------------------------------------------------------- 
    788       INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2 
    789       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab 
    790       LOGICAL                               , INTENT(in   ) ::   before 
    791       ! 
    792       INTEGER  ::   ji, jj, jk 
    793       REAL(wp) ::   zrhoy   
    794       !!---------------------------------------------------------------------- 
    795       ! 
    796       IF( before ) THEN  
    797          DO jk = k1, jpk 
    798             ptab(i1:i2,j1:j2,jk) = e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 
     868      !!---------------------------------------------     
     869      !! 
     870      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 
     871      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 
     872      LOGICAL, INTENT(in) :: before 
     873      INTEGER, INTENT(in) :: nb , ndir 
     874      !! 
     875      INTEGER :: ji,jj,jk 
     876      REAL(wp) :: zrhoy 
     877      ! vertical interpolation: 
     878      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
     879      REAL(wp), DIMENSION(1:jpk) :: h_out 
     880      INTEGER  :: N_in, N_out, iref 
     881      REAL(wp) :: h_diff 
     882      LOGICAL  :: western_side, eastern_side 
     883      !!---------------------------------------------     
     884      ! 
     885      zrhoy = Agrif_rhoy() 
     886      IF (before) THEN  
     887         DO jk=1,jpk 
     888            DO jj=j1,j2 
     889               DO ji=i1,i2 
     890                  ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk))  
     891# if defined key_vertical 
     892                  ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)) 
     893# endif 
     894               END DO 
     895            END DO 
    799896         END DO 
    800897      ELSE 
    801          zrhoy = Agrif_Rhoy() 
     898         zrhoy = Agrif_rhoy() 
     899# if defined key_vertical 
     900! VERTICAL REFINEMENT BEGIN 
     901         western_side  = (nb == 1).AND.(ndir == 1) 
     902         eastern_side  = (nb == 1).AND.(ndir == 2) 
     903 
     904         DO ji=i1,i2 
     905            iref = ji 
     906            IF (western_side) iref = MAX(2,ji) 
     907            IF (eastern_side) iref = MIN(nlci-2,ji) 
     908            DO jj=j1,j2 
     909               N_in = 0 
     910               DO jk=k1,k2 
     911                  IF (ptab(ji,jj,jk,2) == 0) EXIT 
     912                  N_in = N_in + 1 
     913                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 
     914                  h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)  
     915              ENDDO 
     916          
     917              IF (N_in == 0) THEN 
     918                 ua(ji,jj,:) = 0._wp 
     919                 CYCLE 
     920              ENDIF 
     921          
     922              N_out = 0 
     923              DO jk=1,jpk 
     924                 if (umask(iref,jj,jk) == 0) EXIT 
     925                 N_out = N_out + 1 
     926                 h_out(N_out) = e3u_a(iref,jj,jk) 
     927              ENDDO 
     928          
     929              IF (N_out == 0) THEN 
     930                 ua(ji,jj,:) = 0._wp 
     931                 CYCLE 
     932              ENDIF 
     933          
     934              IF (N_in * N_out > 0) THEN 
     935                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     936! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 
     937                 if (h_diff < -1.e4) then 
     938                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 
     939!                    stop 
     940                 endif 
     941              ENDIF 
     942              call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     943            ENDDO 
     944         ENDDO 
     945 
     946# else 
    802947         DO jk = 1, jpkm1 
    803948            DO jj=j1,j2 
    804                ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_n(i1:i2,jj,jk) ) 
    805             END DO 
    806          END DO 
     949               ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) ) 
     950            END DO 
     951         END DO 
     952# endif 
     953 
    807954      ENDIF 
    808955      !  
    809956   END SUBROUTINE interpun 
    810957 
    811  
    812    SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, before ) 
     958   SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir ) 
    813959      !!---------------------------------------------------------------------- 
    814960      !!                  *** ROUTINE interpvn *** 
    815961      !!---------------------------------------------------------------------- 
    816       INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2 
    817       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab 
    818       LOGICAL                               , INTENT(in   ) ::   before 
    819       ! 
    820       INTEGER  ::   ji, jj, jk 
    821       REAL(wp) ::   zrhox   
    822       !!---------------------------------------------------------------------- 
     962      ! 
     963      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 
     964      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 
     965      LOGICAL, INTENT(in) :: before 
     966      INTEGER, INTENT(in) :: nb , ndir 
     967      ! 
     968      INTEGER :: ji,jj,jk 
     969      REAL(wp) :: zrhox 
     970      ! vertical interpolation: 
     971      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
     972      REAL(wp), DIMENSION(1:jpk) :: h_out 
     973      INTEGER  :: N_in, N_out, jref 
     974      REAL(wp) :: h_diff 
     975      LOGICAL  :: northern_side,southern_side 
     976      !!---------------------------------------------     
    823977      !       
    824       IF( before ) THEN       !interpv entre 1 et k2 et interpv2d en jpkp1 
    825          DO jk = k1, jpk 
    826             ptab(i1:i2,j1:j2,jk) = e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) * vn(i1:i2,j1:j2,jk) 
    827          END DO 
    828       ELSE           
    829          zrhox= Agrif_Rhox() 
     978      IF (before) THEN           
     979         DO jk=k1,k2 
     980            DO jj=j1,j2 
     981               DO ji=i1,i2 
     982                  ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk)) 
     983# if defined key_vertical 
     984                  ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk) 
     985# endif 
     986               END DO 
     987            END DO 
     988         END DO 
     989      ELSE        
     990         zrhox = Agrif_rhox() 
     991# if defined key_vertical 
     992 
     993         southern_side = (nb == 2).AND.(ndir == 1) 
     994         northern_side = (nb == 2).AND.(ndir == 2) 
     995 
     996         DO jj=j1,j2 
     997            jref = jj 
     998            IF (southern_side) jref = MAX(2,jj) 
     999            IF (northern_side) jref = MIN(nlcj-2,jj) 
     1000            DO ji=i1,i2 
     1001               N_in = 0 
     1002               DO jk=k1,k2 
     1003                  if (ptab(ji,jj,jk,2) == 0) EXIT 
     1004                  N_in = N_in + 1 
     1005                  tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 
     1006                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 
     1007               END DO 
     1008               IF (N_in == 0) THEN 
     1009                  va(ji,jj,:) = 0._wp 
     1010                  CYCLE 
     1011               ENDIF 
     1012          
     1013               N_out = 0 
     1014               DO jk=1,jpk 
     1015                  if (vmask(ji,jref,jk) == 0) EXIT 
     1016                  N_out = N_out + 1 
     1017                  h_out(N_out) = e3v_a(ji,jref,jk) 
     1018               END DO 
     1019               IF (N_out == 0) THEN 
     1020                 va(ji,jj,:) = 0._wp 
     1021                 CYCLE 
     1022               ENDIF 
     1023               call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 
     1024            END DO 
     1025         END DO 
     1026# else 
    8301027         DO jk = 1, jpkm1 
    831             va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) ) 
    832          END DO 
     1028            va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) ) 
     1029         END DO 
     1030# endif 
    8331031      ENDIF 
    8341032      !         
    8351033   END SUBROUTINE interpvn 
    836     
    8371034 
    8381035   SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir ) 
    … …  
    9551152      !!----------------------------------------------------------------------   
    9561153      IF( before ) THEN 
    957          ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2) 
     1154         IF ( ln_bt_fw ) THEN 
     1155            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2) 
     1156         ELSE 
     1157            ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2) 
     1158         ENDIF 
    9581159      ELSE 
    9591160         western_side  = (nb == 1).AND.(ndir == 1) 
    … …  
    9931194      ! 
    9941195      IF( before ) THEN 
    995          ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2) 
     1196         IF ( ln_bt_fw ) THEN 
     1197            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2) 
     1198         ELSE 
     1199            ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2) 
     1200         ENDIF 
    9961201      ELSE       
    9971202         western_side  = (nb == 1).AND.(ndir == 1) 
    … …  
    11471352 
    11481353 
    1149    SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before ) 
     1354   SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 
    11501355      !!---------------------------------------------------------------------- 
    11511356      !!                  ***  ROUTINE interavm  *** 
    11521357      !!----------------------------------------------------------------------   
    1153       INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2 
    1154       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab 
     1358      INTEGER                              , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, m1, m2 
     1359      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) ::   ptab 
    11551360      LOGICAL                              , INTENT(in   ) ::   before 
     1361      REAL(wp), DIMENSION(k1:k2) :: tabin 
     1362      REAL(wp) :: h_in(k1:k2) 
     1363      REAL(wp) :: h_out(1:jpk) 
     1364      REAL(wp) :: zrhoxy 
     1365      INTEGER  :: N_in, N_out, ji, jj, jk 
    11561366      !!----------------------------------------------------------------------   
    11571367      !       
    1158       IF( before ) THEN   ;   ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2) 
    1159       ELSE                ;   avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 
     1368      zrhoxy = Agrif_rhox()*Agrif_rhoy() 
     1369      IF (before) THEN          
     1370         DO jk=k1,k2 
     1371            DO jj=j1,j2 
     1372              DO ji=i1,i2 
     1373                    ptab(ji,jj,jk,1) = avm_k(ji,jj,jk) 
     1374              END DO 
     1375           END DO 
     1376        END DO 
     1377#ifdef key_vertical          
     1378        DO jk=k1,k2 
     1379           DO jj=j1,j2 
     1380              DO ji=i1,i2 
     1381                 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e1e2t(ji,jj) * e3w_n(ji,jj,jk)  
     1382              END DO 
     1383           END DO 
     1384        END DO 
     1385#else 
     1386      ptab(i1:i2,j1:j2,k1:k2,2) = 0._wp 
     1387#endif 
     1388      ELSE  
     1389#ifdef key_vertical          
     1390         avm_k(i1:i2,j1:j2,1:jpk) = 0. 
     1391         DO jj=j1,j2 
     1392            DO ji=i1,i2 
     1393               N_in = 0 
     1394               DO jk=k1,k2 !k2 = jpk of parent grid 
     1395                  IF (ptab(ji,jj,jk,2) == 0) EXIT 
     1396                  N_in = N_in + 1 
     1397                  tabin(jk) = ptab(ji,jj,jk,1) 
     1398                  h_in(N_in) = ptab(ji,jj,jk,2)/(e1e2t(ji,jj)*zrhoxy) 
     1399               END DO 
     1400               N_out = 0 
     1401               DO jk=1,jpk ! jpk of child grid 
     1402                  IF (wmask(ji,jj,jk) == 0) EXIT  
     1403                  N_out = N_out + 1 
     1404                  h_out(jk) = e3t_n(ji,jj,jk) 
     1405               ENDDO 
     1406               IF (N_in > 0) THEN 
     1407                  CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out) 
     1408               ENDIF 
     1409            ENDDO 
     1410         ENDDO 
     1411#else 
     1412         avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1) 
     1413#endif 
    11601414      ENDIF 
    11611415      ! 

  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90

    r9019 r9031  
    4444      ! 
    4545#if defined SPONGE 
    46       zcoef = REAL( Agrif_NbStepint(), wp ) / Agrif_rhot() 
    47       ! 
     46      !! Assume persistence: 
     47      timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot()) 
     48 
    4849      CALL Agrif_Sponge 
    4950      Agrif_SpecialValue    = 0._wp 
    … …  
    6768      ! 
    6869#if defined SPONGE 
    69       zcoef = REAL( Agrif_NbStepint(), wp ) / Agrif_rhot() 
    70       ! 
    71       Agrif_SpecialValue    = 0._wp 
     70      timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot()) 
     71 
     72      Agrif_SpecialValue=0. 
    7273      Agrif_UseSpecialValue = ln_spc_dyn 
    7374      ! 
    … …  
    189190   END SUBROUTINE Agrif_Sponge 
    190191 
    191  
    192192   SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    193193      !!---------------------------------------------------------------------- 
    … …  
    201201      INTEGER  ::   iku, ikv 
    202202      REAL(wp) :: ztsa, zabe1, zabe2, zbtr 
    203       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ztu, ztv 
    204       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff 
     203      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 
     204      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff 
     205      ! vertical interpolation: 
     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), DIMENSION(k1:k2) :: h_in 
     209      REAL(wp), DIMENSION(1:jpk) :: h_out 
     210      INTEGER :: N_in, N_out 
     211      REAL(wp) :: h_diff 
    205212      !!---------------------------------------------------------------------- 
    206213      ! 
    207214      IF( before ) THEN 
    208          tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 
     215         DO jn = 1, jpts 
     216            DO jk=k1,k2 
     217               DO jj=j1,j2 
     218                  DO ji=i1,i2 
     219                     tabres(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) 
     220                  END DO 
     221               END DO 
     222            END DO 
     223         END DO 
     224 
     225# if defined key_vertical 
     226         DO jk=k1,k2 
     227            DO jj=j1,j2 
     228               DO ji=i1,i2 
     229                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)  
     230               END DO 
     231            END DO 
     232         END DO 
     233# endif 
     234 
    209235      ELSE    
    210236         ! 
    211          tsbdiff(:,:,:,:) = tsb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:)     
     237# if defined key_vertical 
     238         tabres_child(:,:,:,:) = 0. 
     239         DO jj=j1,j2 
     240            DO ji=i1,i2 
     241               N_in = 0 
     242               DO jk=k1,k2 !k2 = jpk of parent grid 
     243                  IF (tabres(ji,jj,jk,n2) == 0) EXIT 
     244                  N_in = N_in + 1 
     245                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 
     246                  h_in(N_in) = tabres(ji,jj,jk,n2) 
     247               END DO 
     248               N_out = 0 
     249               DO jk=1,jpk ! jpk of child grid 
     250                  IF (tmask(ji,jj,jk) == 0) EXIT  
     251                  N_out = N_out + 1 
     252                  h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 
     253               ENDDO 
     254               IF (N_in > 0) THEN 
     255                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     256                  tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for????? 
     257                  DO jn=1,jpts 
     258                     call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
     259                  ENDDO 
     260               ENDIF 
     261            ENDDO 
     262         ENDDO 
     263# endif 
     264 
     265         DO jj=j1,j2 
     266            DO ji=i1,i2 
     267               DO jk=1,jpkm1 
     268# if defined key_vertical 
     269                  tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts) 
     270# else 
     271                  tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts) 
     272# endif 
     273               ENDDO 
     274            ENDDO 
     275         ENDDO 
     276 
    212277         DO jn = 1, jpts             
    213278            DO jk = 1, jpkm1 
    … …  
    256321   END SUBROUTINE interptsn_sponge 
    257322 
    258  
    259    SUBROUTINE interpun_sponge( tabres, i1, i2, j1, j2, k1, k2, before ) 
    260       !!---------------------------------------------------------------------- 
    261       !!                 *** ROUTINE interpun_sponge *** 
    262       !!---------------------------------------------------------------------- 
    263       INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2 
    264       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   tabres 
    265       LOGICAL                               , INTENT(in   ) ::   before 
    266       !! 
    267       INTEGER :: ji, jj, jk 
     323   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before) 
     324      !!--------------------------------------------- 
     325      !!   *** ROUTINE interpun_sponge *** 
     326      !!---------------------------------------------     
     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 
     329      LOGICAL, INTENT(in) :: before 
     330 
     331      INTEGER :: ji,jj,jk,jmax 
     332 
    268333      ! sponge parameters  
    269       INTEGER :: jmax 
    270       REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 
    271       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff 
    272       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff 
    273       !!---------------------------------------------------------------------- 
     334      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff 
     335      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff 
     336      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff 
     337      ! vertical interpolation: 
     338      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
     339      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
     340      REAL(wp), DIMENSION(1:jpk) :: h_out 
     341      INTEGER ::N_in,N_out 
     342      !!---------------------------------------------     
    274343      ! 
    275344      IF( before ) THEN 
    276          tabres = un(i1:i2,j1:j2,:) 
     345         DO jk=1,jpkm1 
     346            DO jj=j1,j2 
     347               DO ji=i1,i2 
     348                  tabres(ji,jj,jk,m1) = ub(ji,jj,jk) 
     349# if defined key_vertical 
     350                  tabres(ji,jj,jk,m2) = e3u_n(ji,jj,jk)*umask(ji,jj,jk) 
     351# endif 
     352               END DO 
     353            END DO 
     354         END DO 
     355 
    277356      ELSE 
    278          ubdiff(i1:i2,j1:j2,:) = ( ub(i1:i2,j1:j2,:) - tabres(:,:,:) )*umask(i1:i2,j1:j2,:) 
     357 
     358# if defined key_vertical 
     359         tabres_child(:,:,:) = 0._wp 
     360         DO jj=j1,j2 
     361            DO ji=i1,i2 
     362               N_in = 0 
     363               DO jk=k1,k2 
     364                  IF (tabres(ji,jj,jk,m2) == 0) EXIT 
     365                  N_in = N_in + 1 
     366                  tabin(jk) = tabres(ji,jj,jk,m1) 
     367                  h_in(N_in) = tabres(ji,jj,jk,m2) 
     368              ENDDO 
     369              ! 
     370              IF (N_in == 0) THEN 
     371                 tabres_child(ji,jj,:) = 0. 
     372                 CYCLE 
     373              ENDIF 
     374          
     375              N_out = 0 
     376              DO jk=1,jpk 
     377                 if (umask(ji,jj,jk) == 0) EXIT 
     378                 N_out = N_out + 1 
     379                 h_out(N_out) = e3u_n(ji,jj,jk) 
     380              ENDDO 
     381          
     382              IF (N_out == 0) THEN 
     383                 tabres_child(ji,jj,:) = 0. 
     384                 CYCLE 
     385              ENDIF 
     386          
     387              IF (N_in * N_out > 0) THEN 
     388                 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     389                 if (h_diff < -1.e4) then 
     390                    print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 
     391                 endif 
     392              ENDIF 
     393              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) 
     394          
     395            ENDDO 
     396         ENDDO 
     397 
     398         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:) 
     399#else 
     400         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 
     401#endif 
     402>>>>>>> .merge-right.r9019 
    279403         ! 
    280404         DO jk = 1, jpkm1                                 ! Horizontal slab 
    … …  
    352476   END SUBROUTINE interpun_sponge 
    353477 
    354  
    355    SUBROUTINE interpvn_sponge( tabres, i1, i2, j1, j2, k1, k2, before, nb, ndir ) 
    356       !!---------------------------------------------------------------------- 
    357       !!                 *** ROUTINE interpvn_sponge *** 
    358       !!---------------------------------------------------------------------- 
    359       INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2 
    360       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   tabres 
    361       LOGICAL                               , INTENT(in   ) ::   before 
    362       INTEGER                               , INTENT(in   ) ::   nb , ndir 
     478   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir) 
     479      !!--------------------------------------------- 
     480      !!   *** ROUTINE interpvn_sponge *** 
     481      !!---------------------------------------------  
     482      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 
     483      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres 
     484      LOGICAL, INTENT(in) :: before 
     485      INTEGER, INTENT(in) :: nb , ndir 
    363486      ! 
    364487      INTEGER ::   ji, jj, jk 
    … …  
    367490      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) ::   vbdiff, rotdiff, hdivdiff 
    368491      !!---------------------------------------------------------------------- 
    369  
     492      INTEGER  ::   ji, jj, jk, imax 
     493      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr, h_diff 
     494      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff 
     495      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff 
     496      ! vertical interpolation: 
     497      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
     498      REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 
     499      REAL(wp), DIMENSION(1:jpk) :: h_out 
     500      INTEGER :: N_in, N_out 
     501      !!---------------------------------------------  
     502>>>>>>> .merge-right.r9019 
     503       
    370504      IF( before ) THEN  
    371          tabres = vn(i1:i2,j1:j2,:) 
     505         DO jk=1,jpkm1 
     506            DO jj=j1,j2 
     507               DO ji=i1,i2 
     508                  tabres(ji,jj,jk,m1) = vb(ji,jj,jk) 
     509# if defined key_vertical 
     510                  tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_n(ji,jj,jk) 
     511# endif 
     512               END DO 
     513            END DO 
     514         END DO 
    372515      ELSE 
    373          ! 
    374          vbdiff(i1:i2,j1:j2,:) = ( vb(i1:i2,j1:j2,:) - tabres(:,:,:) ) * vmask(i1:i2,j1:j2,:) 
     516 
     517# if defined key_vertical 
     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) 
     526                  h_in(N_in) = tabres(ji,jj,jk,m2) 
     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 (vmask(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 
     551         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)   
     552# else 
     553         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 
     554# endif 
    375555         ! 
    376556         DO jk = 1, jpkm1                                 ! Horizontal slab 

  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90

    r9019 r9031  
    11#define TWO_WAY        /* TWO WAY NESTING */ 
    2 #undef DECAL_FEEDBACK /* SEPARATION of INTERFACES*/ 
     2#define DECAL_FEEDBACK /* SEPARATION of INTERFACES*/ 
     3#undef VOL_REFLUX      /* VOLUME REFLUXING*/ 
    34  
    45MODULE agrif_opa_update 
    … …  
    2324   USE in_out_manager ! I/O manager 
    2425   USE lib_mpp        ! MPP library 
     26   USE domvvl         ! Need interpolation routines  
    2527 
    2628   IMPLICIT NONE 
    … …  
    3638CONTAINS 
    3739 
    38    RECURSIVE SUBROUTINE Agrif_Update_Tra( ) 
     40   SUBROUTINE Agrif_Update_Tra( ) 
    3941      !!---------------------------------------------------------------------- 
    4042      !!                   *** ROUTINE Agrif_Update_Tra *** 
    … …  
    4446      ! 
    4547#if defined TWO_WAY   
    46       IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed(), 'nbcline', nbcline 
     48      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed() 
    4749 
    4850      Agrif_UseSpecialValueInUpdate = .TRUE. 
    4951      Agrif_SpecialValueFineGrid    = 0._wp 
    5052      !  
    51       IF (MOD(nbcline,nbclineupdate) == 0) THEN 
    5253# if ! defined DECAL_FEEDBACK 
    53          CALL Agrif_Update_Variable(tsn_id, procname=updateTS) 
     54      CALL Agrif_Update_Variable(tsn_id, procname=updateTS) 
     55! near boundary update: 
     56!      CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS) 
    5457# else 
    55          CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS) 
     58      CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS) 
     59! near boundary update: 
     60!      CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS) 
    5661# endif 
    57       ELSE 
    58 # if ! defined DECAL_FEEDBACK 
    59          CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS) 
    60 # else 
    61          CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS) 
    62 # endif 
    63       ENDIF 
    6462      ! 
    6563      Agrif_UseSpecialValueInUpdate = .FALSE. 
    6664      ! 
    67       IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update: 
    68          CALL Agrif_ChildGrid_To_ParentGrid() 
    69          CALL Agrif_Update_Tra() 
    70          CALL Agrif_ParentGrid_To_ChildGrid() 
    71       ENDIF 
    72       ! 
    7365#endif 
    7466      ! 
    7567   END SUBROUTINE Agrif_Update_Tra 
    7668 
    77  
    78    RECURSIVE SUBROUTINE Agrif_Update_Dyn( ) 
     69   SUBROUTINE Agrif_Update_Dyn( ) 
    7970      !!---------------------------------------------------------------------- 
    8071      !!                   *** ROUTINE Agrif_Update_Dyn *** 
    … …  
    8475      ! 
    8576#if defined TWO_WAY 
    86       IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed(), 'nbcline', nbcline 
     77      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update momentum from grid Number',Agrif_Fixed() 
    8778 
    8879      Agrif_UseSpecialValueInUpdate = .FALSE. 
    8980      Agrif_SpecialValueFineGrid = 0. 
    9081      !      
    91       IF (mod(nbcline,nbclineupdate) == 0) THEN 
    9282# if ! defined DECAL_FEEDBACK 
    93          CALL Agrif_Update_Variable(un_update_id,procname = updateU) 
    94          CALL Agrif_Update_Variable(vn_update_id,procname = updateV) 
     83      CALL Agrif_Update_Variable(un_update_id,procname = updateU) 
     84      CALL Agrif_Update_Variable(vn_update_id,procname = updateV) 
     85! near boundary update: 
     86!      CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU) 
     87!      CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV) 
    9588# else 
    96          CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU) 
    97          CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV) 
     89      CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU) 
     90      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV) 
     91! near boundary update: 
     92!      CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU) 
     93!      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV) 
    9894# endif 
    99       ELSE 
    100 # if ! defined DECAL_FEEDBACK 
    101          CALL Agrif_Update_Variable(un_update_id,locupdate=(/0,1/),procname = updateU) 
    102          CALL Agrif_Update_Variable(vn_update_id,locupdate=(/0,1/),procname = updateV)          
    103 # else 
    104          CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateU) 
    105          CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV) 
    106 # endif 
    107       ENDIF 
    10895 
    10996# if ! defined DECAL_FEEDBACK 
    … …  
    114101      CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)   
    115102# endif 
    116  
     103      ! 
     104# if ! defined DECAL_FEEDBACK 
     105      ! Account for updated thicknesses at boundary edges 
     106      IF (.NOT.ln_linssh) THEN 
     107!         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy) 
     108!         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy) 
     109      ENDIF 
     110# endif 
     111      !  
    117112      IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
    118113         ! Update time integrated transports 
    119          IF (mod(nbcline,nbclineupdate) == 0) THEN 
    120114#  if ! defined DECAL_FEEDBACK 
    121             CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b) 
    122             CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b) 
     115         CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b) 
     116         CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b) 
    123117#  else 
    124             CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b) 
    125             CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b) 
     118         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b) 
     119         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b) 
    126120#  endif 
    127          ELSE 
    128 #  if ! defined DECAL_FEEDBACK 
    129             CALL Agrif_Update_Variable(ub2b_update_id,locupdate=(/0,1/),procname = updateub2b) 
    130             CALL Agrif_Update_Variable(vb2b_update_id,locupdate=(/0,1/),procname = updatevb2b) 
    131 #  else 
    132             CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,1/),locupdate2=(/1,1/),procname = updateub2b) 
    133             CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updatevb2b) 
    134 #  endif 
    135          ENDIF 
    136121      END IF 
    137       ! 
    138       nbcline = nbcline + 1 
     122#endif 
     123      ! 
     124   END SUBROUTINE Agrif_Update_Dyn 
     125 
     126   SUBROUTINE Agrif_Update_ssh( ) 
     127      !!--------------------------------------------- 
     128      !!   *** ROUTINE Agrif_Update_ssh *** 
     129      !!--------------------------------------------- 
     130      !  
     131      IF (Agrif_Root()) RETURN 
     132      ! 
     133#if defined TWO_WAY 
    139134      ! 
    140135      Agrif_UseSpecialValueInUpdate = .TRUE. 
    … …  
    145140      CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH) 
    146141# endif 
     142      ! 
    147143      Agrif_UseSpecialValueInUpdate = .FALSE. 
    148       !  
     144      ! 
     145#  if defined DECAL_FEEDBACK && defined VOL_REFLUX 
     146      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 
     147         ! Refluxing on ssh: 
     148         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu) 
     149         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv) 
     150      END IF 
     151#  endif 
     152      ! 
    149153#endif 
    150154      ! 
    151       ! Do recursive update: 
    152       IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update: 
    153          CALL Agrif_ChildGrid_To_ParentGrid() 
    154          CALL Agrif_Update_Dyn() 
    155          CALL Agrif_ParentGrid_To_ChildGrid() 
    156       ENDIF 
    157       ! 
    158    END SUBROUTINE Agrif_Update_Dyn 
    159  
     155   END SUBROUTINE Agrif_Update_ssh 
     156 
     157 
     158   SUBROUTINE Agrif_Update_Tke( kt ) 
     159      !!--------------------------------------------- 
     160      !!   *** ROUTINE Agrif_Update_Tke *** 
     161      !!--------------------------------------------- 
     162      !! 
     163      INTEGER, INTENT(in) :: kt  
     164      !  
     165      IF (Agrif_Root()) RETURN 
     166      !        
     167      IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN 
     168#  if defined TWO_WAY 
     169 
     170      Agrif_UseSpecialValueInUpdate = .TRUE. 
     171      Agrif_SpecialValueFineGrid = 0. 
     172 
     173      CALL Agrif_Update_Variable( en_id, locupdate=(/0,0/), procname=updateEN  ) 
     174      CALL Agrif_Update_Variable(avt_id, locupdate=(/0,0/), procname=updateAVT ) 
     175      CALL Agrif_Update_Variable(avm_id, locupdate=(/0,0/), procname=updateAVM ) 
     176 
     177      Agrif_UseSpecialValueInUpdate = .FALSE. 
     178 
     179#  endif 
     180       
     181   END SUBROUTINE Agrif_Update_Tke 
     182 
     183 
     184   SUBROUTINE Agrif_Update_vvl( ) 
     185      !!--------------------------------------------- 
     186      !!   *** ROUTINE Agrif_Update_vvl *** 
     187      !!--------------------------------------------- 
     188      ! 
     189      IF (Agrif_Root()) RETURN 
     190      ! 
     191#if defined TWO_WAY   
     192      ! 
     193      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step() 
     194      ! 
     195      Agrif_UseSpecialValueInUpdate = .TRUE. 
     196      Agrif_SpecialValueFineGrid = 0. 
     197      !  
     198      ! No interface separation here, update vertical grid at T points  
     199      ! everywhere over the overlapping regions (one account for refluxing in that case): 
     200      CALL Agrif_Update_Variable(e3t_id, procname=updatee3t)  
     201      ! 
     202      Agrif_UseSpecialValueInUpdate = .FALSE. 
     203      ! 
     204      CALL Agrif_ChildGrid_To_ParentGrid() 
     205      CALL dom_vvl_update_UVF 
     206      CALL Agrif_ParentGrid_To_ChildGrid() 
     207      ! 
     208#endif 
     209      ! 
     210   END SUBROUTINE Agrif_Update_vvl 
     211 
     212   SUBROUTINE dom_vvl_update_UVF 
     213      !!--------------------------------------------- 
     214      !!       *** ROUTINE dom_vvl_update_UVF *** 
     215      !!--------------------------------------------- 
     216      !! 
     217      INTEGER :: jk 
     218      REAL(wp):: zcoef 
     219      !!--------------------------------------------- 
     220 
     221      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', & 
     222                  & Agrif_Fixed(), 'Step', Agrif_Nb_Step() 
     223 
     224      ! Save "old" scale factor (prior update) for subsequent asselin correction 
     225      ! of prognostic variables 
     226      ! ----------------------- 
     227      ! 
     228      e3u_a(:,:,:) = e3u_n(:,:,:) 
     229      e3v_a(:,:,:) = e3v_n(:,:,:) 
     230!      ua(:,:,:) = e3u_b(:,:,:) 
     231!      va(:,:,:) = e3v_b(:,:,:) 
     232      hu_a(:,:) = hu_n(:,:) 
     233      hv_a(:,:) = hv_n(:,:) 
     234 
     235      ! 1) NOW fields 
     236      !-------------- 
     237       
     238         ! Vertical scale factor interpolations 
     239         ! ------------------------------------ 
     240      CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:) ,  'U' ) 
     241      CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:) ,  'V' ) 
     242      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:) ,  'F' ) 
     243 
     244      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 
     245      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 
     246 
     247         ! Update total depths: 
     248         ! -------------------- 
     249      hu_n(:,:) = 0._wp                        ! Ocean depth at U-points 
     250      hv_n(:,:) = 0._wp                        ! Ocean depth at V-points 
     251      DO jk = 1, jpkm1 
     252         hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 
     253         hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 
     254      END DO 
     255      !                                        ! Inverse of the local depth 
     256      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) 
     257      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 
     258 
     259 
     260      ! 2) BEFORE fields: 
     261      !------------------ 
     262!      IF (     (.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_exp)) & 
     263!         & .OR.(.NOT.(lk_agrif_fstep.AND.(neuler==0)).AND.(ln_dynspg_ts    & 
     264!         & .AND.(.NOT.ln_bt_fw)))) THEN 
     265      IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
     266         ! 
     267         ! Vertical scale factor interpolations 
     268         ! ------------------------------------ 
     269         CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:),  'U'  ) 
     270         CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:),  'V'  ) 
     271 
     272         CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 
     273         CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 
     274 
     275         ! Update total depths: 
     276         ! -------------------- 
     277         hu_b(:,:) = 0._wp                     ! Ocean depth at U-points 
     278         hv_b(:,:) = 0._wp                     ! Ocean depth at V-points 
     279         DO jk = 1, jpkm1 
     280            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 
     281            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 
     282         END DO 
     283         !                                     ! Inverse of the local depth 
     284         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 
     285         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 
     286      ENDIF 
     287      ! 
     288   END SUBROUTINE dom_vvl_update_UVF 
     289 
     290#if defined key_vertical 
    160291 
    161292   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    162293      !!---------------------------------------------------------------------- 
    163294      !!           *** ROUTINE updateT *** 
    164       !!---------------------------------------------------------------------- 
    165       INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2 
    166       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres 
    167       LOGICAL                                    , INTENT(in   ) ::   before 
    168       ! 
    169       INTEGER :: ji, jj, jk, jn 
    170       !!---------------------------------------------------------------------- 
    171       ! 
    172       IF( before ) THEN 
    173          DO jn = n1, n2 
    174             DO jk = k1, k2 
    175                DO jj = j1, j2 
    176                   DO ji = i1, i2 
    177                      tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 
     295      !!--------------------------------------------- 
     296      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
     297      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     298      LOGICAL, INTENT(in) :: before 
     299      !! 
     300      INTEGER :: ji,jj,jk,jn 
     301      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 
     302      REAL(wp) :: h_in(k1:k2) 
     303      REAL(wp) :: h_out(1:jpk) 
     304      INTEGER  :: N_in, N_out 
     305      REAL(wp) :: h_diff 
     306      REAL(wp) :: zrho_xy 
     307      REAL(wp) :: tabin(k1:k2,n1:n2) 
     308      !!--------------------------------------------- 
     309      ! 
     310      IF (before) THEN 
     311         AGRIF_SpecialValue = -999._wp 
     312         zrho_xy = Agrif_rhox() * Agrif_rhoy()  
     313         DO jn = n1,n2-1 
     314            DO jk=k1,k2 
     315               DO jj=j1,j2 
     316                  DO ji=i1,i2 
     317                     tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) & 
     318                                           * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp 
    178319                  END DO 
    179320               END DO 
    180321            END DO 
    181322         END DO 
     323         DO jk=k1,k2 
     324            DO jj=j1,j2 
     325               DO ji=i1,i2 
     326                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) & 
     327                                           + (tmask(ji,jj,jk)-1)*999._wp 
     328               END DO 
     329            END DO 
     330         END DO 
    182331      ELSE 
     332         tabres_child(:,:,:,:) = 0. 
     333         AGRIF_SpecialValue = 0._wp 
     334         DO jj=j1,j2 
     335            DO ji=i1,i2 
     336               N_in = 0 
     337               DO jk=k1,k2 !k2 = jpk of child grid 
     338                  IF (tabres(ji,jj,jk,n2) == 0  ) EXIT 
     339                  N_in = N_in + 1 
     340                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     341                  h_in(N_in) = tabres(ji,jj,jk,n2) 
     342               ENDDO 
     343               N_out = 0 
     344               DO jk=1,jpk ! jpk of parent grid 
     345                  IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 
     346                  N_out = N_out + 1 
     347                  h_out(N_out) = e3t_n(ji,jj,jk)  
     348               ENDDO 
     349               IF (N_in > 0) THEN !Remove this? 
     350                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     351                  IF (h_diff < -1.e-4) THEN 
     352                     print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 
     353                     print *,h_in(1:N_in) 
     354                     print *,h_out(1:N_out) 
     355                     STOP 
     356                  ENDIF 
     357                  DO jn=n1,n2-1 
     358                     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) 
     359                  ENDDO 
     360               ENDIF 
     361            ENDDO 
     362         ENDDO 
     363 
    183364         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    184365            ! Add asselin part 
    185             DO jn = n1,n2 
     366            DO jn = n1,n2-1 
     367               DO jk=1,jpk 
     368                  DO jj=j1,j2 
     369                     DO ji=i1,i2 
     370                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
     371                           tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
     372                                 & + atfp * ( tabres_child(ji,jj,jk,jn) & 
     373                                 &          - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
     374                        ENDIF 
     375                     ENDDO 
     376                  ENDDO 
     377               ENDDO 
     378            ENDDO 
     379         ENDIF 
     380         DO jn = n1,n2-1 
     381            DO jk=1,jpk 
     382               DO jj=j1,j2 
     383                  DO ji=i1,i2 
     384                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN  
     385                        tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 
     386                     END IF 
     387                  END DO 
     388               END DO 
     389            END DO 
     390         END DO 
     391      ENDIF 
     392      !  
     393   END SUBROUTINE updateTS 
     394 
     395# else 
     396 
     397   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     398      !!--------------------------------------------- 
     399      !!           *** ROUTINE updateT *** 
     400      !!--------------------------------------------- 
     401      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
     402      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     403      LOGICAL, INTENT(in) :: before 
     404      !! 
     405      INTEGER :: ji,jj,jk,jn 
     406      REAL(wp) :: ztb, ztnu, ztno 
     407      !!--------------------------------------------- 
     408      ! 
     409      IF (before) THEN 
     410         DO jn = 1,jpts 
     411            DO jk=k1,k2 
     412               DO jj=j1,j2 
     413                  DO ji=i1,i2 
     414!> jc tmp 
     415                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk) 
     416!                     tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) 
     417!< jc tmp 
     418                  END DO 
     419               END DO 
     420            END DO 
     421         END DO 
     422      ELSE 
     423!> jc tmp 
     424         DO jn = 1,jpts 
     425            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 
     426                                         & * tmask(i1:i2,j1:j2,k1:k2) 
     427         ENDDO 
     428!< jc tmp 
     429         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     430            ! Add asselin part 
     431            DO jn = 1,jpts 
    186432               DO jk = k1, k2 
    187433                  DO jj = j1, j2 
    188434                     DO ji = i1, i2 
    189435                        IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 
    190                            tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
    191                               &             + atfp * ( tabres(ji,jj,jk,jn) - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
     436                           ztb  = tsb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used 
     437                           ztnu = tabres(ji,jj,jk,jn) 
     438                           ztno = tsn(ji,jj,jk,jn) * e3t_a(ji,jj,jk) 
     439                           tsb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  &  
     440                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk) 
    192441                        ENDIF 
    193442                     END DO 
    … …  
    196445            END DO 
    197446         ENDIF 
    198          DO jn = n1,n2 
     447         DO jn = 1,jpts 
    199448            DO jk=k1,k2 
    200449               DO jj=j1,j2 
    201450                  DO ji=i1,i2 
    202451                     IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN  
    203                         tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) * tmask(ji,jj,jk) 
     452                        tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk) 
    204453                     END IF 
    205454                  END DO 
    … …  
    207456            END DO 
    208457         END DO 
     458         ! 
     459         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     460            tsb(i1:i2,j1:j2,k1:k2,1:jpts)  = tsn(i1:i2,j1:j2,k1:k2,1:jpts) 
     461         ENDIF 
     462         ! 
    209463      ENDIF 
    210464      !  
    211465   END SUBROUTINE updateTS 
    212466 
    213  
    214    SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before ) 
     467# endif 
     468 
     469# if defined key_vertical 
     470 
     471   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    215472      !!--------------------------------------------- 
    216473      !!           *** ROUTINE updateu *** 
    217474      !!--------------------------------------------- 
    218       INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2 
    219       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 
    220       LOGICAL                               , INTENT(in   ) :: before 
     475      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     476      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     477      LOGICAL                                     , INTENT(in   ) :: before 
    221478      ! 
    222479      INTEGER ::   ji, jj, jk 
    223480      REAL(wp)::   zrhoy 
     481! VERTICAL REFINEMENT BEGIN 
     482      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
     483      REAL(wp) :: h_in(k1:k2) 
     484      REAL(wp) :: h_out(1:jpk) 
     485      INTEGER  :: N_in, N_out 
     486      REAL(wp) :: h_diff, excess, thick 
     487      REAL(wp) :: tabin(k1:k2) 
     488! VERTICAL REFINEMENT END 
     489      !!--------------------------------------------- 
     490      !  
     491      IF( before ) THEN 
     492         zrhoy = Agrif_Rhoy() 
     493         AGRIF_SpecialValue = -999._wp 
     494         DO jk=k1,k2 
     495            DO jj=j1,j2 
     496               DO ji=i1,i2 
     497                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk)  & 
     498                                       + (umask(ji,jj,jk)-1)*999._wp 
     499                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)  & 
     500                                       + (umask(ji,jj,jk)-1)*999._wp 
     501               END DO 
     502            END DO 
     503         END DO 
     504      ELSE 
     505         tabres_child(:,:,:) = 0. 
     506         AGRIF_SpecialValue = 0._wp 
     507         DO jj=j1,j2 
     508            DO ji=i1,i2 
     509               N_in = 0 
     510               h_in(:) = 0._wp 
     511               tabin(:) = 0._wp 
     512               DO jk=k1,k2 !k2=jpk of child grid 
     513                  IF( tabres(ji,jj,jk,2) < -900) EXIT 
     514                  N_in = N_in + 1 
     515                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     516                  h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 
     517               ENDDO 
     518               N_out = 0 
     519               DO jk=1,jpk 
     520                  IF (umask(ji,jj,jk) == 0) EXIT 
     521                  N_out = N_out + 1 
     522                  h_out(N_out) = e3u_n(ji,jj,jk) 
     523               ENDDO 
     524               IF (N_in * N_out > 0) THEN 
     525                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     526                  IF (h_diff < -1.e-4) THEN 
     527!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
     528!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 
     529                     excess = 0._wp 
     530                     DO jk=N_in,1,-1 
     531                        thick = MIN(-1*h_diff, h_in(jk)) 
     532                        excess = excess + tabin(jk)*thick*e2u(ji,jj) 
     533                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 
     534                        h_diff = h_diff + thick 
     535                        IF ( h_diff == 0) THEN 
     536                           N_in = jk 
     537                           h_in(jk) = h_in(jk) - thick 
     538                           EXIT 
     539                        ENDIF 
     540                     ENDDO 
     541                  ENDIF 
     542                  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) 
     543                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 
     544               ENDIF 
     545            ENDDO 
     546         ENDDO 
     547 
     548         DO jk=1,jpk 
     549            DO jj=j1,j2 
     550               DO ji=i1,i2 
     551                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     552                     ub(ji,jj,jk) = ub(ji,jj,jk) &  
     553                           & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
     554                  ENDIF 
     555                  ! 
     556                  un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 
     557               END DO 
     558            END DO 
     559         END DO 
     560      ENDIF 
     561      !  
     562   END SUBROUTINE updateu 
     563 
     564#else 
     565 
     566   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     567      !!--------------------------------------------- 
     568      !!           *** ROUTINE updateu *** 
     569      !!--------------------------------------------- 
     570      INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     571      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     572      LOGICAL                                     , INTENT(in   ) :: before 
     573      ! 
     574      INTEGER  :: ji, jj, jk 
     575      REAL(wp) :: zrhoy, zub, zunu, zuno 
    224576      !!--------------------------------------------- 
    225577      !  
    … …  
    227579         zrhoy = Agrif_Rhoy() 
    228580         DO jk = k1, k2 
    229             tabres(i1:i2,j1:j2,jk) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 
     581            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) 
    230582         END DO 
    231583      ELSE 
    … …  
    233585            DO jj=j1,j2 
    234586               DO ji=i1,i2 
    235                   tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) 
     587                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj)  
    236588                  ! 
    237589                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    238                      ub(ji,jj,jk) = ub(ji,jj,jk) &  
    239                            & + atfp * ( tabres(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 
     590                     zub  = ub(ji,jj,jk) * e3u_b(ji,jj,jk)  ! fse3t_b prior update should be used 
     591                     zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk) 
     592                     zunu = tabres(ji,jj,jk,1) 
     593                     ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) &       
     594                                    & * umask(ji,jj,jk) / e3u_b(ji,jj,jk) 
    240595                  ENDIF 
    241596                  ! 
    242                   un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) 
    243                END DO 
    244             END DO 
    245          END DO 
     597                  un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk) 
     598               END DO 
     599            END DO 
     600         END DO 
     601         ! 
     602         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     603            ub(i1:i2,j1:j2,k1:k2)  = un(i1:i2,j1:j2,k1:k2) 
     604         ENDIF 
     605         ! 
    246606      ENDIF 
    247607      !  
    248608   END SUBROUTINE updateu 
    249609 
    250  
    251    SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before ) 
    252       !!---------------------------------------------------------------------- 
    253       !!                      *** ROUTINE updatev *** 
    254       !!---------------------------------------------------------------------- 
    255       INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2 
    256       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 
    257       LOGICAL                               , INTENT(in   ) :: before 
     610# endif 
     611 
     612   SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
     613      !!--------------------------------------------- 
     614      !!           *** ROUTINE correct_u_bdy *** 
     615      !!--------------------------------------------- 
     616      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     617      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     618      LOGICAL                                     , INTENT(in   ) :: before 
     619      INTEGER                                     , INTENT(in)    :: nb, ndir 
    258620      !! 
     621      LOGICAL :: western_side, eastern_side  
     622      ! 
     623      INTEGER  :: jj, jk 
     624      REAL(wp) :: zcor 
     625      !!--------------------------------------------- 
     626      !  
     627      IF( .NOT.before ) THEN 
     628         ! 
     629         western_side  = (nb == 1).AND.(ndir == 1) 
     630         eastern_side  = (nb == 1).AND.(ndir == 2) 
     631         ! 
     632         IF (western_side) THEN 
     633            DO jj=j1,j2 
     634               zcor = un_b(i1-1,jj) * hu_a(i1-1,jj) * r1_hu_n(i1-1,jj) - un_b(i1-1,jj) 
     635               un_b(i1-1,jj) = un_b(i1-1,jj) + zcor 
     636               DO jk=1,jpkm1 
     637                  un(i1-1,jj,jk) = un(i1-1,jj,jk) + zcor * umask(i1-1,jj,jk) 
     638               END DO  
     639            END DO 
     640         ENDIF 
     641         ! 
     642         IF (eastern_side) THEN 
     643            DO jj=j1,j2 
     644               zcor = un_b(i2+1,jj) * hu_a(i2+1,jj) * r1_hu_n(i2+1,jj) - un_b(i2+1,jj) 
     645               un_b(i2+1,jj) = un_b(i2+1,jj) + zcor 
     646               DO jk=1,jpkm1 
     647                  un(i2+1,jj,jk) = un(i2+1,jj,jk) + zcor * umask(i2+1,jj,jk) 
     648               END DO  
     649            END DO 
     650         ENDIF 
     651         ! 
     652      ENDIF 
     653      !  
     654   END SUBROUTINE correct_u_bdy 
     655 
     656# if  defined key_vertical 
     657 
     658   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     659      !!--------------------------------------------- 
     660      !!           *** ROUTINE updatev *** 
     661      !!--------------------------------------------- 
     662      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     663      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     664      LOGICAL                                     , INTENT(in   ) :: before 
     665      ! 
    259666      INTEGER  ::   ji, jj, jk 
    260667      REAL(wp) ::   zrhox 
    261       !!---------------------------------------------------------------------- 
     668! VERTICAL REFINEMENT BEGIN 
     669      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
     670      REAL(wp) :: h_in(k1:k2) 
     671      REAL(wp) :: h_out(1:jpk) 
     672      INTEGER :: N_in, N_out 
     673      REAL(wp) :: h_diff, excess, thick 
     674      REAL(wp) :: tabin(k1:k2) 
     675! VERTICAL REFINEMENT END 
     676      !!---------------------------------------------       
    262677      ! 
    263678      IF( before ) THEN 
    264679         zrhox = Agrif_Rhox() 
    265          DO jk = k1, k2 
    266             DO jj = j1, j2 
    267                DO ji = i1, i2 
    268                   tabres(ji,jj,jk) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
     680         AGRIF_SpecialValue = -999._wp 
     681         DO jk=k1,k2 
     682            DO jj=j1,j2 
     683               DO ji=i1,i2 
     684                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) & 
     685                                       + (vmask(ji,jj,jk)-1)*999._wp 
     686                  tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) & 
     687                                       + (vmask(ji,jj,jk)-1)*999._wp 
    269688               END DO 
    270689            END DO 
    271690         END DO 
    272691      ELSE 
    273          DO jk = k1, k2 
    274             DO jj = j1, j2 
    275                DO ji = i1, i2 
    276                   tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) 
     692         tabres_child(:,:,:) = 0. 
     693         AGRIF_SpecialValue = 0._wp 
     694         DO jj=j1,j2 
     695            DO ji=i1,i2 
     696               N_in = 0 
     697               DO jk=k1,k2 
     698                  IF (tabres(ji,jj,jk,2) < -900) EXIT 
     699                  N_in = N_in + 1 
     700                  tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     701                  h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) 
     702               ENDDO 
     703               N_out = 0 
     704               DO jk=1,jpk 
     705                  IF (vmask(ji,jj,jk) == 0) EXIT 
     706                  N_out = N_out + 1 
     707                  h_out(N_out) = e3v_n(ji,jj,jk) 
     708               ENDDO 
     709               IF (N_in * N_out > 0) THEN 
     710                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     711                  IF (h_diff < -1.e-4) then 
     712!Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
     713!In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 
     714                     excess = 0._wp 
     715                     DO jk=N_in,1,-1 
     716                        thick = MIN(-1*h_diff, h_in(jk)) 
     717                        excess = excess + tabin(jk)*thick*e2u(ji,jj) 
     718                        tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 
     719                        h_diff = h_diff + thick 
     720                        IF ( h_diff == 0) THEN 
     721                           N_in = jk 
     722                           h_in(jk) = h_in(jk) - thick 
     723                           EXIT 
     724                        ENDIF 
     725                     ENDDO 
     726                  ENDIF 
     727                  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) 
     728                  tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 
     729               ENDIF 
     730            ENDDO 
     731         ENDDO 
     732 
     733         DO jk=1,jpk 
     734            DO jj=j1,j2 
     735               DO ji=i1,i2 
    277736                  ! 
    278737                  IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN ! Add asselin part 
    279738                     vb(ji,jj,jk) = vb(ji,jj,jk) &  
    280                         &         + atfp * ( tabres(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
     739                           & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 
    281740                  ENDIF 
    282741                  ! 
    283                   vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) 
     742                  vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) 
    284743               END DO 
    285744            END DO 
    … …  
    288747      !  
    289748   END SUBROUTINE updatev 
     749 
     750# else 
     751 
     752   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) 
     753      !!--------------------------------------------- 
     754      !!           *** ROUTINE updatev *** 
     755      !!--------------------------------------------- 
     756      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     757      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     758      LOGICAL                                     , INTENT(in   ) :: before 
     759      ! 
     760      INTEGER  :: ji, jj, jk 
     761      REAL(wp) :: zrhox, zvb, zvnu, zvno 
     762      !!---------------------------------------------       
     763      ! 
     764      IF (before) THEN 
     765         zrhox = Agrif_Rhox() 
     766         DO jk=k1,k2 
     767            DO jj=j1,j2 
     768               DO ji=i1,i2 
     769                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
     770               END DO 
     771            END DO 
     772         END DO 
     773      ELSE 
     774         DO jk=k1,k2 
     775            DO jj=j1,j2 
     776               DO ji=i1,i2 
     777                  tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) 
     778                  ! 
     779                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
     780                     zvb  = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used 
     781                     zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk) 
     782                     zvnu = tabres(ji,jj,jk,1) 
     783                     vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) &       
     784                                    & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk) 
     785                  ENDIF 
     786                  ! 
     787                  vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk) 
     788               END DO 
     789            END DO 
     790         END DO 
     791         ! 
     792         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     793            vb(i1:i2,j1:j2,k1:k2)  = vn(i1:i2,j1:j2,k1:k2) 
     794         ENDIF 
     795         ! 
     796      ENDIF 
     797      !  
     798   END SUBROUTINE updatev 
     799 
     800# endif 
     801 
     802   SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
     803      !!--------------------------------------------- 
     804      !!           *** ROUTINE correct_u_bdy *** 
     805      !!--------------------------------------------- 
     806      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
     807      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     808      LOGICAL                                     , INTENT(in   ) :: before 
     809      INTEGER                                     , INTENT(in)    :: nb, ndir 
     810      !! 
     811      LOGICAL :: southern_side, northern_side  
     812      ! 
     813      INTEGER  :: ji, jk 
     814      REAL(wp) :: zcor 
     815      !!--------------------------------------------- 
     816      !  
     817      IF( .NOT.before ) THEN 
     818         ! 
     819         southern_side = (nb == 2).AND.(ndir == 1) 
     820         northern_side = (nb == 2).AND.(ndir == 2) 
     821         ! 
     822         IF (southern_side) THEN 
     823            DO ji=i1,i2 
     824               zcor = vn_b(ji,j1-1) * hv_a(ji,j1-1) * r1_hv_n(ji,j1-1) - vn_b(ji,j1-1) 
     825               vn_b(ji,j1-1) = vn_b(ji,j1-1) + zcor 
     826               DO jk=1,jpkm1 
     827                  vn(ji,j1-1,jk) = vn(ji,j1-1,jk) + zcor * vmask(ji,j1-1,jk) 
     828               END DO  
     829            END DO 
     830         ENDIF 
     831         ! 
     832         IF (northern_side) THEN 
     833            DO ji=i1,i2 
     834               zcor = vn_b(ji,j2+1) * hv_a(ji,j2+1) * r1_hv_n(ji,j2+1) - vn_b(ji,j2+1) 
     835               vn_b(ji,j2+1) = vn_b(ji,j2+1) + zcor 
     836               DO jk=1,jpkm1 
     837                  vn(ji,j2+1,jk) = vn(ji,j2+1,jk) + zcor * vmask(ji,j2+1,jk) 
     838               END DO  
     839            END DO 
     840         ENDIF 
     841         ! 
     842      ENDIF 
     843      !  
     844   END SUBROUTINE correct_v_bdy 
    290845 
    291846 
    … …  
    298853      LOGICAL                         , INTENT(in   ) ::   before 
    299854      !!  
    300       INTEGER ::   ji, jj, jk 
    301       REAL(wp)::   zrhoy, zcorr 
     855      INTEGER  :: ji, jj, jk 
     856      REAL(wp) :: zrhoy 
     857      REAL(wp) :: zcorr 
    302858      !!--------------------------------------------- 
    303859      ! 
    … …  
    312868         DO jj=j1,j2 
    313869            DO ji=i1,i2 
    314                tabres(ji,jj) =  tabres(ji,jj) * r1_hu_n(ji,jj) * r1_e2u(ji,jj)   
     870               tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj)   
    315871               !     
    316872               ! Update "now" 3d velocities: 
    … …  
    319875                  spgu(ji,jj) = spgu(ji,jj) + e3u_n(ji,jj,jk) * un(ji,jj,jk) 
    320876               END DO 
    321                spgu(ji,jj) = spgu(ji,jj) * r1_hu_n(ji,jj) 
    322877               ! 
    323                zcorr = tabres(ji,jj) - spgu(ji,jj) 
     878               zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu_n(ji,jj) 
    324879               DO jk=1,jpkm1               
    325880                  un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk)            
    … …  
    329884               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    330885                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    331                      zcorr = tabres(ji,jj) - un_b(ji,jj) 
     886                     zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * r1_hu_b(ji,jj) 
    332887                     ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1) 
    333888                  END IF 
    334                ENDIF              
    335                un_b(ji,jj) = tabres(ji,jj) * umask(ji,jj,1) 
     889               ENDIF     
     890               un_b(ji,jj) = tabres(ji,jj) * r1_hu_n(ji,jj) * umask(ji,jj,1) 
    336891               !        
    337892               ! Correct "before" velocities to hold correct bt component: 
    … …  
    340895                  spgu(ji,jj) = spgu(ji,jj) + e3u_b(ji,jj,jk) * ub(ji,jj,jk) 
    341896               END DO 
    342                spgu(ji,jj) = spgu(ji,jj) * r1_hu_b(ji,jj) 
    343897               ! 
    344                zcorr = ub_b(ji,jj) - spgu(ji,jj) 
     898               zcorr = ub_b(ji,jj) - spgu(ji,jj) * r1_hu_b(ji,jj) 
    345899               DO jk=1,jpkm1               
    346900                  ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk)            
    … …  
    349903            END DO 
    350904         END DO 
     905         ! 
     906         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     907            ub_b(i1:i2,j1:j2)  = un_b(i1:i2,j1:j2) 
     908         ENDIF 
    351909      ENDIF 
    352910      ! 
    … …  
    362920      LOGICAL                         , INTENT(in   ) ::   before 
    363921      !  
    364       INTEGER :: ji, jj, jk 
     922      INTEGER  :: ji, jj, jk 
    365923      REAL(wp) :: zrhox, zcorr 
    366924      !!---------------------------------------------------------------------- 
    … …  
    376934         DO jj=j1,j2 
    377935            DO ji=i1,i2 
    378                tabres(ji,jj) =  tabres(ji,jj) * r1_hv_n(ji,jj) * r1_e1v(ji,jj)   
     936               tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj)   
    379937               !     
    380938               ! Update "now" 3d velocities: 
    … …  
    383941                  spgv(ji,jj) = spgv(ji,jj) + e3v_n(ji,jj,jk) * vn(ji,jj,jk) 
    384942               END DO 
    385                spgv(ji,jj) = spgv(ji,jj) * r1_hv_n(ji,jj) 
    386943               ! 
    387                zcorr = tabres(ji,jj) - spgv(ji,jj) 
     944               zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv_n(ji,jj) 
    388945               DO jk=1,jpkm1               
    389946                  vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk)            
    … …  
    393950               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    394951                  IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 
    395                      zcorr = tabres(ji,jj) - vn_b(ji,jj) 
     952                     zcorr = (tabres(ji,jj) - vn_b(ji,jj) * hv_a(ji,jj)) * r1_hv_b(ji,jj) 
    396953                     vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1) 
    397954                  END IF 
    398955               ENDIF               
    399                vn_b(ji,jj) = tabres(ji,jj) * vmask(ji,jj,1) 
     956               vn_b(ji,jj) = tabres(ji,jj) * r1_hv_n(ji,jj) * vmask(ji,jj,1) 
    400957               !        
    401958               ! Correct "before" velocities to hold correct bt component: 
    … …  
    404961                  spgv(ji,jj) = spgv(ji,jj) + e3v_b(ji,jj,jk) * vb(ji,jj,jk) 
    405962               END DO 
    406                spgv(ji,jj) = spgv(ji,jj) * r1_hv_b(ji,jj) 
    407963               ! 
    408                zcorr = vb_b(ji,jj) - spgv(ji,jj) 
     964               zcorr = vb_b(ji,jj) - spgv(ji,jj) * r1_hv_b(ji,jj) 
    409965               DO jk=1,jpkm1               
    410966                  vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk)            
    … …  
    413969            END DO 
    414970         END DO 
     971         ! 
     972         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     973            vb_b(i1:i2,j1:j2)  = vn_b(i1:i2,j1:j2) 
     974         ENDIF 
     975         ! 
    415976      ENDIF 
    416977      !  
    … …  
    436997         END DO 
    437998      ELSE 
    438          IF( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. .NOT.ln_bt_fw ) ) THEN 
    439             IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    440                DO jj=j1,j2 
    441                   DO ji=i1,i2 
    442                      sshb(ji,jj) =   sshb(ji,jj) & 
    443                            & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 
    444                   END DO 
    445                END DO 
    446             ENDIF 
     999         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     1000            DO jj=j1,j2 
     1001               DO ji=i1,i2 
     1002                  sshb(ji,jj) =   sshb(ji,jj) & 
     1003                        & + atfp * ( tabres(ji,jj) - sshn(ji,jj) ) * tmask(ji,jj,1) 
     1004               END DO 
     1005            END DO 
    4471006         ENDIF 
    4481007         ! 
    … …  
    4521011            END DO 
    4531012         END DO 
     1013         ! 
     1014         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     1015            sshb(i1:i2,j1:j2)  = sshn(i1:i2,j1:j2) 
     1016         ENDIF 
     1017         ! 
     1018 
    4541019      ENDIF 
    4551020      ! 
    … …  
    4651030      LOGICAL                            , INTENT(in) ::   before 
    4661031      !! 
    467       INTEGER ::   ji, jj 
    468       REAL(wp)::   zrhoy 
    469       !!---------------------------------------------------------------------- 
     1032      INTEGER :: ji, jj 
     1033      REAL(wp) :: zrhoy, za1, zcor 
     1034      !!--------------------------------------------- 
    4701035      ! 
    4711036      IF (before) THEN 
    … …  
    4781043         tabres = zrhoy * tabres 
    4791044      ELSE 
     1045         !  
     1046         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2) 
     1047         ! 
     1048         za1 = 1._wp / REAL(Agrif_rhot(), wp) 
    4801049         DO jj=j1,j2 
    4811050            DO ji=i1,i2 
    482                ub2_b(ji,jj) = tabres(ji,jj) / e2u(ji,jj) 
     1051               zcor=tabres(ji,jj) - ub2_b(ji,jj) 
     1052               ! Update time integrated fluxes also in case of multiply nested grids: 
     1053               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor  
     1054               ! Update corrective fluxes: 
     1055               un_bf(ji,jj)  = un_bf(ji,jj) + zcor 
     1056               ! Update half step back fluxes: 
     1057               ub2_b(ji,jj) = tabres(ji,jj) 
    4831058            END DO 
    4841059         END DO 
    … …  
    4871062   END SUBROUTINE updateub2b 
    4881063 
     1064   SUBROUTINE reflux_sshu( tabres, i1, i2, j1, j2, before, nb, ndir ) 
     1065      !!--------------------------------------------- 
     1066      !!          *** ROUTINE reflux_sshu *** 
     1067      !!--------------------------------------------- 
     1068      INTEGER, INTENT(in) :: i1, i2, j1, j2 
     1069      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
     1070      LOGICAL, INTENT(in) :: before 
     1071      INTEGER, INTENT(in) :: nb, ndir 
     1072      !! 
     1073      LOGICAL :: western_side, eastern_side  
     1074      INTEGER :: ji, jj 
     1075      REAL(wp) :: zrhoy, za1, zcor 
     1076      !!--------------------------------------------- 
     1077      ! 
     1078      IF (before) THEN 
     1079         zrhoy = Agrif_Rhoy() 
     1080         DO jj=j1,j2 
     1081            DO ji=i1,i2 
     1082               tabres(ji,jj) = ub2_i_b(ji,jj) * e2u(ji,jj) 
     1083            END DO 
     1084         END DO 
     1085         tabres = zrhoy * tabres 
     1086      ELSE 
     1087         !  
     1088         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e2u(i1:i2,j1:j2) 
     1089         ! 
     1090         western_side  = (nb == 1).AND.(ndir == 1) 
     1091         eastern_side  = (nb == 1).AND.(ndir == 2) 
     1092         ! 
     1093         IF (western_side) THEN 
     1094            DO jj=j1,j2 
     1095               zcor = rdt * r1_e1e2t(i1  ,jj) * e2u(i1,jj) * (ub2_b(i1,jj)-tabres(i1,jj))  
     1096               sshn(i1  ,jj) = sshn(i1  ,jj) + zcor 
     1097               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i1  ,jj) = sshb(i1  ,jj) + atfp * zcor 
     1098            END DO 
     1099         ENDIF 
     1100         IF (eastern_side) THEN 
     1101            DO jj=j1,j2 
     1102               zcor = - rdt * r1_e1e2t(i2+1,jj) * e2u(i2,jj) * (ub2_b(i2,jj)-tabres(i2,jj)) 
     1103               sshn(i2+1,jj) = sshn(i2+1,jj) + zcor 
     1104               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(i2+1,jj) = sshb(i2+1,jj) + atfp * zcor 
     1105            END DO 
     1106         ENDIF 
     1107         ! 
     1108      ENDIF 
     1109      ! 
     1110   END SUBROUTINE reflux_sshu 
    4891111 
    4901112   SUBROUTINE updatevb2b( tabres, i1, i2, j1, j2, before ) 
    … …  
    4961118      LOGICAL                         , INTENT(in   ) ::   before 
    4971119      !! 
    498       INTEGER ::   ji, jj 
    499       REAL(wp)::   zrhox 
    500       !!---------------------------------------------------------------------- 
     1120      INTEGER :: ji, jj 
     1121      REAL(wp) :: zrhox, za1, zcor 
     1122      !!--------------------------------------------- 
    5011123      ! 
    5021124      IF( before ) THEN 
    … …  
    5091131         tabres = zrhox * tabres 
    5101132      ELSE 
     1133         !  
     1134         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2) 
     1135         ! 
     1136         za1 = 1._wp / REAL(Agrif_rhot(), wp) 
    5111137         DO jj=j1,j2 
    5121138            DO ji=i1,i2 
    513                vb2_b(ji,jj) = tabres(ji,jj) / e1v(ji,jj) 
     1139               zcor=tabres(ji,jj) - vb2_b(ji,jj) 
     1140               ! Update time integrated fluxes also in case of multiply nested grids: 
     1141               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor  
     1142               ! Update corrective fluxes: 
     1143               vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor 
     1144               ! Update half step back fluxes: 
     1145               vb2_b(ji,jj) = tabres(ji,jj) 
    5141146            END DO 
    5151147         END DO 
    … …  
    5181150   END SUBROUTINE updatevb2b 
    5191151 
     1152   SUBROUTINE reflux_sshv( tabres, i1, i2, j1, j2, before, nb, ndir ) 
     1153      !!--------------------------------------------- 
     1154      !!          *** ROUTINE reflux_sshv *** 
     1155      !!--------------------------------------------- 
     1156      INTEGER, INTENT(in) :: i1, i2, j1, j2 
     1157      REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres 
     1158      LOGICAL, INTENT(in) :: before 
     1159      INTEGER, INTENT(in) :: nb, ndir 
     1160      !! 
     1161      LOGICAL :: southern_side, northern_side  
     1162      INTEGER :: ji, jj 
     1163      REAL(wp) :: zrhox, za1, zcor 
     1164      !!--------------------------------------------- 
     1165      ! 
     1166      IF (before) THEN 
     1167         zrhox = Agrif_Rhox() 
     1168         DO jj=j1,j2 
     1169            DO ji=i1,i2 
     1170               tabres(ji,jj) = vb2_i_b(ji,jj) * e1v(ji,jj)  
     1171            END DO 
     1172         END DO 
     1173         tabres = zrhox * tabres 
     1174      ELSE 
     1175         !  
     1176         tabres(i1:i2,j1:j2) = tabres(i1:i2,j1:j2) * r1_e1v(i1:i2,j1:j2) 
     1177         ! 
     1178         southern_side = (nb == 2).AND.(ndir == 1) 
     1179         northern_side = (nb == 2).AND.(ndir == 2) 
     1180         ! 
     1181         IF (southern_side) THEN 
     1182            DO ji=i1,i2 
     1183               zcor = rdt * r1_e1e2t(ji,j1  ) * e1v(ji,j1  ) * (vb2_b(ji,j1)-tabres(ji,j1)) 
     1184               sshn(ji,j1  ) = sshn(ji,j1  ) + zcor 
     1185               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j1  ) = sshb(ji,j1) + atfp * zcor 
     1186            END DO 
     1187         ENDIF 
     1188         IF (northern_side) THEN                
     1189            DO ji=i1,i2 
     1190               zcor = - rdt * r1_e1e2t(ji,j2+1) * e1v(ji,j2  ) * (vb2_b(ji,j2)-tabres(ji,j2)) 
     1191               sshn(ji,j2+1) = sshn(ji,j2+1) + zcor 
     1192               IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) sshb(ji,j2+1) = sshb(ji,j2+1) + atfp * zcor 
     1193            END DO 
     1194         ENDIF 
     1195         !  
     1196      ENDIF 
     1197      ! 
     1198   END SUBROUTINE reflux_sshv 
    5201199 
    5211200   SUBROUTINE update_scales( tabres, i1, i2, j1, j2, k1, k2, n1,n2, before ) 
    … …  
    6191298   END SUBROUTINE updateAVM 
    6201299 
     1300   SUBROUTINE updatee3t(ptab_dum, i1, i2, j1, j2, k1, k2, before ) 
     1301      !!--------------------------------------------- 
     1302      !!           *** ROUTINE updatee3t *** 
     1303      !!--------------------------------------------- 
     1304      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ptab_dum 
     1305      INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2 
     1306      LOGICAL, INTENT(in) :: before 
     1307      ! 
     1308      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ptab 
     1309      INTEGER :: ji,jj,jk 
     1310      REAL(wp) :: zcoef 
     1311      !!--------------------------------------------- 
     1312      ! 
     1313      IF (.NOT.before) THEN 
     1314         ! 
     1315         ALLOCATE(ptab(i1:i2,j1:j2,1:jpk))  
     1316         ! 
     1317         ! Update e3t from ssh (z* case only) 
     1318         DO jk = 1, jpkm1 
     1319            DO jj=j1,j2 
     1320               DO ji=i1,i2 
     1321                  ptab(ji,jj,jk) = e3t_0(ji,jj,jk) * (1._wp + sshn(ji,jj) & 
     1322                                     & *ssmask(ji,jj)/(ht_0(ji,jj)-1._wp + ssmask(ji,jj))) 
     1323               END DO 
     1324            END DO 
     1325         END DO 
     1326         ! 
     1327         ! 1) Updates at BEFORE time step: 
     1328         ! ------------------------------- 
     1329         ! 
     1330         ! Save "old" scale factor (prior update) for subsequent asselin correction 
     1331         ! of prognostic variables 
     1332         e3t_a(i1:i2,j1:j2,1:jpkm1) = e3t_n(i1:i2,j1:j2,1:jpkm1) 
     1333 
     1334         ! One should also save e3t_b, but lacking of workspace... 
     1335!         hdivn(i1:i2,j1:j2,1:jpkm1)   = e3t_b(i1:i2,j1:j2,1:jpkm1) 
     1336 
     1337         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0) )) THEN 
     1338            DO jk = 1, jpkm1 
     1339               DO jj=j1,j2 
     1340                  DO ji=i1,i2 
     1341                     e3t_b(ji,jj,jk) =  e3t_b(ji,jj,jk) & 
     1342                           & + atfp * ( ptab(ji,jj,jk) - e3t_n(ji,jj,jk) ) 
     1343                  END DO 
     1344               END DO 
     1345            END DO 
     1346            ! 
     1347            e3w_b  (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_b(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1) 
     1348            gdepw_b(i1:i2,j1:j2,1) = 0.0_wp 
     1349            gdept_b(i1:i2,j1:j2,1) = 0.5_wp * e3w_b(i1:i2,j1:j2,1) 
     1350            ! 
     1351            DO jk = 2, jpk 
     1352               DO jj = j1,j2 
     1353                  DO ji = i1,i2             
     1354                     zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     1355                     e3w_b(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) *        &  
     1356                     &                                        ( e3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )  & 
     1357                     &                                  +            0.5_wp * tmask(ji,jj,jk)   *        & 
     1358                     &                                        ( e3t_b(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) ) 
     1359                     gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 
     1360                     gdept_b(ji,jj,jk) =      zcoef  * ( gdepw_b(ji,jj,jk  ) + 0.5 * e3w_b(ji,jj,jk))  & 
     1361                         &               + (1-zcoef) * ( gdept_b(ji,jj,jk-1) +       e3w_b(ji,jj,jk))  
     1362                  END DO 
     1363               END DO 
     1364            END DO 
     1365            ! 
     1366         ENDIF         
     1367         ! 
     1368         ! 2) Updates at NOW time step: 
     1369         ! ---------------------------- 
     1370         ! 
     1371         ! Update vertical scale factor at T-points: 
     1372         e3t_n(i1:i2,j1:j2,1:jpkm1) = ptab(i1:i2,j1:j2,1:jpkm1) 
     1373         ! 
     1374         ! Update total depth: 
     1375         ht_n(i1:i2,j1:j2) = 0._wp 
     1376         DO jk = 1, jpkm1 
     1377            ht_n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t_n(i1:i2,j1:j2,jk) * tmask(i1:i2,j1:j2,jk) 
     1378         END DO 
     1379         ! 
     1380         ! Update vertical scale factor at W-points and depths: 
     1381         e3w_n (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + e3t_n(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1) 
     1382         gdept_n(i1:i2,j1:j2,1) = 0.5_wp * e3w_n(i1:i2,j1:j2,1) 
     1383         gdepw_n(i1:i2,j1:j2,1) = 0.0_wp 
     1384         gde3w_n(i1:i2,j1:j2,1) = gdept_n(i1:i2,j1:j2,1) - (ht_n(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh 
     1385         ! 
     1386         DO jk = 2, jpk 
     1387            DO jj = j1,j2 
     1388               DO ji = i1,i2             
     1389               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 
     1390               e3w_n(ji,jj,jk)  = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( e3t_n(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) )   & 
     1391               &                                  +            0.5_wp * tmask(ji,jj,jk)   * ( e3t_n(ji,jj,jk  ) - e3t_0(ji,jj,jk  ) ) 
     1392               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 
     1393               gdept_n(ji,jj,jk) =      zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk))  & 
     1394                   &               + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk))  
     1395               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - (ht_n(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh 
     1396               END DO 
     1397            END DO 
     1398         END DO 
     1399         ! 
     1400         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     1401            e3t_b (i1:i2,j1:j2,1:jpk)  = e3t_n (i1:i2,j1:j2,1:jpk) 
     1402            e3w_b (i1:i2,j1:j2,1:jpk)  = e3w_n (i1:i2,j1:j2,1:jpk) 
     1403            gdepw_b(i1:i2,j1:j2,1:jpk) = gdepw_n(i1:i2,j1:j2,1:jpk) 
     1404            gdept_b(i1:i2,j1:j2,1:jpk) = gdept_n(i1:i2,j1:j2,1:jpk) 
     1405         ENDIF 
     1406         ! 
     1407         DEALLOCATE(ptab) 
     1408      ENDIF 
     1409      ! 
     1410   END SUBROUTINE updatee3t 
     1411 
    6211412#else 
    6221413   !!---------------------------------------------------------------------- 

  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_top_interp.F90

    r9019 r9031  
    6464      !!---------------------------------------------------------------------- 
    6565      ! 
    66       IF( before ) THEN          
    67          ptab(i1:i2,j1:j2,k1:k2,n1:n2) = trn(i1:i2,j1:j2,k1:k2,n1:n2) 
     66      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     67      INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out 
     68      REAL(wp) ::   zrhox , zalpha1, zalpha2, zalpha3 
     69      REAL(wp) ::   zalpha4, zalpha5, zalpha6, zalpha7 
     70      LOGICAL :: western_side, eastern_side,northern_side,southern_side 
     71      ! vertical interpolation: 
     72      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 
     73      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 
     74      REAL(wp), DIMENSION(k1:k2) :: h_in 
     75      REAL(wp), DIMENSION(1:jpk) :: h_out(1:jpk) 
     76      REAL(wp) :: h_diff, zrhoxy 
     77 
     78      zrhoxy = Agrif_rhox()*Agrif_rhoy() 
     79      IF (before) THEN          
     80         DO jn = 1,jpts 
     81            DO jk=k1,k2 
     82               DO jj=j1,j2 
     83                 DO ji=i1,i2 
     84                       ptab(ji,jj,jk,jn) = trn(ji,jj,jk,jn) 
     85                 END DO 
     86               END DO 
     87            END DO 
     88         END DO 
     89# if defined key_vertical 
     90        DO jk=k1,k2 
     91           DO jj=j1,j2 
     92              DO ji=i1,i2 
     93                 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)  
     94              END DO 
     95           END DO 
     96        END DO 
     97# endif 
     98 
    6899      ELSE 
    69100         ! 
    70          IF( nbghostcells > 1 ) THEN  ! no smoothing 
    71             tra(i1:i2,j1:j2,k1:k2,n1:n2) = ptab(i1:i2,j1:j2,k1:k2,n1:n2) 
    72          ELSE                         ! smoothing 
    73             ! 
    74             ll_west  = (nb == 1).AND.(ndir == 1)   ;   ll_east  = (nb == 1).AND.(ndir == 2) 
    75             ll_south = (nb == 2).AND.(ndir == 1)   ;   ll_north = (nb == 2).AND.(ndir == 2) 
    76             ! 
    77             zrhox = Agrif_Rhox() 
    78             z1 = ( zrhox - 1. ) * 0.5 
    79             z3 = ( zrhox - 1. ) / ( zrhox + 1. ) 
    80             z6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. ) 
    81             z7 =    - ( zrhox - 1. ) / ( zrhox + 3. ) 
    82             ! 
    83             z2 = 1. - z1 
    84             z4 = 1. - z3 
    85             z5 = 1. - z6 - z7 
    86             ! 
    87             imin = i1   ;   imax = i2 
    88             jmin = j1   ;   jmax = j2 
    89             !  
    90             ! Remove CORNERS 
    91             IF((nbondj == -1).OR.(nbondj == 2))   jmin = 3 
    92             IF((nbondj == +1).OR.(nbondj == 2))   jmax = nlcj-2 
    93             IF((nbondi == -1).OR.(nbondi == 2))   imin = 3 
    94             IF((nbondi == +1).OR.(nbondi == 2))   imax = nlci-2         
    95             ! 
    96             IF( ll_east ) THEN       !== eastern side  ==! 
    97                DO jn = 1, jptra 
    98                   tra(nlci,j1:j2,k1:k2,jn) = z1 * ptab(nlci,j1:j2,k1:k2,jn) + z2 * ptab(nlci-1,j1:j2,k1:k2,jn) 
    99                   DO jk = 1, jpkm1 
    100                      DO jj = jmin,jmax 
    101                         IF( umask(nlci-2,jj,jk) == 0._wp ) THEN 
    102                            tra(nlci-1,jj,jk,jn) = tra(nlci,jj,jk,jn) * tmask(nlci-1,jj,jk) 
    103                         ELSE 
    104                            tra(nlci-1,jj,jk,jn) = ( z4*tra(nlci,jj,jk,jn)+z3*tra(nlci-2,jj,jk,jn) ) * tmask(nlci-1,jj,jk) 
    105                            IF( un(nlci-2,jj,jk) > 0._wp ) THEN 
    106                               tra(nlci-1,jj,jk,jn) = ( z6*tra(nlci-2,jj,jk,jn)+z5*tra(nlci,jj,jk,jn)   &  
    107                                  &                    +z7*tra(nlci-3,jj,jk,jn) ) * tmask(nlci-1,jj,jk) 
    108                            ENDIF 
    109                         ENDIF 
    110                      END DO 
    111                   END DO 
     101         western_side  = (nb == 1).AND.(ndir == 1) 
     102         eastern_side  = (nb == 1).AND.(ndir == 2) 
     103         southern_side = (nb == 2).AND.(ndir == 1) 
     104         northern_side = (nb == 2).AND.(ndir == 2) 
     105 
     106# if defined key_vertical               
     107         DO jj=j1,j2 
     108            DO ji=i1,i2 
     109               iref = ji 
     110               jref = jj 
     111               if(western_side) iref=MAX(2,ji) 
     112               if(eastern_side) iref=MIN(nlci-1,ji) 
     113               if(southern_side) jref=MAX(2,jj) 
     114               if(northern_side) jref=MIN(nlcj-1,jj) 
     115               N_in = 0 
     116               DO jk=k1,k2 !k2 = jpk of parent grid 
     117                  IF (ptab(ji,jj,jk,n2) == 0) EXIT 
     118                  N_in = N_in + 1 
     119                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 
     120                  h_in(N_in) = ptab(ji,jj,jk,n2) 
     121               END DO 
     122               N_out = 0 
     123               DO jk=1,jpk ! jpk of child grid 
     124                  IF (tmask(iref,jref,jk) == 0) EXIT  
     125                  N_out = N_out + 1 
     126                  h_out(jk) = e3t_n(iref,jref,jk) 
    112127               ENDDO 
    113             ENDIF 
    114             !  
    115             IF( ll_north ) THEN        !==  northern side  ==! 
    116                DO jn = 1, jptra 
    117                   tra(i1:i2,nlcj,k1:k2,jn) = z1 * ptab(i1:i2,nlcj,k1:k2,jn) + z2 * ptab(i1:i2,nlcj-1,k1:k2,jn) 
    118                   DO jk = 1, jpkm1 
    119                      DO ji = imin, imax 
    120                         IF( vmask(ji,nlcj-2,jk) == 0._wp ) THEN 
    121                            tra(ji,nlcj-1,jk,jn) = tra(ji,nlcj,jk,jn) * tmask(ji,nlcj-1,jk) 
    122                         ELSE 
    123                            tra(ji,nlcj-1,jk,jn) = ( z4*tra(ji,nlcj,jk,jn)+z3*tra(ji,nlcj-2,jk,jn) ) * tmask(ji,nlcj-1,jk)         
    124                            IF (vn(ji,nlcj-2,jk) > 0._wp ) THEN 
    125                               tra(ji,nlcj-1,jk,jn) = ( z6*tra(ji,nlcj-2,jk,jn)+z5*tra(ji,nlcj,jk,jn)  & 
    126                                  &                    +z7*tra(ji,nlcj-3,jk,jn) ) * tmask(ji,nlcj-1,jk) 
    127                            ENDIF 
    128                         ENDIF 
    129                      END DO 
    130                   END DO 
    131                END DO 
    132             ENDIF 
    133             ! 
    134             IF( ll_west ) THEN         !==  western side  ==!           
    135                DO jn = 1, jptra 
    136                   tra(1,j1:j2,k1:k2,jn) = z1 * ptab(1,j1:j2,k1:k2,jn) + z2 * ptab(2,j1:j2,k1:k2,jn) 
    137                   DO jk = 1, jpkm1 
    138                      DO jj = jmin,jmax 
    139                         IF( umask(2,jj,jk) == 0._wp ) THEN 
    140                            tra(2,jj,jk,jn) = tra(1,jj,jk,jn) * tmask(2,jj,jk) 
    141                         ELSE 
    142                            tra(2,jj,jk,jn) = ( z4*tra(1,jj,jk,jn)+z3*tra(3,jj,jk,jn) ) * tmask(2,jj,jk)         
    143                            IF( un(2,jj,jk) < 0._wp ) THEN 
    144                               tra(2,jj,jk,jn) = ( z6*tra(3,jj,jk,jn)+z5*tra(1,jj,jk,jn)+z7*tra(4,jj,jk,jn) ) * tmask(2,jj,jk) 
    145                            ENDIF 
    146                         ENDIF 
    147                      END DO 
    148                   END DO 
    149                END DO 
    150             ENDIF 
    151             ! 
    152             IF( ll_south ) THEN        !==  southern side  ==! 
    153                DO jn = 1, jptra 
    154                   tra(i1:i2,1,k1:k2,jn) = z1 * ptab(i1:i2,1,k1:k2,jn) + z2 * ptab(i1:i2,2,k1:k2,jn) 
    155                   DO jk = 1, jpk       
    156                      DO ji = imin, imax 
    157                         IF( vmask(ji,2,jk) == 0._wp ) THEN 
    158                            tra(ji,2,jk,jn) = tra(ji,1,jk,jn) * tmask(ji,2,jk) 
    159                         ELSE 
    160                            tra(ji,2,jk,jn) = ( z4*tra(ji,1,jk,jn)+z3*tra(ji,3,jk,jn) ) * tmask(ji,2,jk) 
    161                            IF( vn(ji,2,jk) < 0._wp ) THEN 
    162                               tra(ji,2,jk,jn) = ( z6*tra(ji,3,jk,jn)+z5*tra(ji,1,jk,jn)+z7*tra(ji,4,jk,jn) ) * tmask(ji,2,jk) 
    163                            ENDIF 
     128               IF (N_in > 0) THEN 
     129                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     130                  DO jn=1,jptra 
     131                     call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
     132                  ENDDO 
     133               ENDIF 
     134            ENDDO 
     135         ENDDO 
     136# else 
     137         ptab_child(i1:i2,j1:j2,1:jpk,1:jptra) = ptab(i1:i2,j1:j2,1:jpk,1:jptra) 
     138# endif 
     139 
     140         ! 
     141         zrhox = Agrif_Rhox() 
     142         !  
     143         zalpha1 = ( zrhox - 1. ) * 0.5 
     144         zalpha2 = 1. - zalpha1 
     145         !  
     146         zalpha3 = ( zrhox - 1. ) / ( zrhox + 1. ) 
     147         zalpha4 = 1. - zalpha3 
     148         !  
     149         zalpha6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. ) 
     150         zalpha7 =    - ( zrhox - 1. ) / ( zrhox + 3. ) 
     151         zalpha5 = 1. - zalpha6 - zalpha7 
     152         ! 
     153         imin = i1 
     154         imax = i2 
     155         jmin = j1 
     156         jmax = j2 
     157         !  
     158         ! Remove CORNERS 
     159         IF((nbondj == -1).OR.(nbondj == 2)) jmin = 3 
     160         IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj-2 
     161         IF((nbondi == -1).OR.(nbondi == 2)) imin = 3 
     162         IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci-2         
     163         ! 
     164         IF( eastern_side) THEN 
     165            DO jn = 1, jptra 
     166               tra(nlci,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(nlci,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(nlci-1,j1:j2,1:jpk,jn) 
     167               DO jk = 1, jpkm1 
     168                  DO jj = jmin,jmax 
     169                     IF( umask(nlci-2,jj,jk) == 0.e0 ) THEN 
     170                        tra(nlci-1,jj,jk,jn) = tra(nlci,jj,jk,jn) * tmask(nlci-1,jj,jk) 
     171                     ELSE 
     172                        tra(nlci-1,jj,jk,jn)=(zalpha4*tra(nlci,jj,jk,jn)+zalpha3*tra(nlci-2,jj,jk,jn))*tmask(nlci-1,jj,jk) 
     173                        IF( un(nlci-2,jj,jk) > 0.e0 ) THEN 
     174                           tra(nlci-1,jj,jk,jn)=( zalpha6*tra(nlci-2,jj,jk,jn)+zalpha5*tra(nlci,jj,jk,jn) &  
     175                                 + zalpha7*tra(nlci-3,jj,jk,jn) ) * tmask(nlci-1,jj,jk) 
     176                        ENDIF 
     177                     END DO 
     178                  END DO 
     179               END DO 
     180            ENDDO 
     181         ENDIF 
     182         !  
     183         IF( northern_side ) THEN             
     184            DO jn = 1, jptra 
     185               tra(i1:i2,nlcj,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,nlcj,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,nlcj-1,1:jpk,jn) 
     186               DO jk = 1, jpkm1 
     187                  DO ji = imin,imax 
     188                     IF( vmask(ji,nlcj-2,jk) == 0.e0 ) THEN 
     189                        tra(ji,nlcj-1,jk,jn) = tra(ji,nlcj,jk,jn) * tmask(ji,nlcj-1,jk) 
     190                     ELSE 
     191                        tra(ji,nlcj-1,jk,jn)=(zalpha4*tra(ji,nlcj,jk,jn)+zalpha3*tra(ji,nlcj-2,jk,jn))*tmask(ji,nlcj-1,jk)         
     192                        IF (vn(ji,nlcj-2,jk) > 0.e0 ) THEN 
     193                           tra(ji,nlcj-1,jk,jn)=( zalpha6*tra(ji,nlcj-2,jk,jn)+zalpha5*tra(ji,nlcj,jk,jn)  & 
     194                                 + zalpha7*tra(ji,nlcj-3,jk,jn) ) * tmask(ji,nlcj-1,jk) 
     195                        ENDIF 
     196                     END DO 
     197                  END DO 
     198               END DO 
     199            ENDDO 
     200         ENDIF 
     201         ! 
     202         IF( western_side) THEN             
     203            DO jn = 1, jptra 
     204               tra(1,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(1,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(2,j1:j2,1:jpk,jn) 
     205               DO jk = 1, jpkm1 
     206                  DO jj = jmin,jmax 
     207                     IF( umask(2,jj,jk) == 0.e0 ) THEN 
     208                        tra(2,jj,jk,jn) = tra(1,jj,jk,jn) * tmask(2,jj,jk) 
     209                     ELSE 
     210                        tra(2,jj,jk,jn)=(zalpha4*tra(1,jj,jk,jn)+zalpha3*tra(3,jj,jk,jn))*tmask(2,jj,jk)         
     211                        IF( un(2,jj,jk) < 0.e0 ) THEN 
     212                           tra(2,jj,jk,jn)=(zalpha6*tra(3,jj,jk,jn)+zalpha5*tra(1,jj,jk,jn)+zalpha7*tra(4,jj,jk,jn))*tmask(2,jj,jk) 
     213                        ENDIF 
     214                     END DO 
     215                  END DO 
     216               END DO 
     217            END DO 
     218         ENDIF 
     219         ! 
     220         IF( southern_side ) THEN            
     221            DO jn = 1, jptra 
     222               tra(i1:i2,1,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,1,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,2,1:jpk,jn) 
     223               DO jk=1,jpkm1 
     224                  DO ji=imin,imax 
     225                     IF( vmask(ji,2,jk) == 0.e0 ) THEN 
     226                        tra(ji,2,jk,jn)=tra(ji,1,jk,jn) * tmask(ji,2,jk) 
     227                     ELSE 
     228                        tra(ji,2,jk,jn)=(zalpha4*tra(ji,1,jk,jn)+zalpha3*tra(ji,3,jk,jn))*tmask(ji,2,jk) 
     229                        IF( vn(ji,2,jk) < 0.e0 ) THEN 
     230                           tra(ji,2,jk,jn)=(zalpha6*tra(ji,3,jk,jn)+zalpha5*tra(ji,1,jk,jn)+zalpha7*tra(ji,4,jk,jn))*tmask(ji,2,jk) 
    164231                        ENDIF 
    165232                     END DO 
    … …  
    175242            ! 
    176243         ENDIF 
     244         ! 
     245         ! Treatment of corners 
     246         !  
     247         ! East south 
     248         IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 
     249            tra(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,:) 
     250         ENDIF 
     251         ! East north 
     252         IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 
     253            tra(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,:) 
     254         ENDIF 
     255         ! West south 
     256         IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 
     257            tra(2,2,:,:) = ptab_child(2,2,:,:) 
     258         ENDIF 
     259         ! West north 
     260         IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 
     261            tra(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,:) 
     262         ENDIF 
     263         ! 
    177264      ENDIF 
    178265      ! 

  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_top_sponge.F90

    r9019 r9031  
    4444      ! 
    4545#if defined SPONGE_TOP 
    46       zcoef = REAL( Agrif_NbStepint(), wp ) / Agrif_rhot() 
     46!! Assume persistence  
     47      timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot()) 
    4748      CALL Agrif_sponge 
    4849      Agrif_SpecialValue    = 0._wp 
    … …  
    6869      REAL(wp), DIMENSION(i1:i2,j1:j2)             ::   ztu, ztv 
    6970      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::   trbdiff 
     71      ! vertical interpolation: 
     72      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child 
     73      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 
     74      REAL(wp), DIMENSION(k1:k2) :: h_in 
     75      REAL(wp), DIMENSION(1:jpk) :: h_out 
     76      INTEGER :: N_in, N_out 
     77      REAL(wp) :: h_diff 
    7078      !!---------------------------------------------------------------------- 
    7179      ! 
    7280      IF( before ) THEN 
    73          tabres(i1:i2,j1:j2,k1:k2,n1:n2) = trn(i1:i2,j1:j2,k1:k2,n1:n2) 
     81         DO jn = 1, jptra 
     82            DO jk=k1,k2 
     83               DO jj=j1,j2 
     84                  DO ji=i1,i2 
     85                     tabres(ji,jj,jk,jn) = trb(ji,jj,jk,jn) 
     86                  END DO 
     87               END DO 
     88            END DO 
     89         END DO 
     90 
     91# if defined key_vertical 
     92         DO jk=k1,k2 
     93            DO jj=j1,j2 
     94               DO ji=i1,i2 
     95                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)  
     96               END DO 
     97            END DO 
     98         END DO 
     99# endif 
    74100      ELSE       
    75 !!gm line below use of :,:  versus i1:i2,j1:j2  ....   strange, not wrong.    ===>> to be corrected 
    76          trbdiff(:,:,:,:) = trb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:)       
     101# if defined key_vertical 
     102         tabres_child(:,:,:,:) = 0. 
     103         DO jj=j1,j2 
     104            DO ji=i1,i2 
     105               N_in = 0 
     106               DO jk=k1,k2 !k2 = jpk of parent grid 
     107                  IF (tabres(ji,jj,jk,n2) == 0) EXIT 
     108                  N_in = N_in + 1 
     109                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 
     110                  h_in(N_in) = tabres(ji,jj,jk,n2) 
     111               END DO 
     112               N_out = 0 
     113               DO jk=1,jpk ! jpk of child grid 
     114                  IF (tmask(ji,jj,jk) == 0) EXIT  
     115                  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 
     117               ENDDO 
     118               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????? 
     121                  DO jn=1,jptra 
     122                     call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 
     123                  ENDDO 
     124               ENDIF 
     125            ENDDO 
     126         ENDDO 
     127# endif 
     128 
     129         DO jj=j1,j2 
     130            DO ji=i1,i2 
     131               DO jk=1,jpkm1 
     132# if defined key_vertical 
     133                  trbdiff(ji,jj,jk,1:jptra) = trb(ji,jj,jk,1:jptra) - tabres_child(ji,jj,jk,1:jptra) 
     134# else 
     135                  trbdiff(ji,jj,jk,1:jptra) = trb(ji,jj,jk,1:jptra) - tabres(ji,jj,jk,1:jptra) 
     136# endif 
     137               ENDDO 
     138            ENDDO 
     139         ENDDO 
     140 
    77141         DO jn = 1, jptra 
    78142            DO jk = 1, jpkm1 

  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90

    r9019 r9031  
    2626   PUBLIC Agrif_Update_Trc 
    2727 
    28    INTEGER, PUBLIC ::   nbcline_trc = 0   !: ??? 
    29  
    3028   !!---------------------------------------------------------------------- 
    3129   !! NEMO/NST 4.0 , NEMO Consortium (2017) 
    … …  
    3533CONTAINS 
    3634 
    37    SUBROUTINE Agrif_Update_Trc( kt ) 
     35   SUBROUTINE Agrif_Update_Trc( ) 
    3836      !!---------------------------------------------------------------------- 
    3937      !!                   *** ROUTINE Agrif_Update_Trc *** 
    4038      !!---------------------------------------------------------------------- 
    41       INTEGER, INTENT(in) ::   kt 
    42       !!---------------------------------------------------------------------- 
    43       !  
    44       IF((Agrif_NbStepint() .NE. (Agrif_irhot()-1)).AND.(kt /= 0)) RETURN 
     39      !  
     40      IF (Agrif_Root()) RETURN  
     41      ! 
    4542#if defined TWO_WAY    
    4643      Agrif_UseSpecialValueInUpdate = .TRUE. 
    4744      Agrif_SpecialValueFineGrid    = 0._wp 
    4845      !  
    49       IF( MOD(nbcline_trc,nbclineupdate) == 0 ) THEN 
    5046# if ! defined DECAL_FEEDBACK 
    51          CALL Agrif_Update_Variable(trn_id, procname=updateTRC ) 
     47      CALL Agrif_Update_Variable(trn_id, procname=updateTRC ) 
     48!      CALL Agrif_Update_Variable( trn_id, locupdate=(/0,2/), procname=updateTRC ) 
    5249# else 
    53          CALL Agrif_Update_Variable(trn_id, locupdate=(/1,0/),procname=updateTRC ) 
     50      CALL Agrif_Update_Variable(trn_id, locupdate=(/1,0/),procname=updateTRC ) 
     51!      CALL Agrif_Update_Variable( trn_id, locupdate=(/1,2/), procname=updateTRC ) 
    5452# endif 
     53      ! 
     54      Agrif_UseSpecialValueInUpdate = .FALSE. 
     55      ! 
     56#endif 
     57      ! 
     58   END SUBROUTINE Agrif_Update_Trc 
     59 
     60#ifdef key_vertical 
     61   SUBROUTINE updateTRC( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     62      !!--------------------------------------------- 
     63      !!           *** ROUTINE updateT *** 
     64      !!--------------------------------------------- 
     65      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
     66      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
     67      LOGICAL, INTENT(in) :: before 
     68      !! 
     69      INTEGER :: ji,jj,jk,jn 
     70      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 
     71      REAL(wp) :: h_in(k1:k2) 
     72      REAL(wp) :: h_out(1:jpk) 
     73      INTEGER  :: N_in, N_out 
     74      REAL(wp) :: h_diff 
     75      REAL(wp) :: zrho_xy 
     76      REAL(wp) :: tabin(k1:k2,n1:n2) 
     77      !!--------------------------------------------- 
     78      ! 
     79      IF (before) THEN 
     80         AGRIF_SpecialValue = -999._wp 
     81         zrho_xy = Agrif_rhox() * Agrif_rhoy()  
     82         DO jn = n1,n2-1 
     83            DO jk=k1,k2 
     84               DO jj=j1,j2 
     85                  DO ji=i1,i2 
     86                     tabres(ji,jj,jk,jn) = (trn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) & 
     87                                           * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp 
     88                  END DO 
     89               END DO 
     90            END DO 
     91         END DO 
     92         DO jk=k1,k2 
     93            DO jj=j1,j2 
     94               DO ji=i1,i2 
     95                  tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) & 
     96                                           + (tmask(ji,jj,jk)-1)*999._wp 
     97               END DO 
     98            END DO 
     99         END DO 
    55100      ELSE 
    56 # if ! defined DECAL_FEEDBACK 
    57          CALL Agrif_Update_Variable( trn_id, locupdate=(/0,2/), procname=updateTRC ) 
    58 # else 
    59          CALL Agrif_Update_Variable( trn_id, locupdate=(/1,2/), procname=updateTRC ) 
    60 # endif 
     101         tabres_child(:,:,:,:) = 0. 
     102         AGRIF_SpecialValue = 0._wp 
     103         DO jj=j1,j2 
     104            DO ji=i1,i2 
     105               N_in = 0 
     106               DO jk=k1,k2 !k2 = jpk of child grid 
     107                  IF (tabres(ji,jj,jk,n2) == 0  ) EXIT 
     108                  N_in = N_in + 1 
     109                  tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     110                  h_in(N_in) = tabres(ji,jj,jk,n2) 
     111               ENDDO 
     112               N_out = 0 
     113               DO jk=1,jpk ! jpk of parent grid 
     114                  IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 
     115                  N_out = N_out + 1 
     116                  h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 
     117               ENDDO 
     118               IF (N_in > 0) THEN !Remove this? 
     119                  h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
     120                  IF (h_diff < -1.e-4) THEN 
     121                     print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 
     122                     print *,h_in(1:N_in) 
     123                     print *,h_out(1:N_out) 
     124                     STOP 
     125                  ENDIF 
     126                  DO jn=1,jptra 
     127                     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) 
     128                  ENDDO 
     129               ENDIF 
     130            ENDDO 
     131         ENDDO 
     132 
     133         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
     134            ! Add asselin part 
     135            DO jn = 1,jptra 
     136               DO jk=1,jpk 
     137                  DO jj=j1,j2 
     138                     DO ji=i1,i2 
     139                        IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 
     140                           trb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) &  
     141                                 & + atfp * ( tabres_child(ji,jj,jk,jn) & 
     142                                 &          - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
     143                        ENDIF 
     144                     ENDDO 
     145                  ENDDO 
     146               ENDDO 
     147            ENDDO 
     148         ENDIF 
     149         DO jn = 1,jptra 
     150            DO jk=1,jpk 
     151               DO jj=j1,j2 
     152                  DO ji=i1,i2 
     153                     IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN  
     154                        trn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 
     155                     END IF 
     156                  END DO 
     157               END DO 
     158            END DO 
     159         END DO 
    61160      ENDIF 
    62       ! 
    63       Agrif_UseSpecialValueInUpdate = .FALSE. 
    64       nbcline_trc = nbcline_trc + 1 
    65 #endif 
    66       ! 
    67    END SUBROUTINE Agrif_Update_Trc 
    68  
    69  
    70    SUBROUTINE updateTRC( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    71       !!---------------------------------------------------------------------- 
    72       !!                      *** ROUTINE updateT *** 
     161      !  
     162   END SUBROUTINE updateTRC 
     163 
     164 
     165#else 
     166   SUBROUTINE updateTRC( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     167      !!---------------------------------------------------------------------- 
     168      !!                      *** ROUTINE updateTRC *** 
    73169      !!---------------------------------------------------------------------- 
    74170      INTEGER                                    , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2 
    75       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab 
     171      REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres 
    76172      LOGICAL                                    , INTENT(in   ) ::   before 
    77173      !! 
    78       INTEGER ::   ji, jj, jk, jn 
    79       !!---------------------------------------------------------------------- 
    80       ! 
    81       IF( before ) THEN 
    82          ptab(i1:i2,j1:j2,k1:k2,n1:n2) = trn(i1:i2,j1:j2,k1:k2,n1:n2) 
     174      INTEGER :: ji,jj,jk,jn 
     175      REAL(wp) :: ztb, ztnu, ztno 
     176      !!---------------------------------------------------------------------- 
     177      ! 
     178      ! 
     179      IF (before) THEN 
     180         DO jn = n1,n2 
     181            DO jk=k1,k2 
     182               DO jj=j1,j2 
     183                  DO ji=i1,i2 
     184!> jc tmp 
     185                     tabres(ji,jj,jk,jn) = trn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) / e3t_0(ji,jj,jk) 
     186!                     tabres(ji,jj,jk,jn) = trn(ji,jj,jk,jn)  * e3t_n(ji,jj,jk) 
     187!< jc tmp 
     188                  END DO 
     189               END DO 
     190            END DO 
     191         END DO 
    83192      ELSE 
    84          IF( .NOT.(lk_agrif_fstep.AND.(neuler==0)) ) THEN 
     193!> jc tmp 
     194         DO jn = n1,n2 
     195            tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 
     196                                         & * tmask(i1:i2,j1:j2,k1:k2) 
     197         ENDDO 
     198!< jc tmp 
     199         IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 
    85200            ! Add asselin part 
    86201            DO jn = n1,n2 
    87                DO jk = k1, k2 
    88                   DO jj = j1, j2 
    89                      DO ji = i1, i2 
    90                         IF( ptab(ji,jj,jk,jn) /= 0._wp ) THEN 
    91                            trb(ji,jj,jk,jn) = trb(ji,jj,jk,jn)             &  
    92                               &             + atfp * ( ptab(ji,jj,jk,jn)   & 
    93                                  &                    - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 
     202               DO jk=k1,k2 
     203                  DO jj=j1,j2 
     204                     DO ji=i1,i2 
     205                        IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 
     206                           ztb  = trb(ji,jj,jk,jn) * e3t_b(ji,jj,jk) ! fse3t_b prior update should be used 
     207                           ztnu = tabres(ji,jj,jk,jn) 
     208                           ztno = trn(ji,jj,jk,jn) * e3t_a(ji,jj,jk) 
     209                           trb(ji,jj,jk,jn) = ( ztb + atfp * ( ztnu - ztno) )  &  
     210                                     &        * tmask(ji,jj,jk) / e3t_b(ji,jj,jk) 
    94211                        ENDIF 
    95                      END DO 
    96                   END DO 
    97                END DO 
    98             END DO 
     212                     ENDDO 
     213                  ENDDO 
     214               ENDDO 
     215            ENDDO 
    99216         ENDIF 
    100          DO jn = n1, n2 
    101             DO jk = k1, k2 
    102                DO jj = j1, j2 
    103                   DO ji = i1, i2 
    104                      IF( ptab(ji,jj,jk,jn) /= 0._wp ) THEN  
    105                         trn(ji,jj,jk,jn) = ptab(ji,jj,jk,jn) * tmask(ji,jj,jk) 
     217         DO jn = n1,n2 
     218            DO jk=k1,k2 
     219               DO jj=j1,j2 
     220                  DO ji=i1,i2 
     221                     IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN  
     222                        trn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / e3t_n(ji,jj,jk) 
    106223                     END IF 
    107224                  END DO 
    … …  
    109226            END DO 
    110227         END DO 
     228         ! 
     229         IF  ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 
     230            trb(i1:i2,j1:j2,k1:k2,n1:n2)  = trn(i1:i2,j1:j2,k1:k2,n1:n2) 
     231         ENDIF 
     232         ! 
    111233      ENDIF 
    112234      !  
    113235   END SUBROUTINE updateTRC 
     236#endif 
    114237 
    115238#else 

  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_user.F90

    r9019 r9031  
     1#define UPD_HIGH   /* MIX HIGH UPDATE */ 
    12#if defined key_agrif 
    23!!---------------------------------------------------------------------- 
    … …  
    8889# endif 
    8990   ! 
     91   IF ( Agrif_Level().EQ.Agrif_MaxLevel() ) CALL agrif_Update_ini() 
     92 
     93   Agrif_UseSpecialValueInUpdate = .FALSE.      
     94 
    9095END SUBROUTINE Agrif_initvalues 
    9196 
    … …  
    149154   CALL Agrif_Set_bc(e2v_id,(/0,ind1-1/)) 
    150155 
    151    ! 5. Update type 
     156   ! 4. Update type 
    152157   !---------------  
    153    CALL Agrif_Set_Updatetype( e1u_id, update1=Agrif_Update_Copy   , update2=Agrif_Update_Average ) 
    154    CALL Agrif_Set_Updatetype( e2v_id, update1=Agrif_Update_Average, update2=Agrif_Update_Copy    ) 
    155  
    156 ! High order updates 
    157 !   CALL Agrif_Set_Updatetype( e1u_id, update1=Agrif_Update_Average       , update2=Agrif_Update_Full_Weighting ) 
    158 !   CALL Agrif_Set_Updatetype( e2v_id, update1=Agrif_Update_Full_Weighting, update2=Agrif_Update_Average        ) 
    159     ! 
     158# if defined UPD_HIGH 
     159   CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Full_Weighting) 
     160   CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Full_Weighting, update2=Agrif_Update_Average) 
     161#else 
     162   CALL Agrif_Set_Updatetype(e1u_id,update1 = Agrif_Update_Copy, update2=Agrif_Update_Average) 
     163   CALL Agrif_Set_Updatetype(e2v_id,update1 = Agrif_Update_Average, update2=Agrif_Update_Copy) 
     164#endif 
     165 
    160166END SUBROUTINE agrif_declare_var_dom 
    161167 
    … …  
    182188   ! 
    183189   LOGICAL :: check_namelist 
    184    CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3 
     190   CHARACTER(len=15) :: cl_check1, cl_check2, cl_check3, cl_check4  
    185191   !!---------------------------------------------------------------------- 
    186192 
    … …  
    212218   Agrif_UseSpecialValue = .TRUE. 
    213219   CALL Agrif_Bc_variable(sshn_id,calledweight=1., procname=interpsshn ) 
     220   hbdy_w(:) = 0.e0 ; hbdy_e(:) = 0.e0 ; hbdy_n(:) = 0.e0 ; hbdy_s(:) = 0.e0 
     221   ssha(:,:) = 0.e0 
    214222 
    215223   IF ( ln_dynspg_ts ) THEN 
    … …  
    219227      CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b) 
    220228      CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b) 
    221       ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 ; hbdy_w(:) =0.e0 
    222       ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 ; hbdy_e(:) =0.e0  
    223       ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 ; hbdy_n(:) =0.e0  
    224       ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 ; hbdy_s(:) =0.e0 
     229      ubdy_w(:) = 0.e0 ; vbdy_w(:) = 0.e0 
     230      ubdy_e(:) = 0.e0 ; vbdy_e(:) = 0.e0 
     231      ubdy_n(:) = 0.e0 ; vbdy_n(:) = 0.e0 
     232      ubdy_s(:) = 0.e0 ; vbdy_s(:) = 0.e0 
    225233   ENDIF 
    226234 
    … …  
    241249         WRITE(cl_check2,*)  NINT(rdt) 
    242250         WRITE(cl_check3,*)  NINT(Agrif_Parent(rdt)/Agrif_Rhot()) 
    243          CALL ctl_stop( 'incompatible time step between ocean grids',   & 
     251         CALL ctl_stop( 'Incompatible time step between ocean grids',   & 
    244252               &               'parent grid value : '//cl_check1    ,   &  
    245253               &               'child  grid value : '//cl_check2    ,   &  
    … …  
    252260         WRITE(cl_check1,*)  (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1 
    253261         WRITE(cl_check2,*)   Agrif_Parent(nitend)   *Agrif_IRhot() 
    254          CALL ctl_warn( 'incompatible run length between grids'               ,   & 
    255                &              ' nit000 on fine grid will be change to : '//cl_check1,   & 
    256                &              ' nitend on fine grid will be change to : '//cl_check2    ) 
     262         CALL ctl_warn( 'Incompatible run length between grids'                      ,   & 
     263               &               'nit000 on fine grid will be changed to : '//cl_check1,   & 
     264               &               'nitend on fine grid will be changed to : '//cl_check2    ) 
    257265         nit000 = (Agrif_Parent(nit000)-1)*Agrif_IRhot() + 1 
    258266         nitend =  Agrif_Parent(nitend)   *Agrif_IRhot() 
    259267      ENDIF 
    260  
    261       ! Check coordinates 
    262      !SF  IF( ln_zps ) THEN 
    263      !SF     ! check parameters for partial steps  
    264      !SF     IF( Agrif_Parent(e3zps_min) .NE. e3zps_min ) THEN 
    265      !SF        WRITE(*,*) 'incompatible e3zps_min between grids' 
    266      !SF        WRITE(*,*) 'parent grid :',Agrif_Parent(e3zps_min) 
    267      !SF        WRITE(*,*) 'child grid  :',e3zps_min 
    268      !SF        WRITE(*,*) 'those values should be identical' 
    269      !SF        STOP 
    270      !SF     ENDIF 
    271      !SF     IF( Agrif_Parent(e3zps_rat) /= e3zps_rat ) THEN 
    272      !SF        WRITE(*,*) 'incompatible e3zps_rat between grids' 
    273      !SF        WRITE(*,*) 'parent grid :',Agrif_Parent(e3zps_rat) 
    274      !SF        WRITE(*,*) 'child grid  :',e3zps_rat 
    275      !SF        WRITE(*,*) 'those values should be identical'                   
    276      !SF        STOP 
    277      !SF     ENDIF 
    278      !SF  ENDIF 
    279268 
    280269      ! Check free surface scheme 
    281270      IF ( ( Agrif_Parent(ln_dynspg_ts ).AND.ln_dynspg_exp ).OR.& 
    282271         & ( Agrif_Parent(ln_dynspg_exp).AND.ln_dynspg_ts ) ) THEN 
    283          WRITE(*,*) 'incompatible free surface scheme between grids' 
    284          WRITE(*,*) 'parent grid ln_dynspg_ts  :', Agrif_Parent(ln_dynspg_ts ) 
    285          WRITE(*,*) 'parent grid ln_dynspg_exp :', Agrif_Parent(ln_dynspg_exp) 
    286          WRITE(*,*) 'child grid  ln_dynspg_ts  :', ln_dynspg_ts 
    287          WRITE(*,*) 'child grid  ln_dynspg_exp :', ln_dynspg_exp 
    288          WRITE(*,*) 'those logicals should be identical'                   
     272         WRITE(cl_check1,*)  Agrif_Parent( ln_dynspg_ts ) 
     273         WRITE(cl_check2,*)  ln_dynspg_ts 
     274         WRITE(cl_check3,*)  Agrif_Parent( ln_dynspg_exp ) 
     275         WRITE(cl_check4,*)  ln_dynspg_exp 
     276         CALL ctl_stop( 'Incompatible free surface scheme between grids' ,  & 
     277               &               'parent grid ln_dynspg_ts  :'//cl_check1  ,  &  
     278               &               'child  grid ln_dynspg_ts  :'//cl_check2  ,  & 
     279               &               'parent grid ln_dynspg_exp :'//cl_check3  ,  & 
     280               &               'child  grid ln_dynspg_exp :'//cl_check4  ,  & 
     281               &               'those logicals should be identical' )                  
     282         STOP 
     283      ENDIF 
     284 
     285      ! Check if identical linear free surface option 
     286      IF ( ( Agrif_Parent(ln_linssh ).AND.(.NOT.ln_linssh )).OR.& 
     287         & ( (.NOT.Agrif_Parent(ln_linssh)).AND.ln_linssh ) ) THEN 
     288         WRITE(cl_check1,*)  Agrif_Parent(ln_linssh ) 
     289         WRITE(cl_check2,*)  ln_linssh 
     290         CALL ctl_stop( 'Incompatible linearized fs option between grids',  & 
     291               &               'parent grid ln_linssh  :'//cl_check1     ,  & 
     292               &               'child  grid ln_linssh  :'//cl_check2     ,  & 
     293               &               'those logicals should be identical' )                   
    289294         STOP 
    290295      ENDIF 
    … …  
    313318   ENDIF 
    314319   !  
    315    ! Do update at initialisation because not done before writing restarts 
    316    ! This would indeed change boundary conditions values at initial time 
    317    ! hence produce restartability issues. 
    318    ! Note that update below is recursive (with lk_agrif_doupd=T): 
    319    !  
    320 ! JC: I am not sure if Agrif_MaxLevel() is the "relative" 
    321 !     or the absolute maximum nesting level...TBC                         
    322    IF ( Agrif_Level().EQ.Agrif_MaxLevel() ) THEN  
    323       ! NB: Do tracers first, dynamics after because nbcline incremented in dynamics 
    324       CALL Agrif_Update_tra() 
    325       CALL Agrif_Update_dyn() 
    326    ENDIF 
    327    ! 
    328    Agrif_UseSpecialValueInUpdate = .FALSE. 
    329    nbcline = 0 
    330    lk_agrif_doupd = .FALSE. 
    331    ! 
    332320END SUBROUTINE Agrif_InitValues_cont 
    333321 
     322RECURSIVE SUBROUTINE Agrif_Update_ini( ) 
     323   !!---------------------------------------------------------------------- 
     324   !!                 *** ROUTINE agrif_Update_ini *** 
     325   !! 
     326   !! ** Purpose :: Recursive update done at initialization 
     327   !!---------------------------------------------------------------------- 
     328   USE dom_oce 
     329   USE agrif_opa_update 
     330#if defined key_top 
     331   USE agrif_top_update 
     332#endif 
     333   ! 
     334   IMPLICIT NONE 
     335   !!---------------------------------------------------------------------- 
     336   ! 
     337   IF (Agrif_Root()) RETURN 
     338   ! 
     339   IF (.NOT.ln_linssh) CALL Agrif_Update_vvl() 
     340   CALL Agrif_Update_tra() 
     341#if defined key_top 
     342   CALL Agrif_Update_Trc() 
     343#endif 
     344   CALL Agrif_Update_dyn() 
     345! JC remove update because this precludes from perfect restartability 
     346!!   CALL Agrif_Update_tke(0) 
     347 
     348   CALL Agrif_ChildGrid_To_ParentGrid() 
     349   CALL Agrif_Update_ini() 
     350   CALL Agrif_ParentGrid_To_ChildGrid() 
     351 
     352END SUBROUTINE agrif_update_ini 
    334353 
    335354SUBROUTINE agrif_declare_var 
    … …  
    355374   ind2 = 1 + nbghostcells 
    356375   ind3 = 2 + nbghostcells 
     376# if defined key_vertical 
     377   CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_id) 
     378   CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_sponge_id) 
     379 
     380   CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id) 
     381   CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_interp_id) 
     382   CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_update_id) 
     383   CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_update_id) 
     384   CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_sponge_id) 
     385   CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_sponge_id) 
     386# else 
    357387   CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_id) 
    358388   CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_sponge_id) 
    359389 
    360    CALL agrif_declare_variable((/1,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_interp_id) 
    361    CALL agrif_declare_variable((/2,1,0/),(/ind3,ind2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_interp_id) 
    362    CALL agrif_declare_variable((/1,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_update_id) 
    363    CALL agrif_declare_variable((/2,1,0/),(/ind3,ind2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_update_id) 
    364    CALL agrif_declare_variable((/1,2,0/),(/ind2,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_sponge_id) 
    365    CALL agrif_declare_variable((/2,1,0/),(/ind3,ind2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_sponge_id) 
     390   CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_interp_id) 
     391   CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_interp_id) 
     392   CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_update_id) 
     393   CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_update_id) 
     394   CALL agrif_declare_variable((/1,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_sponge_id) 
     395   CALL agrif_declare_variable((/2,1,0,0/),(/ind3,ind2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_sponge_id) 
     396# endif 
    366397 
    367398   CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),e3t_id) 
    … …  
    383414      CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/), en_id) 
    384415      CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avt_id) 
    385       CALL agrif_declare_variable((/2,2,0/),(/ind3,ind3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avm_id) 
     416# if defined key_vertical 
     417   CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/jpi,jpj,jpk,2/),avm_id) 
     418# else 
     419   CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),avm_id) 
     420# endif 
    386421   ENDIF 
    387422 
    … …  
    433468   IF( ln_zdftke )   CALL Agrif_Set_bc( avm_id, (/0,ind1/) ) 
    434469 
    435    ! 5. Update type 
     470   ! 4. Update type 
    436471   !---------------  
     472   CALL Agrif_Set_Updatetype(scales_t_id, update = AGRIF_Update_Average) 
     473 
     474# if defined UPD_HIGH 
     475   CALL Agrif_Set_Updatetype(tsn_id, update = Agrif_Update_Full_Weighting) 
     476   CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting) 
     477   CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average) 
     478 
     479   CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting) 
     480   CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average) 
     481   CALL Agrif_Set_Updatetype(sshn_id, update = Agrif_Update_Full_Weighting) 
     482   CALL Agrif_Set_Updatetype(e3t_id, update = Agrif_Update_Full_Weighting) 
     483 
     484   IF( ln_zdftke) THEN 
     485      CALL Agrif_Set_Updatetype( en_id, update = AGRIF_Update_Full_Weighting) 
     486      CALL Agrif_Set_Updatetype(avt_id, update = AGRIF_Update_Full_Weighting) 
     487      CALL Agrif_Set_Updatetype(avm_id, update = AGRIF_Update_Full_Weighting) 
     488   ENDIF 
     489 
     490#else 
    437491   CALL Agrif_Set_Updatetype(tsn_id, update = AGRIF_Update_Average) 
    438  
    439    CALL Agrif_Set_Updatetype(scales_t_id, update = AGRIF_Update_Average) 
    440  
    441492   CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average) 
    442493   CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy) 
    443494 
    444    CALL Agrif_Set_Updatetype(sshn_id, update = AGRIF_Update_Average) 
    445  
    446495   CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Copy, update2 = Agrif_Update_Average) 
    447496   CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Copy) 
     497   CALL Agrif_Set_Updatetype(sshn_id, update = AGRIF_Update_Average) 
     498   CALL Agrif_Set_Updatetype(e3t_id, update = AGRIF_Update_Average) 
    448499 
    449500   IF( ln_zdftke) THEN 
    … …  
    453504   ENDIF 
    454505 
    455 ! High order updates 
    456 !   CALL Agrif_Set_Updatetype(tsn_id, update = Agrif_Update_Full_Weighting) 
    457 !   CALL Agrif_Set_Updatetype(un_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting) 
    458 !   CALL Agrif_Set_Updatetype(vn_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average) 
    459 ! 
    460 !   CALL Agrif_Set_Updatetype(ub2b_update_id,update1 = Agrif_Update_Average, update2 = Agrif_Update_Full_Weighting) 
    461 !   CALL Agrif_Set_Updatetype(vb2b_update_id,update1 = Agrif_Update_Full_Weighting, update2 = Agrif_Update_Average) 
    462 !   CALL Agrif_Set_Updatetype(sshn_id, update = Agrif_Update_Full_Weighting) 
    463  
     506#endif 
    464507   ! 
    465508END SUBROUTINE agrif_declare_var 
    … …  
    500543      CALL ctl_stop('rhot * nn_fsbc(parent) /= N * nn_fsbc(child), therefore nn_fsbc(child) should be set to 1 or nn_fsbc(parent)') 
    501544   ENDIF 
    502  
    503    ! stop if update frequency is different from nn_fsbc 
    504    IF( nbclineupdate > nn_fsbc )  CALL ctl_stop('With ice model on child grid, nn_cln_update should be set to 1 or nn_fsbc') 
    505  
    506  
    507545   ! First Interpolations (using "after" ice subtime step => lim_nbstep=1) 
    508546   !---------------------------------------------------------------------- 
    … …  
    645683      ENDIF 
    646684 
    647       ! Check coordinates 
    648       IF( ln_zps ) THEN 
    649          ! check parameters for partial steps  
    650          IF( Agrif_Parent(e3zps_min) .NE. e3zps_min ) THEN 
    651             WRITE(*,*) 'incompatible e3zps_min between grids' 
    652             WRITE(*,*) 'parent grid :',Agrif_Parent(e3zps_min) 
    653             WRITE(*,*) 'child grid  :',e3zps_min 
    654             WRITE(*,*) 'those values should be identical' 
    655             STOP 
    656          ENDIF 
    657          IF( Agrif_Parent(e3zps_rat) .NE. e3zps_rat ) THEN 
    658             WRITE(*,*) 'incompatible e3zps_rat between grids' 
    659             WRITE(*,*) 'parent grid :',Agrif_Parent(e3zps_rat) 
    660             WRITE(*,*) 'child grid  :',e3zps_rat 
    661             WRITE(*,*) 'those values should be identical'                   
    662             STOP 
    663          ENDIF 
    664685      ENDIF 
    665686      ! Check passive tracer cell 
    … …  
    668689      ENDIF 
    669690   ENDIF 
    670  
    671    CALL Agrif_Update_trc(0) 
    672    ! 
    673    Agrif_UseSpecialValueInUpdate = .FALSE. 
    674    nbcline_trc = 0 
    675691   ! 
    676692END SUBROUTINE Agrif_InitValues_cont_top 
    … …  
    698714   ind2 = 1 + nbghostcells 
    699715   ind3 = 2 + nbghostcells 
    700    CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_id) 
    701    CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_sponge_id) 
     716# if defined key_vertical 
     717   CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_id) 
     718   CALL agrif_declare_variable((/2,2,0,0/),(/ind3,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_sponge_id) 
     719# else 
     720   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_id) 
     721   CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_sponge_id) 
     722# endif 
    702723 
    703724   ! 2. Type of interpolation 
    … …  
    711732   CALL Agrif_Set_bc(trn_sponge_id,(/-nn_sponge_len*Agrif_irhox()-1,0/)) 
    712733 
    713    ! 5. Update type 
     734   ! 4. Update type 
    714735   !---------------  
     736# if defined UPD_HIGH 
     737   CALL Agrif_Set_Updatetype(trn_id, update = Agrif_Update_Full_Weighting) 
     738#else 
    715739   CALL Agrif_Set_Updatetype(trn_id, update = AGRIF_Update_Average) 
    716  
    717 !   Higher order update 
    718 !   CALL Agrif_Set_Updatetype(tsn_id, update = Agrif_Update_Full_Weighting) 
    719  
     740#endif 
    720741   ! 
    721742END SUBROUTINE agrif_declare_var_top 
    … …  
    748769   INTEGER  ::   ios                 ! Local integer output status for namelist read 
    749770   INTEGER  ::   iminspon 
    750    NAMELIST/namagrif/ nn_cln_update, rn_sponge_tra, rn_sponge_dyn, ln_spc_dyn, ln_chk_bathy 
     771   NAMELIST/namagrif/ rn_sponge_tra, rn_sponge_dyn, ln_spc_dyn, ln_chk_bathy 
    751772   !!-------------------------------------------------------------------------------------- 
    752773   ! 
    … …  
    765786      WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    766787      WRITE(numout,*) '   Namelist namagrif : set AGRIF parameters' 
    767       WRITE(numout,*) '      baroclinic update frequency       nn_cln_update = ', nn_cln_update 
    768788      WRITE(numout,*) '      sponge coefficient for tracers    rn_sponge_tra = ', rn_sponge_tra, ' s' 
    769789      WRITE(numout,*) '      sponge coefficient for dynamics   rn_sponge_tra = ', rn_sponge_dyn, ' s' 
    … …  
    774794   ! 
    775795   ! convert DOCTOR namelist name into OLD names 
    776    nbclineupdate = nn_cln_update 
    777796   visc_tra      = rn_sponge_tra 
    778797   visc_dyn      = rn_sponge_dyn 
    … …  
    791810SUBROUTINE Agrif_InvLoc( indloc, nprocloc, i, indglob ) 
    792811   !!---------------------------------------------------------------------- 
    793    !!                     *** ROUTINE Agrif_detect *** 
     812   !!                     *** ROUTINE Agrif_InvLoc *** 
    794813   !!---------------------------------------------------------------------- 
    795814   USE dom_oce 
Note: See TracChangeset for help on using the changeset viewer.