Skip to content
Snippets Groups Projects
dynspg_ts.F90 77.8 KiB
Newer Older
Guillaume Samson's avatar
Guillaume Samson committed
            CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b'  , vb2_i_b(:,:) )
         ENDIF
#endif
      ENDIF
      !
   END SUBROUTINE ts_rst


   SUBROUTINE dyn_spg_ts_init
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE dyn_spg_ts_init  ***
      !!
      !! ** Purpose : Set time splitting options
      !!----------------------------------------------------------------------
      INTEGER  ::   ji ,jj              ! dummy loop indices
      REAL(wp) ::   zxr2, zyr2, zcmax   ! local scalar
      REAL(wp), DIMENSION(jpi,jpj) ::   zcu
      !!----------------------------------------------------------------------
      !
      ! Max courant number for ext. grav. waves
      !
      DO_2D( 0, 0, 0, 0 )
         zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj)
         zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj)
         zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) )
      END_2D
      !
#if defined key_agrif
     ! Discard points that are not stepped by 2d mode:
     zcu(Nis0:Nie0,Njs0:Nje0) = zcu(Nis0:Nie0,Njs0:Nje0) &
                              & * (1._wp - tmask_upd(Nis0:Nie0,Njs0:Nje0))
#endif
      !
Guillaume Samson's avatar
Guillaume Samson committed
      zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) )
      CALL mpp_max( 'dynspg_ts', zcmax )

      ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax
      IF( ln_bt_auto )   nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax)
      
Guillaume Samson's avatar
Guillaume Samson committed
      zcmax = zcmax * rDt_e
      ! Print results
      IF(lwp) WRITE(numout,*)
      IF(lwp) WRITE(numout,*) 'dyn_spg_ts_init : split-explicit free surface'
      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
      IF( ln_bt_auto ) THEN
         IF(lwp) WRITE(numout,*) '     ln_ts_auto =.true. Automatically set nn_e '
         IF(lwp) WRITE(numout,*) '     Max. courant number allowed: ', rn_bt_cmax
      ELSE
         IF(lwp) WRITE(numout,*) '     ln_ts_auto=.false.: Use nn_e in namelist   nn_e = ', nn_e
      ENDIF

      IF(ln_bt_av) THEN
         IF(lwp) WRITE(numout,*) '     ln_bt_av =.true.  ==> Time averaging over nn_e time steps is on '
      ELSE
         IF(lwp) WRITE(numout,*) '     ln_bt_av =.false. => No time averaging of barotropic variables '
      ENDIF
      !
      !
      IF(ln_bt_fw) THEN
         IF(lwp) WRITE(numout,*) '     ln_bt_fw=.true.  => Forward integration of barotropic variables '
      ELSE
         IF(lwp) WRITE(numout,*) '     ln_bt_fw =.false.=> Centred integration of barotropic variables '
      ENDIF
      !
#if defined key_agrif
      ! Restrict the use of Agrif to the forward case only
      IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() )   CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )
