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 14086 for NEMO/trunk/src/NST/agrif_oce_update.F90 – NEMO

Ignore:
Timestamp:
2020-12-04T12:37:14+01:00 (3 years ago)
Author:
cetlod
Message:

Adding AGRIF branches into the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed

  • NEMO/trunk/src/NST/agrif_oce_update.F90

    r14064 r14086  
    1 #undef DECAL_FEEDBACK     /* SEPARATION of INTERFACES */ 
     1#undef DECAL_FEEDBACK    /* SEPARATION of INTERFACES */ 
    22#undef DECAL_FEEDBACK_2D  /* SEPARATION of INTERFACES (Barotropic mode) */ 
    33#undef VOL_REFLUX         /* VOLUME REFLUXING*/ 
    … …  
    2121   USE zdf_oce        ! vertical physics: ocean variables  
    2222   USE agrif_oce 
     23   USE dom_oce 
    2324   ! 
    2425   USE in_out_manager ! I/O manager 
    … …  
    3435 
    3536   PUBLIC   Agrif_Update_Tra, Agrif_Update_Dyn, Agrif_Update_vvl, Agrif_Update_ssh 
    36    PUBLIC   Update_Scales 
     37   PUBLIC   Update_Scales, Agrif_Check_parent_bat 
    3738 
    3839   !! * Substitutions 
    … …  
    5455      IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update tracers  from grid Number',Agrif_Fixed() 
    5556 
    56 #if defined key_vertical 
    57 ! Effect of this has to be carrefully checked  
    58 ! depending on what the nesting tools ensure for 
    59 ! volume conservation: 
     57      l_vremap                      = ln_vert_remap 
     58      Agrif_UseSpecialValueInUpdate = .NOT.l_vremap 
     59      Agrif_SpecialValueFineGrid    = 0._wp 
     60      !  
     61# if ! defined DECAL_FEEDBACK 
     62      CALL Agrif_Update_Variable(ts_update_id, procname=updateTS) 
     63! near boundary update: 
     64!      CALL Agrif_Update_Variable(ts_update_id,locupdate=(/0,2/), procname=updateTS) 
     65# else 
     66      CALL Agrif_Update_Variable(ts_update_id, locupdate=(/1,0/),procname=updateTS) 
     67! near boundary update: 
     68!      CALL Agrif_Update_Variable(ts_update_id,locupdate=(/1,2/), procname=updateTS) 
     69# endif 
     70      ! 
    6071      Agrif_UseSpecialValueInUpdate = .FALSE. 
    61 #else 
    62       Agrif_UseSpecialValueInUpdate = .TRUE. 
    63 #endif 
    64       Agrif_SpecialValueFineGrid    = 0._wp 
    65       !  
    66 # if ! defined DECAL_FEEDBACK 
    67       CALL Agrif_Update_Variable(tsn_id, procname=updateTS) 
    68 ! near boundary update: 
    69 !      CALL Agrif_Update_Variable(tsn_id,locupdate=(/0,2/), procname=updateTS) 
    70 # else 
    71       CALL Agrif_Update_Variable(tsn_id, locupdate=(/1,0/),procname=updateTS) 
    72 ! near boundary update: 
    73 !      CALL Agrif_Update_Variable(tsn_id,locupdate=(/1,2/), procname=updateTS) 
    74 # endif 
    75       ! 
    76       Agrif_UseSpecialValueInUpdate = .FALSE. 
     72      l_vremap                      = .FALSE. 
    7773      ! 
    7874      ! 
    … …  
    9086      Agrif_UseSpecialValueInUpdate = .FALSE. 
    9187      Agrif_SpecialValueFineGrid    = 0._wp 
    92  
    93       use_sign_north = .TRUE. 
    94       sign_north     = -1._wp 
    95  
    96       !      
     88      l_vremap                      = ln_vert_remap 
     89      use_sign_north                = .TRUE. 
     90      sign_north                    = -1._wp      
     91! 
     92# if ! defined DECAL_FEEDBACK_2D 
     93      CALL Agrif_Update_Variable(unb_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateU2d) 
     94      CALL Agrif_Update_Variable(vnb_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateV2d) 
     95# else 
     96      CALL Agrif_Update_Variable(unb_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/1+nn_shift_bar,-2/),procname = updateU2d) 
     97      CALL Agrif_Update_Variable(vnb_update_id,locupdate1=(/1+nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateV2d)   
     98# endif 
     99      !  
     100      IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
     101         ! Update time integrated transports 
     102#  if ! defined DECAL_FEEDBACK_2D 
     103         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updateub2b) 
     104         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updatevb2b) 
     105#  else 
     106         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/  nn_shift_bar,-2/),locupdate2=(/1+nn_shift_bar,-2/),procname = updateub2b) 
     107         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1+nn_shift_bar,-2/),locupdate2=(/  nn_shift_bar,-2/),procname = updatevb2b) 
     108#  endif 
     109      END IF 
     110 
    97111# if ! defined DECAL_FEEDBACK 
    98112      CALL Agrif_Update_Variable(un_update_id,procname = updateU) 
    … …  
    108122!      CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/1,1/),locupdate2=(/0,1/),procname = updateV) 
    109123# endif 
    110  
    111 # if ! defined DECAL_FEEDBACK_2D 
    112       CALL Agrif_Update_Variable(e1u_id,procname = updateU2d) 
    113       CALL Agrif_Update_Variable(e2v_id,procname = updateV2d)   
    114 # else 
    115       CALL Agrif_Update_Variable(e1u_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateU2d) 
    116       CALL Agrif_Update_Variable(e2v_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updateV2d)   
    117 # endif 
    118       ! 
    119 # if ! defined DECAL_FEEDBACK_2D 
    120       ! Account for updated thicknesses at boundary edges 
    121       IF (.NOT.ln_linssh) THEN 
    122 !         CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy) 
    123 !         CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy) 
    124       ENDIF 
    125 # endif 
    126       !  
    127       IF ( ln_dynspg_ts .AND. ln_bt_fw ) THEN 
    128          ! Update time integrated transports 
    129 #  if ! defined DECAL_FEEDBACK_2D 
    130          CALL Agrif_Update_Variable(ub2b_update_id,procname = updateub2b) 
    131          CALL Agrif_Update_Variable(vb2b_update_id,procname = updatevb2b) 
    132 #  else 
    133          CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0,-1/),locupdate2=(/1,-2/),procname = updateub2b) 
    134          CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1,-2/),locupdate2=(/0,-1/),procname = updatevb2b) 
    135 #  endif 
    136       END IF 
    137124      ! 
    138125      use_sign_north = .FALSE. 
     126      l_vremap = .FALSE. 
    139127      ! 
    140128   END SUBROUTINE Agrif_Update_Dyn 
    … …  
    150138      Agrif_SpecialValueFineGrid = 0._wp 
    151139# if ! defined DECAL_FEEDBACK_2D 
    152       CALL Agrif_Update_Variable(sshn_id,procname = updateSSH) 
     140      CALL Agrif_Update_Variable(sshn_id,locupdate=(/  nn_shift_bar,-2/), procname = updateSSH)  
    153141# else 
    154       CALL Agrif_Update_Variable(sshn_id,locupdate=(/1,0/),procname = updateSSH) 
     142      CALL Agrif_Update_Variable(sshn_id,locupdate=(/1+nn_shift_bar,-2/),procname = updateSSH) 
    155143# endif 
    156144      ! 
    … …  
    158146      ! 
    159147#  if defined VOL_REFLUX 
    160       IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN 
     148      IF ( ln_dynspg_ts.AND.ln_bt_fw ) THEN  
    161149         use_sign_north = .TRUE. 
    162150         sign_north = -1._wp 
    163151         ! Refluxing on ssh: 
    164152#  if defined DECAL_FEEDBACK_2D 
    165          CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/0, 0/),locupdate2=(/1, 1/),procname = reflux_sshu) 
    166          CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1, 1/),locupdate2=(/0, 0/),procname = reflux_sshv) 
     153         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/nn_shift_bar,nn_shift_bar/),locupdate2=(/1+nn_shift_bar,1+nn_shift_bar/),procname = reflux_sshu) 
     154         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/1+nn_shift_bar,1+nn_shift_bar/),locupdate2=(/nn_shift_bar,nn_shift_bar/),procname = reflux_sshv) 
    167155#  else 
    168          CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1,-1/),locupdate2=(/ 0, 0/),procname = reflux_sshu) 
    169          CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ 0, 0/),locupdate2=(/-1,-1/),procname = reflux_sshv) 
     156         CALL Agrif_Update_Variable(ub2b_update_id,locupdate1=(/-1+nn_shift_bar,-1+nn_shift_bar/),locupdate2=(/nn_shift_bar, nn_shift_bar/),procname = reflux_sshu) 
     157         CALL Agrif_Update_Variable(vb2b_update_id,locupdate1=(/ nn_shift_bar, nn_shift_bar/),locupdate2=(/-1+nn_shift_bar,-1+nn_shift_bar/),procname = reflux_sshv) 
    170158#  endif 
    171159         use_sign_north = .FALSE. 
    … …  
    320308#endif 
    321309 
    322 #if defined key_vertical 
    323310 
    324311   SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    … …  
    340327      ! 
    341328      IF (before) THEN 
    342 !jc_alt 
    343 !         AGRIF_SpecialValue = -999._wp 
    344          DO jn = n1,n2-1 
     329         IF ( l_vremap ) THEN 
     330            DO jn = n1,n2-1 
     331               DO jk=k1,k2 
     332                  DO jj=j1,j2 
     333                     DO ji=i1,i2 
     334                        tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 
     335                     END DO 
     336                  END DO 
     337               END DO 
     338            END DO 
    345339            DO jk=k1,k2 
    346340               DO jj=j1,j2 
    347341                  DO ji=i1,i2 
    348 !jc_alt 
    349 !                     tabres(ji,jj,jk,jn) = (ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) ) & 
    350 !                                         &  * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1._wp) * 999._wp 
    351                      tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Kmm_a) 
     342                     tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
    352343                  END DO 
    353344               END DO 
    354345            END DO 
    355          END DO 
    356          DO jk=k1,k2 
     346         ELSE 
     347            DO jn = 1,jpts 
     348               DO jk=k1,k2 
     349                  DO jj=j1,j2 
     350                     DO ji=i1,i2 
     351                        tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) 
     352                     END DO 
     353                  END DO 
     354               END DO 
     355            END DO 
     356 
     357         ENDIF 
     358      ELSE 
     359         IF ( l_vremap ) THEN 
     360            tabres_child(:,:,:,:) = 0._wp 
     361            AGRIF_SpecialValue = 0._wp 
    357362            DO jj=j1,j2 
    358363               DO ji=i1,i2 
    359 !jc_alt 
    360 !                  tabres(ji,jj,jk,n2) =      tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) & 
    361 !                                      &   + (tmask(ji,jj,jk) - 1._wp) * 999._wp 
    362                   tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
    363                END DO 
    364             END DO 
    365          END DO 
    366       ELSE 
    367          tabres_child(:,:,:,:) = 0._wp 
    368          AGRIF_SpecialValue = 0._wp 
    369          DO jj=j1,j2 
    370             DO ji=i1,i2 
    371                N_in = 0 
    372                DO jk=k1,k2 !k2 = jpk of child grid 
    373 ! jc_alt 
    374 !                  IF (tabres(ji,jj,jk,n2) < -900._wp  ) EXIT 
    375                   IF (tabres(ji,jj,jk,n2) == 0._wp  ) EXIT 
    376                   N_in = N_in + 1 
    377                   tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
    378                   h_in(N_in) = tabres(ji,jj,jk,n2) 
     364                  N_in = 0 
     365                  DO jk=k1,k2 !k2 = jpk of child grid 
     366                     IF (tabres(ji,jj,jk,n2) <= 1.e-6_wp  ) EXIT 
     367                     N_in = N_in + 1 
     368                     tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 
     369                     h_in(N_in) = tabres(ji,jj,jk,n2) 
     370                  ENDDO 
     371                  N_out = 0 
     372                  DO jk=1,jpk ! jpk of parent grid 
     373                     IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF 
     374                     N_out = N_out + 1 
     375                     h_out(N_out) = e3t(ji,jj,jk,Kmm_a)  
     376                  ENDDO 
     377                  IF (N_in*N_out > 0) THEN !Remove this? 
     378                     CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 
     379                  ENDIF 
    379380               ENDDO 
    380                N_out = 0 
    381                DO jk=1,jpk ! jpk of parent grid 
    382                   IF (tmask(ji,jj,jk) == 0 ) EXIT ! TODO: Will not work with ISF 
    383                   N_out = N_out + 1 
    384                   h_out(N_out) = e3t(ji,jj,jk,Kmm_a)  
    385                ENDDO 
    386                IF (N_in*N_out > 0) THEN !Remove this? 
    387                   CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 
    388                ENDIF 
    389381            ENDDO 
    390          ENDDO 
    391  
    392          IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 
    393             ! Add asselin part 
     382 
     383            IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 
     384               ! Add asselin part 
     385               DO jn = 1,jpts 
     386                  DO jk = 1, jpkm1 
     387                     DO jj = j1, j2 
     388                        DO ji = i1, i2 
     389                           IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 
     390                              ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 
     391                              ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a) 
     392                              ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 
     393                              ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  &  
     394                                        &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 
     395                           ENDIF 
     396                        END DO 
     397                     END DO 
     398                  END DO 
     399               END DO 
     400            ENDIF 
    394401            DO jn = 1,jpts 
    395402               DO jk = 1, jpkm1 
    396403                  DO jj = j1, j2 
    397404                     DO ji = i1, i2 
    398                         IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN 
    399                            ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 
    400                            ztnu = tabres_child(ji,jj,jk,jn) * e3t(ji,jj,jk,Kmm_a) 
    401                            ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 
    402                            ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  &  
    403                                      &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 
    404                         ENDIF 
     405                        IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN  
     406                           ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) 
     407                        END IF 
    405408                     END DO 
    406409                  END DO 
    407410               END DO 
    408411            END DO 
    409          ENDIF 
    410          DO jn = 1,jpts 
    411             DO jk = 1, jpkm1 
    412                DO jj = j1, j2 
    413                   DO ji = i1, i2 
    414                      IF( tabres_child(ji,jj,jk,jn) /= 0._wp ) THEN  
    415                         ts(ji,jj,jk,jn,Kmm_a) = tabres_child(ji,jj,jk,jn) 
    416                      END IF 
     412         ELSE 
     413            DO jn = 1,jpts 
     414               tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 
     415                                            & * tmask(i1:i2,j1:j2,k1:k2) 
     416            ENDDO 
     417  
     418            IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 
     419               ! Add asselin part 
     420               DO jn = 1,jpts 
     421                  DO jk = k1, k2 
     422                     DO jj = j1, j2 
     423                        DO ji = i1, i2 
     424                           IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 
     425                              ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 
     426                              ztnu = tabres(ji,jj,jk,jn) 
     427                              ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 
     428                              ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  &  
     429                                        &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 
     430                           ENDIF 
     431                        END DO 
     432                     END DO 
    417433                  END DO 
    418434               END DO 
    419             END DO 
    420          END DO 
    421          ! 
     435            ENDIF 
     436            DO jn = 1,jpts 
     437               DO jk=k1,k2 
     438                  DO jj=j1,j2 
     439                     DO ji=i1,i2 
     440                        IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN  
     441                           ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a) 
     442                        END IF 
     443                     END DO 
     444                  END DO 
     445               END DO 
     446            END DO 
     447            ! 
     448         ENDIF 
    422449         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    423450            ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kbb_a)  = ts(i1:i2,j1:j2,1:jpkm1,1:jpts,Kmm_a) 
    424          ENDIF 
     451         ENDIF          
    425452      ENDIF 
    426453      !  
    427454   END SUBROUTINE updateTS 
    428455 
    429 # else 
    430  
    431    SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    432       !!--------------------------------------------- 
    433       !!           *** ROUTINE updateT *** 
    434       !!--------------------------------------------- 
    435       INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 
    436       REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    437       LOGICAL, INTENT(in) :: before 
    438       !! 
    439       INTEGER :: ji,jj,jk,jn 
    440       REAL(wp) :: ztb, ztnu, ztno 
    441       !!--------------------------------------------- 
    442       ! 
    443       IF (before) THEN 
    444          DO jn = 1,jpts 
    445             DO jk=k1,k2 
    446                DO jj=j1,j2 
    447                   DO ji=i1,i2 
    448 !> jc tmp 
    449                      tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) / e3t_0(ji,jj,jk) 
    450 !                     tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)  * e3t(ji,jj,jk,Kmm_a) 
    451 !< jc tmp 
    452                   END DO 
    453                END DO 
    454             END DO 
    455          END DO 
    456       ELSE 
    457 !> jc tmp 
    458          DO jn = 1,jpts 
    459             tabres(i1:i2,j1:j2,k1:k2,jn) =  tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 
    460                                          & * tmask(i1:i2,j1:j2,k1:k2) 
    461          ENDDO 
    462 !< jc tmp 
    463          IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN 
    464             ! Add asselin part 
    465             DO jn = 1,jpts 
    466                DO jk = k1, k2 
    467                   DO jj = j1, j2 
    468                      DO ji = i1, i2 
    469                         IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN 
    470                            ztb  = ts(ji,jj,jk,jn,Kbb_a) * e3t(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 
    471                            ztnu = tabres(ji,jj,jk,jn) 
    472                            ztno = ts(ji,jj,jk,jn,Kmm_a) * e3t(ji,jj,jk,Krhs_a) 
    473                            ts(ji,jj,jk,jn,Kbb_a) = ( ztb + rn_atfp * ( ztnu - ztno) )  &  
    474                                      &        * tmask(ji,jj,jk) / e3t(ji,jj,jk,Kbb_a) 
    475                         ENDIF 
    476                      END DO 
    477                   END DO 
    478                END DO 
    479             END DO 
    480          ENDIF 
    481          DO jn = 1,jpts 
    482             DO jk=k1,k2 
    483                DO jj=j1,j2 
    484                   DO ji=i1,i2 
    485                      IF( tabres(ji,jj,jk,jn) /= 0._wp ) THEN  
    486                         ts(ji,jj,jk,jn,Kmm_a) = tabres(ji,jj,jk,jn) / e3t(ji,jj,jk,Kmm_a) 
    487                      END IF 
    488                   END DO 
    489                END DO 
    490             END DO 
    491          END DO 
    492          ! 
    493          IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    494             ts(i1:i2,j1:j2,k1:k2,1:jpts,Kbb_a)  = ts(i1:i2,j1:j2,k1:k2,1:jpts,Kmm_a) 
    495          ENDIF 
    496          ! 
    497       ENDIF 
    498       !  
    499    END SUBROUTINE updateTS 
    500  
    501 # endif 
    502  
    503 # if defined key_vertical 
    504456 
    505457   SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    … …  
    513465      INTEGER ::   ji, jj, jk 
    514466      REAL(wp)::   zrhoy, zub, zunu, zuno 
     467      REAL(wp), DIMENSION(jpi,jpj) ::   zpgu   ! 2D workspace 
    515468! VERTICAL REFINEMENT BEGIN 
    516469      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
    517470      REAL(wp) :: h_in(k1:k2) 
    518471      REAL(wp) :: h_out(1:jpk) 
    519       INTEGER  :: N_in, N_out 
    520       REAL(wp) :: h_diff, excess, thick 
     472      INTEGER  :: N_in, N_out, N_in_save, N_out_save 
     473      REAL(wp) :: zhmin, zd 
    521474      REAL(wp) :: tabin(k1:k2) 
    522475! VERTICAL REFINEMENT END 
    … …  
    525478      IF( before ) THEN 
    526479         zrhoy = Agrif_Rhoy() 
    527 !jc_alt 
    528 !         AGRIF_SpecialValue = -999._wp 
    529480         DO jk=k1,k2 
     481            tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) &  
     482                                     &   * umask(i1:i2,j1:j2,jk) * uu(i1:i2,j1:j2,jk,Kmm_a)   
     483         END DO 
     484 
     485         IF ( l_vremap ) THEN 
     486            DO jk=k1,k2 
     487               tabres(i1:i2,j1:j2,jk,2) = zrhoy * umask(i1:i2,j1:j2,jk) * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) 
     488            END DO 
     489         ENDIF 
     490 
     491      ELSE 
     492 
     493         tabres_child(:,:,:) = 0._wp 
     494 
     495         IF ( l_vremap ) THEN 
     496 
    530497            DO jj=j1,j2 
    531498               DO ji=i1,i2 
    532 !jc_alt 
    533 !                  tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a)  & 
    534 !                                     &  + (umask(ji,jj,jk)-1._wp)*999._wp 
    535                   tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) * uu(ji,jj,jk,Kmm_a)   
    536 !jc_alt 
    537 !                  tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)  & 
    538 !                                     &  + (umask(ji,jj,jk)-1._wp)*999._wp 
    539                   tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 
    540                END DO 
    541             END DO 
    542          END DO 
    543       ELSE 
    544          tabres_child(:,:,:) = 0. 
    545          AGRIF_SpecialValue = 0._wp 
    546          DO jj=j1,j2 
    547             DO ji=i1,i2 
    548                N_in = 0 
    549                h_in(:) = 0._wp 
    550                tabin(:) = 0._wp 
    551                DO jk=k1,k2 !k2=jpk of child grid 
    552 !jc_alt 
    553 !                  IF( tabres(ji,jj,jk,2) < -900._wp) EXIT 
    554                   IF( tabres(ji,jj,jk,2) == 0.) EXIT 
    555                   N_in = N_in + 1 
    556                   tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
    557                   h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 
     499                  N_in = 0 
     500                  h_in(:) = 0._wp 
     501                  tabin(:) = 0._wp 
     502                  DO jk=k1,k2 !k2=jpk of child grid 
     503                     IF( tabres(ji,jj,jk,2)*r1_e2u(ji,jj) <= 1.e-6_wp ) EXIT 
     504                     N_in = N_in + 1 
     505                     tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     506                     h_in(N_in) = tabres(ji,jj,jk,2) * r1_e2u(ji,jj) 
     507                  ENDDO 
     508                  N_out = 0 
     509                  DO jk=1,jpk 
     510                     IF (umask(ji,jj,jk) == 0._wp) EXIT 
     511                     N_out = N_out + 1 
     512                     h_out(N_out) = e3u(ji,jj,jk,Kmm_a) 
     513                  ENDDO 
     514                  IF (N_in * N_out > 0) THEN 
     515                     ! Deal with potentially different depths at velocity points: 
     516                     N_in_save  = N_in 
     517                     N_out_save = N_out 
     518                     IF ( ABS(sum(h_out(1:N_out))-sum(h_in(1:N_in))) > 1.e-6_wp ) THEN 
     519                        zhmin = MIN(sum(h_out(1:N_out)), sum(h_in(1:N_in))) 
     520                        zd = 0._wp 
     521                        DO jk=1, N_in_save 
     522                           IF ( (zd +  h_in(jk)) > zhmin-1.e-6) THEN 
     523                              N_in = jk 
     524                              h_in(jk) = zhmin - zd 
     525                              EXIT  
     526                           ENDIF 
     527                           zd = zd + h_in(jk) 
     528                        END DO 
     529                        zd = 0._wp 
     530                        DO jk=1, N_out_save 
     531                           IF ( (zd +  h_out(jk)) > zhmin-1.e-6) THEN 
     532                              N_out = jk 
     533                              h_out(jk) = zhmin - zd 
     534                              EXIT  
     535                           ENDIF 
     536                           zd = zd + h_out(jk) 
     537                        END DO 
     538                     END IF 
     539                     CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
     540                     IF (N_out < N_out_save) tabres_child(ji,jj,N_out+1:N_out_save) = tabres_child(ji,jj,N_out) 
     541                  ENDIF 
    558542               ENDDO 
    559                N_out = 0 
    560                DO jk=1,jpk 
    561                   IF (umask(ji,jj,jk) == 0) EXIT 
    562                   N_out = N_out + 1 
    563                   h_out(N_out) = e3u(ji,jj,jk,Kmm_a) 
    564                ENDDO 
    565                IF (N_in * N_out > 0) THEN 
    566                   h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    567                   excess = 0._wp 
    568                   IF (h_diff < -1.e-4) THEN 
    569 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid.  
    570 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 
    571                      DO jk=N_in,1,-1 
    572                         thick = MIN(-1*h_diff, h_in(jk)) 
    573                         excess = excess + tabin(jk)*thick*e2u(ji,jj) 
    574                         tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 
    575                         h_diff = h_diff + thick 
    576                         IF ( h_diff == 0) THEN 
    577                            N_in = jk 
    578                            h_in(jk) = h_in(jk) - thick 
    579                            EXIT 
    580                         ENDIF 
    581                      ENDDO 
    582                   ENDIF 
    583                   CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
    584                   tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 
    585                ENDIF 
    586543            ENDDO 
    587          ENDDO 
     544 
     545         ELSE 
     546            DO jk=1,jpk 
     547               DO jj=j1,j2 
     548                  DO ji=i1,i2 
     549                     tabres_child(ji,jj,jk) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) /  e3u(ji,jj,jk,Kmm_a) 
     550                  END DO 
     551               END DO 
     552            END DO 
     553         ENDIF 
    588554         ! 
    589555         DO jk=1,jpk 
    … …  
    603569         END DO 
    604570         ! 
     571         ! Correct now and before transports: 
     572         DO jj=j1,j2 
     573            DO ji=i1,i2 
     574               zpgu(ji,jj) = 0._wp 
     575               DO jk=1,jpkm1 
     576                  zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a) 
     577               END DO 
     578               ! 
     579               DO jk=1,jpkm1               
     580                  uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + & 
     581                      &  (uu_b(ji,jj,Kmm_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kmm_a)) * umask(ji,jj,jk)            
     582               END DO 
     583               ! 
     584               zpgu(ji,jj) = 0._wp 
     585               DO jk=1,jpkm1 
     586                  zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) 
     587               END DO 
     588               ! 
     589               DO jk=1,jpkm1               
     590                  uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + & 
     591                      &  (uu_b(ji,jj,Kbb_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kbb_a)) * umask(ji,jj,jk)            
     592               END DO 
     593               ! 
     594            END DO 
     595         END DO 
     596         ! 
    605597         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    606598            uu(i1:i2,j1:j2,1:jpkm1,Kbb_a)  = uu(i1:i2,j1:j2,1:jpkm1,Kmm_a) 
    … …  
    611603   END SUBROUTINE updateu 
    612604 
    613 #else 
    614  
    615    SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    616       !!--------------------------------------------- 
    617       !!           *** ROUTINE updateu *** 
    618       !!--------------------------------------------- 
    619       INTEGER                               , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    620       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    621       LOGICAL                                     , INTENT(in   ) :: before 
    622       ! 
    623       INTEGER  :: ji, jj, jk 
    624       REAL(wp) :: zrhoy, zub, zunu, zuno 
    625       !!--------------------------------------------- 
    626       !  
    627       IF( before ) THEN 
    628          zrhoy = Agrif_Rhoy() 
    629          DO jk = k1, k2 
    630             tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u(i1:i2,j1:j2,jk,Kmm_a) * uu(i1:i2,j1:j2,jk,Kmm_a) 
    631          END DO 
    632       ELSE 
    633          DO jk=k1,k2 
    634             DO jj=j1,j2 
    635                DO ji=i1,i2 
    636                   tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj)  
    637                   ! 
    638                   IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part 
    639                      zub  = uu(ji,jj,jk,Kbb_a) * e3u(ji,jj,jk,Kbb_a)  ! fse3t_b prior update should be used 
    640                      zuno = uu(ji,jj,jk,Kmm_a) * e3u(ji,jj,jk,Krhs_a) 
    641                      zunu = tabres(ji,jj,jk,1) 
    642                      uu(ji,jj,jk,Kbb_a) = ( zub + rn_atfp * ( zunu - zuno) ) &       
    643                                     & * umask(ji,jj,jk) / e3u(ji,jj,jk,Kbb_a) 
    644                   ENDIF 
    645                   ! 
    646                   uu(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u(ji,jj,jk,Kmm_a) 
    647                END DO 
    648             END DO 
    649          END DO 
    650          ! 
    651          IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    652             uu(i1:i2,j1:j2,k1:k2,Kbb_a)  = uu(i1:i2,j1:j2,k1:k2,Kmm_a) 
    653          ENDIF 
    654          ! 
    655       ENDIF 
    656       !  
    657    END SUBROUTINE updateu 
    658  
    659 # endif 
    660  
    661    SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
    662       !!--------------------------------------------- 
    663       !!           *** ROUTINE correct_u_bdy *** 
     605 
     606   SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
     607      !!--------------------------------------------- 
     608      !!           *** ROUTINE updatev *** 
    664609      !!--------------------------------------------- 
    665610      INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    666611      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    667612      LOGICAL                                     , INTENT(in   ) :: before 
    668       INTEGER                                     , INTENT(in)    :: nb, ndir 
    669       !! 
    670       LOGICAL :: western_side, eastern_side  
    671       ! 
    672       INTEGER  :: jj, jk 
    673       REAL(wp) :: zcor 
    674       !!--------------------------------------------- 
    675       !  
    676       IF( .NOT.before ) THEN 
    677          ! 
    678          western_side  = (nb == 1).AND.(ndir == 1) 
    679          eastern_side  = (nb == 1).AND.(ndir == 2) 
    680          ! 
    681          IF (western_side) THEN 
    682             DO jj=j1,j2 
    683                zcor = uu_b(i1-1,jj,Kmm_a) * hu(i1-1,jj,Krhs_a) * r1_hu(i1-1,jj,Kmm_a) - uu_b(i1-1,jj,Kmm_a) 
    684                uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor 
    685                DO jk=1,jpkm1 
    686                   uu(i1-1,jj,jk,Kmm_a) = uu(i1-1,jj,jk,Kmm_a) + zcor * umask(i1-1,jj,jk) 
    687                END DO  
    688             END DO 
    689          ENDIF 
    690          ! 
    691          IF (eastern_side) THEN 
    692             DO jj=j1,j2 
    693                zcor = uu_b(i2+1,jj,Kmm_a) * hu(i2+1,jj,Krhs_a) * r1_hu(i2+1,jj,Kmm_a) - uu_b(i2+1,jj,Kmm_a) 
    694                uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor 
    695                DO jk=1,jpkm1 
    696                   uu(i2+1,jj,jk,Kmm_a) = uu(i2+1,jj,jk,Kmm_a) + zcor * umask(i2+1,jj,jk) 
    697                END DO  
    698             END DO 
    699          ENDIF 
    700          ! 
    701       ENDIF 
    702       !  
    703    END SUBROUTINE correct_u_bdy 
    704  
    705 # if  defined key_vertical 
    706  
    707    SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 
    708       !!--------------------------------------------- 
    709       !!           *** ROUTINE updatev *** 
    710       !!--------------------------------------------- 
    711       INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    712       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    713       LOGICAL                                     , INTENT(in   ) :: before 
    714613      ! 
    715614      INTEGER  ::   ji, jj, jk 
    716615      REAL(wp) ::   zrhox, zvb, zvnu, zvno 
     616      REAL(wp), DIMENSION(jpi,jpj) ::   zpgv   ! 2D workspace 
    717617! VERTICAL REFINEMENT BEGIN 
    718618      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 
    719619      REAL(wp) :: h_in(k1:k2) 
    720620      REAL(wp) :: h_out(1:jpk) 
    721       INTEGER :: N_in, N_out 
    722       REAL(wp) :: h_diff, excess, thick 
     621      INTEGER  :: N_in, N_out, N_in_save, N_out_save 
     622      REAL(wp) :: zhmin, zd 
    723623      REAL(wp) :: tabin(k1:k2) 
    724624! VERTICAL REFINEMENT END 
    … …  
    727627      IF( before ) THEN 
    728628         zrhox = Agrif_Rhox() 
    729 !jc_alt 
    730 !         AGRIF_SpecialValue = -999._wp 
    731629         DO jk=k1,k2 
     630            tabres(i1:i2,j1:j2,jk,1) = zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Kmm_a) &  
     631                                     &   * vmask(i1:i2,j1:j2,jk) * vv(i1:i2,j1:j2,jk,Kmm_a)   
     632         END DO 
     633 
     634         IF ( l_vremap ) THEN 
     635            DO jk=k1,k2 
     636               tabres(i1:i2,j1:j2,jk,2) = zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Kmm_a) * vmask(i1:i2,j1:j2,jk) 
     637            END DO 
     638         ENDIF 
     639 
     640      ELSE 
     641 
     642         tabres_child(:,:,:) = 0._wp 
     643 
     644         IF ( l_vremap ) THEN 
     645 
    732646            DO jj=j1,j2 
    733647               DO ji=i1,i2 
    734 !jc_alt 
    735 !                  tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a) & 
    736 !                                     & + (vmask(ji,jj,jk)-1._wp) * 999._wp 
    737                   tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) * vv(ji,jj,jk,Kmm_a)  
    738 !jc_alt 
    739 !                  tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) & 
    740 !                                     & + (vmask(ji,jj,jk)-1._wp) * 999._wp 
    741                   tabres(ji,jj,jk,2) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) 
    742                END DO 
    743             END DO 
    744          END DO 
    745       ELSE 
    746          tabres_child(:,:,:) = 0. 
    747          AGRIF_SpecialValue = 0._wp 
    748          DO jj=j1,j2 
    749             DO ji=i1,i2 
    750                N_in = 0 
    751                DO jk=k1,k2 
    752 !jc_alt 
    753 !                  IF (tabres(ji,jj,jk,2) < -900._wp) EXIT 
    754                   IF (tabres(ji,jj,jk,2) == 0) EXIT 
    755                   N_in = N_in + 1 
    756                   tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
    757                   h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) 
     648                  N_in = 0 
     649                  DO jk=k1,k2 
     650                     IF (tabres(ji,jj,jk,2)* r1_e1v(ji,jj) <= 1.e-6_wp) EXIT 
     651                     N_in = N_in + 1 
     652                     tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 
     653                     h_in(N_in) = tabres(ji,jj,jk,2) * r1_e1v(ji,jj) 
     654                  ENDDO 
     655                  N_out = 0 
     656                  DO jk=1,jpk 
     657                     IF (vmask(ji,jj,jk) == 0) EXIT 
     658                     N_out = N_out + 1 
     659                     h_out(N_out) = e3v(ji,jj,jk,Kmm_a) 
     660                  ENDDO 
     661                  IF (N_in * N_out > 0) THEN 
     662                     ! Deal with potentially different depths at velocity points: 
     663                     N_in_save  = N_in 
     664                     N_out_save = N_out 
     665                     IF ( ABS(sum(h_out(1:N_out))-sum(h_in(1:N_in))) > 1.e-6_wp ) THEN 
     666                        zhmin = MIN(sum(h_out(1:N_out)), sum(h_in(1:N_in))) 
     667                        zd = 0._wp 
     668                        DO jk=1, N_in_save 
     669                           IF ( (zd +  h_in(jk)) > zhmin-1.e-6) THEN 
     670                              N_in = jk 
     671                              h_in(jk) = zhmin - zd 
     672                              EXIT  
     673                           ENDIF 
     674                           zd = zd + h_in(jk) 
     675                        END DO 
     676                        zd = 0._wp 
     677                        DO jk=1, N_out_save 
     678                           IF ( (zd +  h_out(jk)) > zhmin-1.e-6) THEN 
     679                              N_out = jk 
     680                              h_out(jk) = zhmin - zd 
     681                              EXIT  
     682                           ENDIF 
     683                           zd = zd + h_out(jk) 
     684                        END DO 
     685                     END IF 
     686                     CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
     687                     IF (N_out < N_out_save) tabres_child(ji,jj,N_out+1:N_out_save) = tabres_child(ji,jj,N_out) 
     688                  ENDIF 
    758689               ENDDO 
    759                N_out = 0 
    760                DO jk=1,jpk 
    761                   IF (vmask(ji,jj,jk) == 0) EXIT 
    762                   N_out = N_out + 1 
    763                   h_out(N_out) = e3v(ji,jj,jk,Kmm_a) 
    764                ENDDO 
    765                IF (N_in * N_out > 0) THEN 
    766                   h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 
    767                   excess = 0._wp 
    768                   IF (h_diff < -1.e-4) then 
    769 !Even if bathy at T points match it's possible for the V points to be deeper in the child grid.  
    770 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 
    771                      DO jk=N_in,1,-1 
    772                         thick = MIN(-1*h_diff, h_in(jk)) 
    773                         excess = excess + tabin(jk)*thick*e2u(ji,jj) 
    774                         tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 
    775                         h_diff = h_diff + thick 
    776                         IF ( h_diff == 0) THEN 
    777                            N_in = jk 
    778                            h_in(jk) = h_in(jk) - thick 
    779                            EXIT 
    780                         ENDIF 
    781                      ENDDO 
    782                   ENDIF 
    783                   CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 
    784                   tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 
    785                ENDIF 
    786690            ENDDO 
    787          ENDDO 
     691 
     692         ELSE 
     693            DO jk=1,jpk 
     694               DO jj=j1,j2 
     695                  DO ji=i1,i2 
     696                     tabres_child(ji,jj,jk) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) /  e3v(ji,jj,jk,Kmm_a) 
     697                  END DO 
     698               END DO 
     699            END DO 
     700         ENDIF 
    788701         ! 
    789702         DO jk=1,jpkm1 
    … …  
    803716         END DO 
    804717         ! 
     718         ! Correct now and before transports: 
     719         DO jj=j1,j2 
     720            DO ji=i1,i2 
     721               zpgv(ji,jj) = 0._wp 
     722               DO jk=1,jpkm1 
     723                  zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 
     724               END DO 
     725               ! 
     726               DO jk=1,jpkm1               
     727                  vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + & 
     728                      &  (vv_b(ji,jj,Kmm_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kmm_a)) * vmask(ji,jj,jk)            
     729               END DO 
     730               ! 
     731               zpgv(ji,jj) = 0._wp 
     732               DO jk=1,jpkm1 
     733                  zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) 
     734               END DO 
     735               ! 
     736               DO jk=1,jpkm1               
     737                  vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + & 
     738                      &  (vv_b(ji,jj,Kbb_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kbb_a)) * vmask(ji,jj,jk)            
     739               END DO 
     740               ! 
     741            END DO 
     742         END DO 
     743         ! 
    805744         IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    806745            vv(i1:i2,j1:j2,1:jpkm1,Kbb_a)  = vv(i1:i2,j1:j2,1:jpkm1,Kmm_a) 
    … …  
    810749      !  
    811750   END SUBROUTINE updatev 
    812  
    813 # else 
    814  
    815    SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) 
    816       !!--------------------------------------------- 
    817       !!           *** ROUTINE updatev *** 
    818       !!--------------------------------------------- 
    819       INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    820       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    821       LOGICAL                                     , INTENT(in   ) :: before 
    822       ! 
    823       INTEGER  :: ji, jj, jk 
    824       REAL(wp) :: zrhox, zvb, zvnu, zvno 
    825       !!---------------------------------------------       
    826       ! 
    827       IF (before) THEN 
    828          zrhox = Agrif_Rhox() 
    829          DO jk=k1,k2 
    830             DO jj=j1,j2 
    831                DO ji=i1,i2 
    832                   tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 
    833                END DO 
    834             END DO 
    835          END DO 
    836       ELSE 
    837          DO jk=k1,k2 
    838             DO jj=j1,j2 
    839                DO ji=i1,i2 
    840                   tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) 
    841                   ! 
    842                   IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) THEN ! Add asselin part 
    843                      zvb  = vv(ji,jj,jk,Kbb_a) * e3v(ji,jj,jk,Kbb_a) ! fse3t_b prior update should be used 
    844                      zvno = vv(ji,jj,jk,Kmm_a) * e3v(ji,jj,jk,Krhs_a) 
    845                      zvnu = tabres(ji,jj,jk,1) 
    846                      vv(ji,jj,jk,Kbb_a) = ( zvb + rn_atfp * ( zvnu - zvno) ) &       
    847                                     & * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kbb_a) 
    848                   ENDIF 
    849                   ! 
    850                   vv(ji,jj,jk,Kmm_a) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v(ji,jj,jk,Kmm_a) 
    851                END DO 
    852             END DO 
    853          END DO 
    854          ! 
    855          IF  ((l_1st_euler).AND.(Agrif_Nb_Step()==0) ) THEN 
    856             vv(i1:i2,j1:j2,k1:k2,Kbb_a)  = vv(i1:i2,j1:j2,k1:k2,Kmm_a) 
    857          ENDIF 
    858          ! 
    859       ENDIF 
    860       !  
    861    END SUBROUTINE updatev 
    862  
    863 # endif 
    864  
    865    SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 
    866       !!--------------------------------------------- 
    867       !!           *** ROUTINE correct_v_bdy *** 
    868       !!--------------------------------------------- 
    869       INTEGER                                     , INTENT(in   ) :: i1, i2, j1, j2, k1, k2, n1, n2 
    870       REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 
    871       LOGICAL                                     , INTENT(in   ) :: before 
    872       INTEGER                                     , INTENT(in)    :: nb, ndir 
    873       !! 
    874       LOGICAL :: southern_side, northern_side  
    875       ! 
    876       INTEGER  :: ji, jk 
    877       REAL(wp) :: zcor 
    878       !!--------------------------------------------- 
    879       !  
    880       IF( .NOT.before ) THEN 
    881          ! 
    882          southern_side = (nb == 2).AND.(ndir == 1) 
    883          northern_side = (nb == 2).AND.(ndir == 2) 
    884          ! 
    885          IF (southern_side) THEN 
    886             DO ji=i1,i2 
    887                zcor = vv_b(ji,j1-1,Kmm_a) * hv(ji,j1-1,Krhs_a) * r1_hv(ji,j1-1,Kmm_a) - vv_b(ji,j1-1,Kmm_a) 
    888                vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor 
    889                DO jk=1,jpkm1 
    890                   vv(ji,j1-1,jk,Kmm_a) = vv(ji,j1-1,jk,Kmm_a) + zcor * vmask(ji,j1-1,jk) 
    891                END DO  
    892             END DO 
    893          ENDIF 
    894          ! 
    895          IF (northern_side) THEN 
    896             DO ji=i1,i2 
    897                zcor = vv_b(ji,j2+1,Kmm_a) * hv(ji,j2+1,Krhs_a) * r1_hv(ji,j2+1,Kmm_a) - vv_b(ji,j2+1,Kmm_a) 
    898                vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor 
    899                DO jk=1,jpkm1 
    900                   vv(ji,j2+1,jk,Kmm_a) = vv(ji,j2+1,jk,Kmm_a) + zcor * vmask(ji,j2+1,jk) 
    901                END DO  
    902             END DO 
    903          ENDIF 
    904          ! 
    905       ENDIF 
    906       !  
    907    END SUBROUTINE correct_v_bdy 
    908751 
    909752 
    … …  
    935778               tabres(ji,jj) =  tabres(ji,jj) * r1_e2u(ji,jj)   
    936779               !     
    937                ! Update "now" 3d velocities: 
    938                zpgu(ji,jj) = 0._wp 
    939                DO jk=1,jpkm1 
    940                   zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kmm_a) * uu(ji,jj,jk,Kmm_a) 
    941                END DO 
    942                ! 
    943                zcorr = (tabres(ji,jj) - zpgu(ji,jj)) * r1_hu(ji,jj,Kmm_a) 
    944                DO jk=1,jpkm1               
    945                   uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk)            
    946                END DO 
    947                ! 
    948780               ! Update barotropic velocities: 
    949781               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    … …  
    955787               uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1) 
    956788               !        
    957                ! Correct "before" velocities to hold correct bt component: 
    958                zpgu(ji,jj) = 0.e0 
    959                DO jk=1,jpkm1 
    960                   zpgu(ji,jj) = zpgu(ji,jj) + e3u(ji,jj,jk,Kbb_a) * uu(ji,jj,jk,Kbb_a) 
    961                END DO 
    962                ! 
    963                zcorr = uu_b(ji,jj,Kbb_a) - zpgu(ji,jj) * r1_hu(ji,jj,Kbb_a) 
    964                DO jk=1,jpkm1               
    965                   uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk)            
    966                END DO 
    967                ! 
    968789            END DO 
    969790         END DO 
    … …  
    1002823            DO ji=i1,i2 
    1003824               tabres(ji,jj) =  tabres(ji,jj) * r1_e1v(ji,jj)   
    1004                !     
    1005                ! Update "now" 3d velocities: 
    1006                zpgv(ji,jj) = 0.e0 
    1007                DO jk=1,jpkm1 
    1008                   zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kmm_a) * vv(ji,jj,jk,Kmm_a) 
    1009                END DO 
    1010                ! 
    1011                zcorr = (tabres(ji,jj) - zpgv(ji,jj)) * r1_hv(ji,jj,Kmm_a) 
    1012                DO jk=1,jpkm1               
    1013                   vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk)            
    1014                END DO 
    1015                ! 
    1016825               ! Update barotropic velocities: 
    1017826               IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 
    … …  
    1023832               vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1) 
    1024833               !        
    1025                ! Correct "before" velocities to hold correct bt component: 
    1026                zpgv(ji,jj) = 0.e0 
    1027                DO jk=1,jpkm1 
    1028                   zpgv(ji,jj) = zpgv(ji,jj) + e3v(ji,jj,jk,Kbb_a) * vv(ji,jj,jk,Kbb_a) 
    1029                END DO 
    1030                ! 
    1031                zcorr = vv_b(ji,jj,Kbb_a) - zpgv(ji,jj) * r1_hv(ji,jj,Kbb_a) 
    1032                DO jk=1,jpkm1               
    1033                   vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk)            
    1034                END DO 
    1035                ! 
    1036834            END DO 
    1037835         END DO 
    … …  
    1120918               ub2_i_b(ji,jj) = ub2_i_b(ji,jj) + za1 * zcor  
    1121919               ! Update corrective fluxes: 
    1122                un_bf(ji,jj)  = un_bf(ji,jj) + zcor 
     920               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler))) un_bf(ji,jj)  = un_bf(ji,jj) + zcor 
    1123921               ! Update half step back fluxes: 
    1124922               ub2_b(ji,jj) = tabres(ji,jj) 
    … …  
    12081006               vb2_i_b(ji,jj) = vb2_i_b(ji,jj) + za1 * zcor  
    12091007               ! Update corrective fluxes: 
    1210                vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor 
     1008               IF (.NOT.(lk_agrif_fstep.AND.(l_1st_euler)))  vn_bf(ji,jj)  = vn_bf(ji,jj) + zcor 
    12111009               ! Update half step back fluxes: 
    12121010               vb2_b(ji,jj) = tabres(ji,jj) 
    … …  
    14791277#endif 
    14801278 
     1279   SUBROUTINE Agrif_Check_parent_bat( ) 
     1280      !!---------------------------------------------------------------------- 
     1281      !!                   *** ROUTINE Agrif_Check_parent_bat *** 
     1282      !!---------------------------------------------------------------------- 
     1283      !  
     1284      IF (( .NOT.ln_agrif_2way ).OR.(.NOT.ln_chk_bathy).OR.(Agrif_Root())) RETURN 
     1285      ! 
     1286      Agrif_UseSpecialValueInUpdate = .FALSE. 
     1287      ! 
     1288      IF(lwp) WRITE(numout,*) ' ' 
     1289      IF(lwp) WRITE(numout,*) 'AGRIF: Check parent volume at Level:', Agrif_Level() 
     1290      ! 
     1291# if ! defined DECAL_FEEDBACK 
     1292      CALL Agrif_Update_Variable(batupd_id,procname = update_bat) 
     1293# else 
     1294      CALL Agrif_Update_Variable(batupd_id,locupdate=(/1,0/),procname = update_bat) 
     1295# endif 
     1296      ! 
     1297      kindic_agr = Agrif_Parent(kindic_agr) 
     1298      CALL mpp_sum( 'Agrif_Check_parent_bat', kindic_agr ) 
     1299 
     1300      IF( kindic_agr /= 0 ) THEN 
     1301         CALL ctl_stop('==> Averaged Bathymetry does not match parent volume')  
     1302      ELSE 
     1303         IF(lwp) WRITE(numout,*) '==> Averaged Bathymetry matches parent '  
     1304         IF(lwp) WRITE(numout,*) '' 
     1305      ENDIF 
     1306      ! 
     1307   END SUBROUTINE Agrif_Check_parent_bat 
     1308 
     1309   SUBROUTINE update_bat(ptab, i1, i2, j1, j2, before ) 
     1310      !!--------------------------------------------- 
     1311      !!           *** ROUTINE update_bat *** 
     1312      !!--------------------------------------------- 
     1313      REAL(wp), DIMENSION(i1:i2,j1:j2) :: ptab 
     1314      INTEGER, INTENT(in) :: i1, i2, j1, j2 
     1315      LOGICAL, INTENT(in) :: before 
     1316      INTEGER :: ji, jj 
     1317      ! 
     1318      !!--------------------------------------------- 
     1319      ! 
     1320      IF( before ) THEN 
     1321         ptab(i1:i2,j1:j2) = ht_0(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1) 
     1322      ELSE 
     1323         kindic_agr = 0 
     1324         ! 
     1325         DO jj=j1,j2 
     1326            DO ji=i1,i2 
     1327               IF ( (ssmask(ji,jj).NE.0._wp).AND.& 
     1328               & (ABS(ptab(ji,jj)-ht_0(ji,jj)).GE.1.e-6) ) THEN  
     1329                  kindic_agr = kindic_agr + 1  
     1330               ENDIF 
     1331            END DO 
     1332         END DO 
     1333         ! 
     1334      ENDIF 
     1335      ! 
     1336   END SUBROUTINE update_bat 
     1337 
    14811338#else 
    14821339   !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.