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 5038 for branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 – NEMO

Ignore:
Timestamp:
2015-01-20T15:26:13+01:00 (9 years ago)
Author:
jamesharle
Message:

Merging branch with HEAD of the trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90

    r4792 r5038  
    1212   !!             -   ! 2010-06  (C. Ethe, G. Madec)  merge TRA-TRC 
    1313   !!             -   ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level 
     14   !!             -   ! 2013-04  (F. Roquet, G. Madec)  use of eosbn2 instead of local hard coded alpha and beta 
    1415   !!---------------------------------------------------------------------- 
    1516#if   defined key_trabbl   ||   defined key_esopa 
     
    2829   USE phycst         ! physical constant 
    2930   USE eosbn2         ! equation of state 
    30    USE trdmod_oce     ! trends: ocean variables 
     31   USE trd_oce     ! trends: ocean variables 
    3132   USE trdtra         ! trends: active tracers 
    32    USE iom            ! IOM server 
     33   ! 
     34   USE iom            ! IOM library                
    3335   USE in_out_manager ! I/O manager 
    3436   USE lbclnk         ! ocean lateral boundary conditions 
     
    3638   USE wrk_nemo       ! Memory Allocation 
    3739   USE timing         ! Timing 
    38  
     40   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    3941 
    4042   IMPLICIT NONE 
     
    5759   REAL(wp), PUBLIC ::   rn_gambbl   !: lateral coeff. for bottom boundary layer scheme [s] 
    5860 
    59    LOGICAL , PUBLIC ::   l_bbl                  !: flag to compute bbl diffu. flux coef and transport 
     61   LOGICAL , PUBLIC ::   l_bbl       !: flag to compute bbl diffu. flux coef and transport 
    6062 
    6163   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer 
     
    8486         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     & 
    8587         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          & 
    86          &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) , STAT= tra_bbl_alloc                ) 
     88         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc ) 
    8789         ! 
    8890      IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc ) 
     
    104106      !!---------------------------------------------------------------------- 
    105107      INTEGER, INTENT( in ) ::   kt   ! ocean time-step 
    106       !! 
     108      ! 
    107109      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds 
    108110      !!---------------------------------------------------------------------- 
     
    110112      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl') 
    111113      ! 
    112       IF( l_trdtra )   THEN                        !* Save ta and sa trends 
     114      IF( l_trdtra )   THEN                         !* Save ta and sa trends 
    113115         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    114116         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) 
     
    116118      ENDIF 
    117119 
    118       IF( l_bbl )  CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl) 
    119  
    120       IF( nn_bbl_ldf == 1 ) THEN                   !* Diffusive bbl 
     120      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl) 
     121 
     122      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl 
    121123         ! 
    122124         CALL tra_bbl_dif( tsb, tsa, jpts ) 
    123125         IF( ln_ctl )  & 
    124126         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, & 
    125          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     127            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    126128         ! lateral boundary conditions ; just need for outputs 
    127129         CALL lbc_lnk( ahu_bbl, 'U', 1. )     ;     CALL lbc_lnk( ahv_bbl, 'V', 1. ) 
     
    131133      END IF 
    132134 
    133       IF( nn_bbl_adv /= 0 ) THEN                !* Advective bbl 
     135      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl 
    134136         ! 
    135137         CALL tra_bbl_adv( tsb, tsa, jpts ) 
    136138         IF(ln_ctl)   & 
    137139         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   & 
    138          &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
     140            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' ) 
    139141         ! lateral boundary conditions ; just need for outputs 
    140142         CALL lbc_lnk( utr_bbl, 'U', 1. )     ;   CALL lbc_lnk( vtr_bbl, 'V', 1. ) 
     
    147149         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:) 
    148150         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:) 
    149          CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbl, ztrdt ) 
    150          CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_bbl, ztrds ) 
     151         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt ) 
     152         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds ) 
    151153         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds ) 
    152154      ENDIF 
     
    164166      !!                advection terms. 
    165167      !! 
    166       !! ** Method  : 
    167       !!        * diffusive bbl (nn_bbl_ldf=1) : 
     168      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) : 
    168169      !!        When the product grad( rho) * grad(h) < 0 (where grad is an 
    169170      !!      along bottom slope gradient) an additional lateral 2nd order 
     
    179180      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 
    180181      !!---------------------------------------------------------------------- 
    181       ! 
    182182      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers 
    183183      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields 
     
    196196      DO jn = 1, kjpt                                     ! tracer loop 
    197197         !                                                ! =========== 
    198 #  if defined key_vectopt_loop 
    199          DO jj = 1, 1   ! vector opt. (forced unrolling) 
    200             DO ji = 1, jpij 
    201 #else 
    202198         DO jj = 1, jpj 
    203199            DO ji = 1, jpi 
    204 #endif 
    205                ik = mbkt(ji,jj)                        ! bottom T-level index 
    206                zptb(ji,jj) = ptb(ji,jj,ik,jn)              ! bottom before T and S 
     200               ik = mbkt(ji,jj)                              ! bottom T-level index 
     201               zptb(ji,jj) = ptb(ji,jj,ik,jn)       ! bottom before T and S 
    207202            END DO 
    208203         END DO 
    209          !                                                ! Compute the trend 
    210 #  if defined key_vectopt_loop 
    211          DO jj = 1, 1   ! vector opt. (forced unrolling) 
    212             DO ji = jpi+1, jpij-jpi-1 
    213 #  else 
    214          DO jj = 2, jpjm1 
     204         !                
     205         DO jj = 2, jpjm1                                    ! Compute the trend 
    215206            DO ji = 2, jpim1 
    216 #  endif 
    217                ik = mbkt(ji,jj)                            ! bottom T-level index 
     207               ik = mbkt(ji,jj)                              ! bottom T-level index 
    218208               zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik) 
    219209               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         & 
     
    264254      DO jn = 1, kjpt                                            ! tracer loop 
    265255         !                                                       ! =========== 
    266 # if defined key_vectopt_loop 
    267          DO jj = 1, 1 
    268             DO ji = 1, jpij-jpi-1   ! vector opt. (forced unrolling) 
    269 # else 
    270256         DO jj = 1, jpjm1 
    271257            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west 
    272 # endif 
    273258               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection 
    274259                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf) 
     
    333318      !!                advection terms. 
    334319      !! 
    335       !! ** Method  : 
    336       !!        * diffusive bbl (nn_bbl_ldf=1) : 
     320      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) : 
    337321      !!        When the product grad( rho) * grad(h) < 0 (where grad is an 
    338322      !!      along bottom slope gradient) an additional lateral 2nd order 
     
    342326      !!      a downslope velocity of 20 cm/s if the condition for slope 
    343327      !!      convection is satified) 
    344       !!        * advective bbl (nn_bbl_adv=1 or 2) : 
     328      !!              * advective bbl (nn_bbl_adv=1 or 2) : 
    345329      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity 
    346330      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation 
     
    353337      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430. 
    354338      !!---------------------------------------------------------------------- 
    355       ! 
    356339      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index 
    357       INTEGER         , INTENT(in   ) ::   kit000          ! first time step index 
     340      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index 
    358341      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator) 
    359342      !! 
    360343      INTEGER  ::   ji, jj                    ! dummy loop indices 
    361344      INTEGER  ::   ik                        ! local integers 
    362       INTEGER  ::   iis , iid , ijs , ijd     !   -       - 
    363       INTEGER  ::   ikus, ikud, ikvs, ikvd    !   -       - 
    364       REAL(wp) ::   zsign, zsigna, zgbbl      ! local scalars 
    365       REAL(wp) ::   zgdrho, zt, zs, zh        !   -      - 
    366       !! 
    367       REAL(wp) ::   fsalbt, fsbeta, pft, pfs, pfh   ! statement function 
    368       REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb, ztb, zsb, zdep 
    369       !!----------------------- zv_bbl----------------------------------------------- 
    370       ! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients 
    371       ! ================            pft :  potential temperature in degrees celcius 
    372       !                             pfs :  salinity anomaly (s-35) in psu 
    373       !                             pfh :  depth in meters 
    374       ! nn_eos = 0  (Jackett and McDougall 1994 formulation) 
    375       fsalbt( pft, pfs, pfh ) =                                              &   ! alpha/beta 
    376          ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    & 
    377                                    - 0.203814e-03 ) * pft                    & 
    378                                    + 0.170907e-01 ) * pft                    & 
    379                                    + 0.665157e-01                            & 
    380          +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   & 
    381          +  ( ( - 0.302285e-13 * pfh                                         & 
    382                 - 0.251520e-11 * pfs                                         & 
    383                 + 0.512857e-12 * pft * pft          ) * pfh                  & 
    384                                      - 0.164759e-06   * pfs                  & 
    385              +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  & 
    386                                      + 0.380374e-04 ) * pfh 
    387       fsbeta( pft, pfs, pfh ) =                                              &   ! beta 
    388          ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft                      & 
    389                                  - 0.301985e-05 ) * pft                      & 
    390                                  + 0.785567e-03                              & 
    391          + (     0.515032e-08 * pfs                                          & 
    392                + 0.788212e-08 * pft - 0.356603e-06 ) * pfs                   & 
    393                +(  (   0.121551e-17 * pfh                                    & 
    394                      - 0.602281e-15 * pfs                                    & 
    395                      - 0.175379e-14 * pft + 0.176621e-12 ) * pfh             & 
    396                                           + 0.408195e-10   * pfs             & 
    397                  + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft             & 
    398                                           - 0.121555e-07 ) * pfh 
    399       !!---------------------------------------------------------------------- 
    400  
     345      INTEGER  ::   iis, iid, ikus, ikud      !   -       - 
     346      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       - 
     347      REAL(wp) ::   za, zb, zgdrho            ! local scalars 
     348      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      - 
     349      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace 
     350      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace 
     351      !!---------------------------------------------------------------------- 
    401352      ! 
    402353      IF( nn_timing == 1 )  CALL timing_start( 'bbl') 
    403354      ! 
    404       CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 
    405       ! 
    406  
    407355      IF( kt == kit000 )  THEN 
    408356         IF(lwp)  WRITE(numout,*) 
     
    410358         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~' 
    411359      ENDIF 
    412  
    413       !                                        !* bottom temperature, salinity, velocity and depth 
    414 #if defined key_vectopt_loop 
    415       DO jj = 1, 1   ! vector opt. (forced unrolling) 
    416          DO ji = 1, jpij 
    417 #else 
     360      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity) 
    418361      DO jj = 1, jpj 
    419362         DO ji = 1, jpi 
    420 #endif 
    421             ik = mbkt(ji,jj)                        ! bottom T-level index 
    422             ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S 
    423             zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1) 
    424             zdep(ji,jj) = gdept_0(ji,jj,ik)         ! bottom T-level reference depth 
     363            ik = mbkt(ji,jj)                             ! bottom T-level index 
     364            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S 
     365            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal) 
    425366            ! 
    426             zub(ji,jj) = un(ji,jj,mbku(ji,jj))      ! bottom velocity 
    427             zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) 
     367            zdep(ji,jj) = fsdept(ji,jj,ik)               ! bottom T-level reference depth 
     368            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity 
     369            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj)) 
    428370         END DO 
    429371      END DO 
    430  
     372      ! 
     373      CALL eos_rab( zts, zdep, zab ) 
     374      ! 
    431375      !                                   !-------------------! 
    432376      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   ! 
    433377         !                                !-------------------! 
    434378         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 ) 
    435             DO ji = 1, jpim1 
    436                !                                                ! i-direction 
    437                zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth 
    438                zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 
    439                zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 
    440                !                                                         ! masked bbl i-gradient of density 
    441                zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    & 
    442                   &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1) 
     379            DO ji = 1, fs_jpim1   ! vector opt. 
     380               !                                                   ! i-direction 
     381               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point 
     382               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 
     383               !                                                         ! 2*masked bottom density gradient 
     384               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    & 
     385                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1) 
    443386               ! 
    444                zsign          = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope ) 
    445                ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)                  ! masked diffusive flux coeff. 
     387               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope ) 
     388               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff. 
    446389               ! 
    447                !                                                ! j-direction 
    448                zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                ! T, S anomalie, and depth 
    449                zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0 
    450                zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) 
    451                !                                                         ! masked bbl j-gradient of density 
    452                zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    & 
    453                   &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1) 
     390               !                                                   ! j-direction 
     391               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point 
     392               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) 
     393               !                                                         ! 2*masked bottom density gradient 
     394               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    & 
     395                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1) 
    454396               ! 
    455                zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope ) 
     397               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope ) 
    456398               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj) 
    457                ! 
    458399            END DO 
    459400         END DO 
     
    469410            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0 
    470411               DO ji = 1, fs_jpim1   ! vector opt. 
    471                   !                                               ! i-direction 
    472                   zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth 
    473                   zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 
    474                   zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 
    475                   !                                                           ! masked bbl i-gradient of density 
    476                   zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    & 
    477                      &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1) 
    478                   ! 
    479                   zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope 
    480                   zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope 
    481                   ! 
    482                   !                                                           ! bbl velocity 
     412                  !                                                  ! i-direction 
     413                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point 
     414                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 
     415                  !                                                          ! 2*masked bottom density gradient  
     416                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    & 
     417                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1) 
     418                  ! 
     419                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope 
     420                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope 
     421                  ! 
     422                  !                                                          ! bbl velocity 
    483423                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj) 
    484424                  ! 
    485                   !                                               ! j-direction 
    486                   zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth 
    487                   zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0 
    488                   zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) ) 
    489                   !                                                           ! masked bbl j-gradient of density 
    490                   zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    & 
    491                      &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1) 
    492                   zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )    ! sign of j-gradient * j-slope 
    493                   zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )    ! sign of u * i-slope 
    494                   ! 
    495                   !                                                           ! bbl velocity 
     425                  !                                                  ! j-direction 
     426                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point 
     427                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) 
     428                  !                                                          ! 2*masked bottom density gradient 
     429                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    & 
     430                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1) 
     431                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope 
     432                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope 
     433                  ! 
     434                  !                                                          ! bbl transport 
    496435                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj) 
    497436               END DO 
     
    502441            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down 
    503442               DO ji = 1, fs_jpim1   ! vector opt. 
    504                   !                                         ! i-direction 
     443                  !                                                  ! i-direction 
    505444                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf) 
    506                   iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) ) 
    507                   ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj) 
    508                   ! 
    509                   !                                             ! mid-depth density anomalie (up-slope minus down-slope) 
    510                   zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth 
    511                   zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0 
    512                   zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) ) 
    513                   zgdrho =    fsbeta( zt, zs, zh )                                    & 
    514                      &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    & 
    515                      &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask(ji,jj,1) 
    516                   zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep 
    517                   ! 
    518                   !                                             ! bbl transport (down-slope direction) 
     445                  iid  = ji + MAX( 0, mgrhu(ji,jj) ) 
     446                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) ) 
     447                  ! 
     448                  ikud = mbku_d(ji,jj) 
     449                  ikus = mbku(ji,jj) 
     450                  ! 
     451                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point 
     452                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal) 
     453                  !                                                          !   masked bottom density gradient 
     454                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    & 
     455                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1) 
     456                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep 
     457                  ! 
     458                  !                                                          ! bbl transport (down-slope direction) 
    519459                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) ) 
    520460                  ! 
    521                   !                                         ! j-direction 
     461                  !                                                  ! j-direction 
    522462                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf) 
    523                   ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) ) 
    524                   ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj) 
    525                   ! 
    526                   !                                             ! mid-depth density anomalie (up-slope minus down-slope) 
    527                   zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth 
    528                   zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0 
    529                   zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) ) 
    530                   zgdrho =    fsbeta( zt, zs, zh )                                    & 
    531                      &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    & 
    532                      &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask(ji,jj,1) 
    533                   zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep 
    534                   ! 
    535                   !                                             ! bbl transport (down-slope direction) 
     463                  ijd  = jj + MAX( 0, mgrhv(ji,jj) ) 
     464                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) ) 
     465                  ! 
     466                  ikvd = mbkv_d(ji,jj) 
     467                  ikvs = mbkv(ji,jj) 
     468                  ! 
     469                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point 
     470                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal) 
     471                  !                                                          !   masked bottom density gradient 
     472                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    & 
     473                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1) 
     474                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep 
     475                  ! 
     476                  !                                                          ! bbl transport (down-slope direction) 
    536477                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) ) 
    537478               END DO 
     
    541482      ENDIF 
    542483      ! 
    543       CALL wrk_dealloc( jpi, jpj, zub, zvb, ztb, zsb, zdep ) 
    544       ! 
    545484      IF( nn_timing == 1 )  CALL timing_stop( 'bbl') 
    546485      ! 
     
    558497      !!---------------------------------------------------------------------- 
    559498      INTEGER ::   ji, jj               ! dummy loop indices 
    560       INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integer 
    561       INTEGER  ::   ios                 ! Local integer output status for namelist read 
     499      INTEGER ::   ii0, ii1, ij0, ij1   ! local integer 
     500      INTEGER ::   ios                  !   -      - 
    562501      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk 
    563502      !! 
     
    598537      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)' 
    599538 
    600       IF( nn_eos /= 0 )   CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' ) 
    601  
    602539      !                             !* vertical index of  "deep" bottom u- and v-points 
    603540      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv) 
     
    607544         END DO 
    608545      END DO 
    609       ! convert into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
     546      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 
    610547      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    611548      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 ) 
    612549 
    613                                      !* sign of grad(H) at u- and v-points 
    614       mgrhu(jpi,:) = 0.    ;    mgrhu(:,jpj) = 0.   ;    mgrhv(jpi,:) = 0.    ;    mgrhv(:,jpj) = 0. 
     550                                        !* sign of grad(H) at u- and v-points 
     551      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0 
    615552      DO jj = 1, jpjm1 
    616553         DO ji = 1, jpim1 
     
    621558 
    622559      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point 
    623          DO ji = 1, jpim1           ! minimum of top & bottom e3u_0 (e3v_0) 
     560         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0) 
    624561            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) ) 
    625562            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) ) 
Note: See TracChangeset for help on using the changeset viewer.