#endif
      !
      IF(lwp) WRITE(numout,*)    '     Time filter choice, nn_bt_flt: ', nn_bt_flt
      SELECT CASE ( nn_bt_flt )
         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '           Dirac'
         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = nn_e'
         CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '           Boxcar: width = 2*nn_e' 
         CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' )
      END SELECT
      !
      IF(lwp) WRITE(numout,*) ' '
      IF(lwp) WRITE(numout,*) '     nn_e = ', nn_e
      IF(lwp) WRITE(numout,*) '     Barotropic time step [s] is :', rDt_e
      IF(lwp) WRITE(numout,*) '     Maximum Courant number is   :', zcmax
      !
      IF(lwp) WRITE(numout,*)    '     Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha
      IF ((ln_bt_av.AND.nn_bt_flt/=0).AND.(rn_bt_alpha>0._wp)) THEN
         CALL ctl_stop( 'dynspg_ts ERROR: if rn_bt_alpha > 0, remove temporal averaging' )
      ENDIF
      !
      IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN
         CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' )
      ENDIF
      IF( zcmax>0.9_wp ) THEN
         CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' )          
      ENDIF
      !
      !                             ! Allocate time-splitting arrays
      IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
      !
      !                             ! read restart when needed
      CALL ts_rst( nit000, 'READ' )
      !
   END SUBROUTINE dyn_spg_ts_init

   
   SUBROUTINE dyn_cor_2D_init( Kmm )
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE dyn_cor_2D_init  ***
      !!
      !! ** Purpose : Set time splitting options
      !! Set arrays to remove/compute coriolis trend.
      !! Do it once during initialization if volume is fixed, else at each long time step.
      !! Note that these arrays are also used during barotropic loop. These are however frozen
      !! although they should be updated in the variable volume case. Not a big approximation.
      !! To remove this approximation, copy lines below inside barotropic loop
      !! and update depths at T- points (ht) at each barotropic time step
      !!
      !! Compute zwz = f / ( height of the water colomn )
      !!----------------------------------------------------------------------
      INTEGER,  INTENT(in)         ::  Kmm  ! Time index
      INTEGER  ::   ji ,jj, jk              ! dummy loop indices
      REAL(wp) ::   z1_ht
      !!----------------------------------------------------------------------
      !
      SELECT CASE( nvor_scheme )
      CASE( np_EEN, np_ENE, np_ENS , np_MIX )   !=  schemes using the same e3f definition
         SELECT CASE( nn_e3f_typ )                  !* ff_f/e3 at F-point
         CASE ( 0 )                                   ! original formulation  (masked averaging of e3t divided by 4)
            DO_2D( 0, 0, 0, 0 )
               zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1)   &
                    &       + ht(ji,jj  ) + ht(ji+1,jj  ) ) * 0.25_wp  
               IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
            END_2D
         CASE ( 1 )                                   ! new formulation  (masked averaging of e3t divided by the sum of mask)
            DO_2D( 0, 0, 0, 0 )
               zwz(ji,jj) =     (    ht(ji,jj+1) +     ht(ji+1,jj+1)      &
                    &            +   ht(ji,jj  ) +     ht(ji+1,jj  )  )   &
                    &    / ( MAX(ssmask(ji,jj+1) + ssmask(ji+1,jj+1)      &
                    &          + ssmask(ji,jj  ) + ssmask(ji+1,jj  ) , 1._wp  )   )
               IF( zwz(ji,jj) /= 0._wp )   zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)
            END_2D
         END SELECT
         CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )
      END SELECT
      !
      SELECT CASE( nvor_scheme )
      CASE( np_EEN )
         !
         ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp
         DO_2D( 0, 1, 0, 1 )
            ftne(ji,jj) = zwz(ji-1,jj  ) + zwz(ji  ,jj  ) + zwz(ji  ,jj-1)
            ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj  ) + zwz(ji  ,jj  )
            ftse(ji,jj) = zwz(ji  ,jj  ) + zwz(ji  ,jj-1) + zwz(ji-1,jj-1)
            ftsw(ji,jj) = zwz(ji  ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj  )
         END_2D
         !
      CASE( np_EET )                            != EEN scheme using e3t energy conserving scheme
         ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp
         DO_2D( 0, 1, 0, 1 )
            z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) )
            ftne(ji,jj) = ( ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) ) * z1_ht
            ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) + ff_f(ji  ,jj  ) ) * z1_ht
            ftse(ji,jj) = ( ff_f(ji  ,jj  ) + ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht
            ftsw(ji,jj) = ( ff_f(ji  ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj  ) ) * z1_ht
         END_2D
         !
      END SELECT
      
   END SUBROUTINE dyn_cor_2d_init


   SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV,    zu_trd, zv_trd   )
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE dyn_cor_2d  ***
      !!
      !! ** Purpose : Compute u and v coriolis trends
      !!----------------------------------------------------------------------
      INTEGER  ::   ji ,jj                             ! dummy loop indices
      REAL(wp) ::   zx1, zx2, zy1, zy2, z1_hu, z1_hv   !   -      -
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   )  :: pht, phu, phv, punb, pvnb, zhV
      REAL(dp), DIMENSION(jpi,jpj), INTENT(in   )  :: zhU
Guillaume Samson's avatar
Guillaume Samson committed
      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: zu_trd, zv_trd
      !!----------------------------------------------------------------------
      SELECT CASE( nvor_scheme )
      CASE( np_ENT )                ! enstrophy conserving scheme (f-point)
         DO_2D( 0, 0, 0, 0 )
            z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) )
            z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) )
            zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu                    &
               &               * (  e1e2t(ji+1,jj)*pht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) )   &
               &                  + e1e2t(ji  ,jj)*pht(ji  ,jj)*ff_t(ji  ,jj) * ( pvnb(ji  ,jj) + pvnb(ji  ,jj-1) )   )
               !
            zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv                    &
               &               * (  e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) )   & 
               &                  + e1e2t(ji,jj  )*pht(ji,jj  )*ff_t(ji,jj  ) * ( punb(ji,jj  ) + punb(ji-1,jj  ) )   ) 
         END_2D
         !         
      CASE( np_ENE , np_MIX )        ! energy conserving scheme (t-point) ENE or MIX
         DO_2D( 0, 0, 0, 0 )
            zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj)
            zy2 = ( zhV(ji,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
            zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj)
            zx2 = ( zhU(ji  ,jj) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
            ! energy conserving formulation for planetary vorticity term
            zu_trd(ji,jj) =   r1_4 * ( zwz(ji  ,jj-1) * zy1 + zwz(ji,jj) * zy2 )
            zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj  ) * zx1 + zwz(ji,jj) * zx2 )
         END_2D
         !
      CASE( np_ENS )                ! enstrophy conserving scheme (f-point)
         DO_2D( 0, 0, 0, 0 )
            zy1 =   r1_8 * ( zhV(ji  ,jj-1) + zhV(ji+1,jj-1) &
              &            + zhV(ji  ,jj  ) + zhV(ji+1,jj  ) ) * r1_e1u(ji,jj)
            zx1 = - r1_8 * ( zhU(ji-1,jj  ) + zhU(ji-1,jj+1) &
              &            + zhU(ji  ,jj  ) + zhU(ji  ,jj+1) ) * r1_e2v(ji,jj)
            zu_trd(ji,jj)  = zy1 * ( zwz(ji  ,jj-1) + zwz(ji,jj) )
            zv_trd(ji,jj)  = zx1 * ( zwz(ji-1,jj  ) + zwz(ji,jj) )
         END_2D
         !
      CASE( np_EET , np_EEN )      ! energy & enstrophy scheme (using e3t or e3f)         
         DO_2D( 0, 0, 0, 0 )
            zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * (  ftne(ji,jj  ) * zhV(ji  ,jj  ) &
             &                                         + ftnw(ji+1,jj) * zhV(ji+1,jj  ) &
             &                                         + ftse(ji,jj  ) * zhV(ji  ,jj-1) &
             &                                         + ftsw(ji+1,jj) * zhV(ji+1,jj-1) )
            zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * (  ftsw(ji,jj+1) * zhU(ji-1,jj+1) &
             &                                         + ftse(ji,jj+1) * zhU(ji  ,jj+1) &
             &                                         + ftnw(ji,jj  ) * zhU(ji-1,jj  ) &
             &                                         + ftne(ji,jj  ) * zhU(ji  ,jj  ) )
         END_2D
         !
      END SELECT
      !
   END SUBROUTINE dyn_cor_2D


   SUBROUTINE wad_tmsk( pssh, ptmsk )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE wad_lmt  ***
      !!                    
      !! ** Purpose :   set wetting & drying mask at tracer points 
      !!              for the current barotropic sub-step 
      !!
      !! ** Method  :   ??? 
      !!
      !! ** Action  :  ptmsk : wetting & drying t-mask
      !!----------------------------------------------------------------------
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pssh    !
      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   ptmsk   !
      !
      INTEGER  ::   ji, jj   ! dummy loop indices
      !!----------------------------------------------------------------------
      !
      IF( ln_wd_dl_rmp ) THEN     
         DO_2D( 1, 1, 1, 1 )
            IF    ( pssh(ji,jj) + ht_0(ji,jj) >  2._wp * rn_wdmin1 ) THEN 
               !           IF    ( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin2 ) THEN 
               ptmsk(ji,jj) = 1._wp
            ELSEIF( pssh(ji,jj) + ht_0(ji,jj) >          rn_wdmin1 ) THEN
               ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) -  rn_wdmin1 )*r_rn_wdmin1) )
            ELSE 
               ptmsk(ji,jj) = 0._wp
            ENDIF
         END_2D
      ELSE  
         DO_2D( 1, 1, 1, 1 )
            IF ( pssh(ji,jj) + ht_0(ji,jj) >  rn_wdmin1 ) THEN   ;   ptmsk(ji,jj) = 1._wp
            ELSE                                                 ;   ptmsk(ji,jj) = 0._wp
            ENDIF
         END_2D
      ENDIF
      !
   END SUBROUTINE wad_tmsk


   SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE wad_lmt  ***
      !!                    
      !! ** Purpose :   set wetting & drying mask at tracer points 
      !!              for the current barotropic sub-step 
      !!
      !! ** Method  :   ??? 
      !!
      !! ** Action  :  ptmsk : wetting & drying t-mask
      !!----------------------------------------------------------------------
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pTmsk              ! W & D t-mask
      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout)  :: phV, pu, pv! ocean velocities and transports
      REAL(dp), DIMENSION(jpi,jpj), INTENT(inout)  :: phU! ocean velocities and transports
Guillaume Samson's avatar
Guillaume Samson committed
      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pUmsk, pVmsk       ! W & D u- and v-mask
      !
      INTEGER  ::   ji, jj   ! dummy loop indices
      !!----------------------------------------------------------------------
      !
      DO_2D( 1, 0, 1, 1 )   ! not jpi-column
         IF ( phU(ji,jj) > 0._wp ) THEN   ;   pUmsk(ji,jj) = pTmsk(ji  ,jj) 
         ELSE                             ;   pUmsk(ji,jj) = pTmsk(ji+1,jj)  
         ENDIF
         phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj)
         pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj)
      END_2D
      !
      DO_2D( 1, 1, 1, 0 )   ! not jpj-row
         IF ( phV(ji,jj) > 0._wp ) THEN   ;   pVmsk(ji,jj) = pTmsk(ji,jj  )
         ELSE                             ;   pVmsk(ji,jj) = pTmsk(ji,jj+1)  
         ENDIF
         phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 
         pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj)
      END_2D
      !
   END SUBROUTINE wad_Umsk


   SUBROUTINE wad_spg( pshn, zcpx, zcpy )
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE  wad_sp  ***
      !!
      !! ** Purpose : 
      !!----------------------------------------------------------------------
      INTEGER  ::   ji ,jj               ! dummy loop indices
      LOGICAL  ::   ll_tmp1, ll_tmp2
      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: pshn
      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
      !!----------------------------------------------------------------------
      DO_2D( 0, 0, 0, 0 )
         ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji+1,jj) ) >                &
              &      MAX( -ht_0(ji,jj)               , -ht_0(ji+1,jj) ) .AND.            &
              &      MAX(  pshn(ji,jj) + ht_0(ji,jj) ,  pshn(ji+1,jj) + ht_0(ji+1,jj) )  &
              &                                                         > rn_wdmin1 + rn_wdmin2
         ll_tmp2 = ( ABS( pshn(ji+1,jj)            -  pshn(ji  ,jj))  > 1.E-12 ).AND.( &
              &      MAX(   pshn(ji,jj)              ,  pshn(ji+1,jj) ) >                &
              &      MAX(  -ht_0(ji,jj)              , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
         IF(ll_tmp1) THEN
            zcpx(ji,jj) = 1.0_wp
         ELSEIF(ll_tmp2) THEN
            ! no worries about  pshn(ji+1,jj) -  pshn(ji  ,jj) = 0, it won't happen ! here
            zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) &
                 &           / (pshn(ji+1,jj) - pshn(ji  ,jj)) )
            zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)
         ELSE
            zcpx(ji,jj) = 0._wp
         ENDIF
         !
         ll_tmp1 = MIN(  pshn(ji,jj)               ,  pshn(ji,jj+1) ) >                &
              &      MAX( -ht_0(ji,jj)               , -ht_0(ji,jj+1) ) .AND.            &
              &      MAX(  pshn(ji,jj) + ht_0(ji,jj) ,  pshn(ji,jj+1) + ht_0(ji,jj+1) )  &
              &                                                       > rn_wdmin1 + rn_wdmin2
         ll_tmp2 = ( ABS( pshn(ji,jj)              -  pshn(ji,jj+1))  > 1.E-12 ).AND.( &
              &      MAX(   pshn(ji,jj)              ,  pshn(ji,jj+1) ) >                &
              &      MAX(  -ht_0(ji,jj)              , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
         
         IF(ll_tmp1) THEN
            zcpy(ji,jj) = 1.0_wp
         ELSE IF(ll_tmp2) THEN
            ! no worries about  pshn(ji,jj+1) -  pshn(ji,jj  ) = 0, it won't happen ! here
            zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) &
                 &           / (pshn(ji,jj+1) - pshn(ji,jj  )) )
            zcpy(ji,jj) = MAX(  0._wp , MIN( zcpy(ji,jj) , 1.0_wp )  )
         ELSE
            zcpy(ji,jj) = 0._wp
         ENDIF
      END_2D
            
   END SUBROUTINE wad_spg
     

   SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dyn_drg_init  ***
      !!                    
      !! ** Purpose : - add the baroclinic top/bottom drag contribution to 
      !!              the baroclinic part of the barotropic RHS
      !!              - compute the barotropic drag coefficients
      !!
      !! ** Method  :   computation done over the INNER domain only 
      !!----------------------------------------------------------------------
      INTEGER                             , INTENT(in   ) ::  Kbb, Kmm           ! ocean time level indices
      REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in   ) ::  puu, pvv           ! ocean velocities and RHS of momentum equation
Guillaume Samson's avatar
Guillaume Samson committed
      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(in   ) ::  puu_b, pvv_b       ! barotropic velocities at main time levels
      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(inout) ::  pu_RHSi, pv_RHSi   ! baroclinic part of the barotropic RHS
      REAL(wp), DIMENSION(jpi,jpj)        , INTENT(  out) ::  pCdU_u , pCdU_v    ! barotropic drag coefficients
      !
      INTEGER  ::   ji, jj   ! dummy loop indices
      INTEGER  ::   ikbu, ikbv, iktu, iktv
      REAL(wp) ::   zztmp
      REAL(wp), DIMENSION(jpi,jpj) ::   zu_i, zv_i
      !!----------------------------------------------------------------------
      !
      !                    !==  Set the barotropic drag coef.  ==!
      !
      IF( ln_isfcav.OR.ln_drgice_imp ) THEN          ! top+bottom friction (ocean cavities)
         
         DO_2D( 0, 0, 0, 0 )
            pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )
            pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )
         END_2D
      ELSE                          ! bottom friction only
         DO_2D( 0, 0, 0, 0 )
            pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
            pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
         END_2D
      ENDIF
      !
      !                    !==  BOTTOM stress contribution from baroclinic velocities  ==!
      !
      IF( ln_bt_fw ) THEN                 ! FORWARD integration: use NOW bottom baroclinic velocities
         
         DO_2D( 0, 0, 0, 0 )
            ikbu = mbku(ji,jj)       
            ikbv = mbkv(ji,jj)    
            zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
            zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
         END_2D
      ELSE                                ! CENTRED integration: use BEFORE bottom baroclinic velocities
         
         DO_2D( 0, 0, 0, 0 )
            ikbu = mbku(ji,jj)       
            ikbv = mbkv(ji,jj)    
            zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb)
            zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb)
         END_2D
      ENDIF
      !
      IF( ln_wd_il ) THEN      ! W/D : use the "clipped" bottom friction   !!gm   explain WHY, please !
         zztmp = -1._wp / rDt_e
         DO_2D( 0, 0, 0, 0 )
            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) *  wdrampu(ji,jj) * MAX(                                 & 
                 &                              r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp  )
            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) *  wdrampv(ji,jj) * MAX(                                 & 
                 &                              r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp  )
         END_2D
      ELSE                    ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
         
         DO_2D( 0, 0, 0, 0 )
            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj)
            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj)
         END_2D
      END IF
      !
      !                    !==  TOP stress contribution from baroclinic velocities  ==!   (no W/D case)
      !
      IF( ln_isfcav.OR.ln_drgice_imp ) THEN
         !
         IF( ln_bt_fw ) THEN                ! FORWARD integration: use NOW top baroclinic velocity
            
            DO_2D( 0, 0, 0, 0 )
               iktu = miku(ji,jj)
               iktv = mikv(ji,jj)
               zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
               zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
            END_2D
         ELSE                                ! CENTRED integration: use BEFORE top baroclinic velocity
            
            DO_2D( 0, 0, 0, 0 )
               iktu = miku(ji,jj)
               iktv = mikv(ji,jj)
               zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb)
               zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
            END_2D
         ENDIF
         !
         !                    ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
         
         DO_2D( 0, 0, 0, 0 )
            pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj)
            pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj)
         END_2D
         !
      ENDIF
      !
   END SUBROUTINE dyn_drg_init

   SUBROUTINE ts_bck_interp( jn, ll_init,       &   ! <== in
      &                      za0, za1, za2, za3 )   ! ==> out
      !!----------------------------------------------------------------------
      INTEGER ,INTENT(in   ) ::   jn                   ! index of sub time step
      LOGICAL ,INTENT(in   ) ::   ll_init              !
      REAL(wp),INTENT(  out)  :: za0, za2, za3! Half-step back interpolation coefficient
      REAL(dp),INTENT(  out)  :: za1! Half-step back interpolation coefficient
Guillaume Samson's avatar
Guillaume Samson committed
      !
      REAL(wp) ::   zepsilon, zgamma                   !   -      -
      !!----------------------------------------------------------------------
      !                             ! set Half-step back interpolation coefficient
      IF    ( jn==1 .AND. ll_init ) THEN   !* Forward-backward
         za0 = 1._wp                        
         za1 = 0._wp                           
         za2 = 0._wp
         za3 = 0._wp
      ELSEIF( jn==2 .AND. ll_init ) THEN   !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12
         za0 = 1.0833333333333_wp                 ! za0 = 1-gam-eps
         za1 =-0.1666666666666_wp                 ! za1 = gam
         za2 = 0.0833333333333_wp                 ! za2 = eps
         za3 = 0._wp              
      ELSE                                 !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 
         IF( rn_bt_alpha == 0._wp ) THEN      ! Time diffusion  
            za0 = 0.614_wp                        ! za0 = 1/2 +   gam + 2*eps
            za1 = 0.285_wp                        ! za1 = 1/2 - 2*gam - 3*eps
            za2 = 0.088_wp                        ! za2 = gam
            za3 = 0.013_wp                        ! za3 = eps
         ELSE                                 ! no time diffusion
            zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha
            zgamma   = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha
            za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon
            za1 = 1._wp - za0 - zgamma - zepsilon
            za2 = zgamma
            za3 = zepsilon
         ENDIF 
      ENDIF
   END SUBROUTINE ts_bck_interp


   !!======================================================================
END MODULE dynspg_ts