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 6225 for branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90 – NEMO

Ignore:
Timestamp:
2016-01-08T10:35:19+01:00 (8 years ago)
Author:
jamesharle
Message:

Update MPP_BDY_UPDATE branch to be consistent with head of trunk

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90

    r3294 r6225  
    88   !!             -   !  2004-03  (C. Ethe)  adapted for passive tracers 
    99   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA  
     10   !!            3.6  !  2014-11  (P. Mathiot) Add zps_hde_isf (needed to open a cavity) 
    1011   !!====================================================================== 
    1112    
     
    2728   PRIVATE 
    2829 
    29    PUBLIC   zps_hde    ! routine called by step.F90 
     30   PUBLIC   zps_hde     ! routine called by step.F90 
     31   PUBLIC   zps_hde_isf ! routine called by step.F90 
    3032 
    3133   !! * Substitutions 
    32 #  include "domzgr_substitute.h90" 
    3334#  include "vectopt_loop_substitute.h90" 
    3435   !!---------------------------------------------------------------------- 
     
    4041 
    4142   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   & 
    42                                  prd, pgru, pgrv    ) 
     43      &                          prd, pgru, pgrv    ) 
    4344      !!---------------------------------------------------------------------- 
    4445      !!                     ***  ROUTINE zps_hde  *** 
     
    7475      !!          Idem for di(s) and dj(s)           
    7576      !! 
    76       !!      For rho, we call eos_insitu_2d which will compute rd~(t~,s~) at  
    77       !!      the good depth zh from interpolated T and S for the different 
    78       !!      formulation of the equation of state (eos). 
     77      !!      For rho, we call eos which will compute rd~(t~,s~) at the right 
     78      !!      depth zh from interpolated T and S for the different formulations 
     79      !!      of the equation of state (eos). 
    7980      !!      Gradient formulation for rho : 
    80       !!          di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~ 
    81       !! 
    82       !! ** Action  : - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 
    83       !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points  
    84       !!---------------------------------------------------------------------- 
    85       ! 
     81      !!          di(rho) = rd~ - rd(i,j,k)   or  rd(i+1,j,k) - rd~ 
     82      !! 
     83      !! ** Action  : compute for top interfaces 
     84      !!              - pgtu, pgtv: horizontal gradient of tracer at u- & v-points 
     85      !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points 
     86      !!---------------------------------------------------------------------- 
    8687      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index 
    8788      INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers 
     
    8990      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts  
    9091      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields 
    91       REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad. of prd at u- & v-pts  
    92       ! 
    93       INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
    94       INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points 
    95       REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars 
    96       REAL(wp), POINTER, DIMENSION(:,:  ) ::  zri, zrj, zhi, zhj 
    97       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zti, ztj    ! interpolated value of tracer 
    98       !!---------------------------------------------------------------------- 
    99       ! 
    100       IF( nn_timing == 1 )  CALL timing_start( 'zps_hde') 
    101       ! 
    102       CALL wrk_alloc( jpi, jpj,       zri, zrj, zhi, zhj )  
    103       CALL wrk_alloc( jpi, jpj, kjpt, zti, ztj           )  
     92      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad of prd at u- & v-pts (bottom) 
     93      ! 
     94      INTEGER  ::   ji, jj, jn                  ! Dummy loop indices 
     95      INTEGER  ::   iku, ikv, ikum1, ikvm1      ! partial step level (ocean bottom level) at u- and v-points 
     96      REAL(wp) ::   ze3wu, ze3wv, zmaxu, zmaxv  ! local scalars 
     97      REAL(wp), DIMENSION(jpi,jpj)      ::   zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
     98      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::   zti, ztj             !  
     99      !!---------------------------------------------------------------------- 
     100      ! 
     101      IF( nn_timing == 1 )   CALL timing_start( 'zps_hde') 
     102      ! 
     103      pgtu(:,:,:)=0._wp   ;   zti (:,:,:)=0._wp   ;   zhi (:,:  )=0._wp 
     104      pgtv(:,:,:)=0._wp   ;   ztj (:,:,:)=0._wp   ;   zhj (:,:  )=0._wp 
    104105      ! 
    105106      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
    106107         ! 
    107 # if defined key_vectopt_loop 
    108          jj = 1 
    109          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled) 
    110 # else 
    111          DO jj = 1, jpjm1 
    112             DO ji = 1, jpim1 
    113 # endif 
     108         DO jj = 1, jpjm1 
     109            DO ji = 1, jpim1 
    114110               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
    115111               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
    116                ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    117                ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
     112!!gm BUG ? when applied to before fields, e3w_b should be used.... 
     113               ze3wu = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku) 
     114               ze3wv = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv) 
    118115               ! 
    119116               ! i- direction 
    120117               IF( ze3wu >= 0._wp ) THEN      ! case 1 
    121                   zmaxu =  ze3wu / fse3w(ji+1,jj,iku) 
    122                   ! interpolated values of tracers 
    123                   zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     118                  zmaxu =  ze3wu / e3w_n(ji+1,jj,iku) 
     119                  ! interpolated values of tracers 
     120                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
    124121                  ! gradient of  tracers 
    125122                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
    126123               ELSE                           ! case 2 
    127                   zmaxu = -ze3wu / fse3w(ji,jj,iku) 
    128                   ! interpolated values of tracers 
    129                   zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
     124                  zmaxu = -ze3wu / e3w_n(ji,jj,iku) 
     125                  ! interpolated values of tracers 
     126                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
    130127                  ! gradient of tracers 
    131128                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     
    134131               ! j- direction 
    135132               IF( ze3wv >= 0._wp ) THEN      ! case 1 
    136                   zmaxv =  ze3wv / fse3w(ji,jj+1,ikv) 
    137                   ! interpolated values of tracers 
    138                   ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
     133                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikv) 
     134                  ! interpolated values of tracers 
     135                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
    139136                  ! gradient of tracers 
    140137                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
    141138               ELSE                           ! case 2 
    142                   zmaxv =  -ze3wv / fse3w(ji,jj,ikv) 
    143                   ! interpolated values of tracers 
    144                   ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
     139                  zmaxv =  -ze3wv / e3w_n(ji,jj,ikv) 
     140                  ! interpolated values of tracers 
     141                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
    145142                  ! gradient of tracers 
    146143                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
    147144               ENDIF 
    148 # if ! defined key_vectopt_loop 
    149             END DO 
    150 # endif 
     145            END DO 
     146         END DO 
     147         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     148         ! 
     149      END DO 
     150      !                 
     151      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part) 
     152         pgru(:,:) = 0._wp 
     153         pgrv(:,:) = 0._wp                ! depth of the partial step level 
     154         DO jj = 1, jpjm1 
     155            DO ji = 1, jpim1 
     156               iku = mbku(ji,jj) 
     157               ikv = mbkv(ji,jj) 
     158               ze3wu  = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku) 
     159               ze3wv  = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv) 
     160               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)     ! i-direction: case 1 
     161               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)     ! -     -      case 2 
     162               ENDIF 
     163               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)     ! j-direction: case 1 
     164               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)     ! -     -      case 2 
     165               ENDIF 
     166            END DO 
     167         END DO 
     168         ! 
     169         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj  
     170         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj  
     171         ! 
     172         DO jj = 1, jpjm1                 ! Gradient of density at the last level  
     173            DO ji = 1, jpim1 
     174               iku = mbku(ji,jj) 
     175               ikv = mbkv(ji,jj) 
     176               ze3wu  = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku) 
     177               ze3wv  = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv) 
     178               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1 
     179               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2 
     180               ENDIF 
     181               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1 
     182               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2 
     183               ENDIF 
     184            END DO 
     185         END DO 
     186         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
     187         ! 
     188      END IF 
     189      ! 
     190      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde') 
     191      ! 
     192   END SUBROUTINE zps_hde 
     193   ! 
     194   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  & 
     195      &                          prd, pgru, pgrv, pgrui, pgrvi ) 
     196      !!---------------------------------------------------------------------- 
     197      !!                     ***  ROUTINE zps_hde_isf  *** 
     198      !!                     
     199      !! ** Purpose :   Compute the horizontal derivative of T, S and rho 
     200      !!      at u- and v-points with a linear interpolation for z-coordinate 
     201      !!      with partial steps for top (ice shelf) and bottom. 
     202      !! 
     203      !! ** Method  :   In z-coord with partial steps, scale factors on last  
     204      !!      levels are different for each grid point, so that T, S and rd  
     205      !!      points are not at the same depth as in z-coord. To have horizontal 
     206      !!      gradients again, we interpolate T and S at the good depth : 
     207      !!      For the bottom case: 
     208      !!      Linear interpolation of T, S    
     209      !!         Computation of di(tb) and dj(tb) by vertical interpolation: 
     210      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~ 
     211      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~ 
     212      !!         This formulation computes the two cases: 
     213      !!                 CASE 1                   CASE 2   
     214      !!         k-1  ___ ___________   k-1   ___ ___________ 
     215      !!                    Ti  T~                  T~  Ti+1 
     216      !!                  _____                        _____ 
     217      !!         k        |   |Ti+1     k           Ti |   | 
     218      !!                  |   |____                ____|   | 
     219      !!              ___ |   |   |           ___  |   |   | 
     220      !!                   
     221      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then 
     222      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1) 
     223      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  ) 
     224      !!          or 
     225      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then 
     226      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i ) 
     227      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) ) 
     228      !!          Idem for di(s) and dj(s)           
     229      !! 
     230      !!      For rho, we call eos which will compute rd~(t~,s~) at the right 
     231      !!      depth zh from interpolated T and S for the different formulations 
     232      !!      of the equation of state (eos). 
     233      !!      Gradient formulation for rho : 
     234      !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~ 
     235      !! 
     236      !!      For the top case (ice shelf): As for the bottom case but upside down 
     237      !! 
     238      !! ** Action  : compute for top and bottom interfaces 
     239      !!              - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points 
     240      !!              - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points 
     241      !!---------------------------------------------------------------------- 
     242      INTEGER                              , INTENT(in   )           ::  kt           ! ocean time-step index 
     243      INTEGER                              , INTENT(in   )           ::  kjpt         ! number of tracers 
     244      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta          ! 4D tracers fields 
     245      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv   ! hor. grad. of ptra at u- & v-pts  
     246      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 
     247      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd          ! 3D density anomaly fields 
     248      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv   ! hor. grad of prd at u- & v-pts (bottom) 
     249      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) 
     250      ! 
     251      INTEGER  ::   ji, jj, jn      ! Dummy loop indices 
     252      INTEGER  ::   iku, ikv, ikum1, ikvm1,ikup1, ikvp1   ! partial step level (ocean bottom level) at u- and v-points 
     253      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv             ! temporary scalars 
     254      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos 
     255      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !  
     256      !!---------------------------------------------------------------------- 
     257      ! 
     258      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf') 
     259      ! 
     260      pgtu (:,:,:) = 0._wp   ;   pgtv (:,:,:) =0._wp 
     261      pgtui(:,:,:) = 0._wp   ;   pgtvi(:,:,:) =0._wp 
     262      zti  (:,:,:) = 0._wp   ;   ztj  (:,:,:) =0._wp 
     263      zhi  (:,:  ) = 0._wp   ;   zhj  (:,:  ) =0._wp 
     264      ! 
     265      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==! 
     266         ! 
     267         DO jj = 1, jpjm1 
     268            DO ji = 1, jpim1 
     269 
     270               iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points 
     271               ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1 
     272               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku) 
     273               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv) 
     274               ! 
     275               ! i- direction 
     276               IF( ze3wu >= 0._wp ) THEN      ! case 1 
     277                  zmaxu =  ze3wu / e3w_n(ji+1,jj,iku) 
     278                  ! interpolated values of tracers 
     279                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) ) 
     280                  ! gradient of  tracers 
     281                  pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     282               ELSE                           ! case 2 
     283                  zmaxu = -ze3wu / e3w_n(ji,jj,iku) 
     284                  ! interpolated values of tracers 
     285                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) ) 
     286                  ! gradient of tracers 
     287                  pgtu(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     288               ENDIF 
     289               ! 
     290               ! j- direction 
     291               IF( ze3wv >= 0._wp ) THEN      ! case 1 
     292                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikv) 
     293                  ! interpolated values of tracers 
     294                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) ) 
     295                  ! gradient of tracers 
     296                  pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     297               ELSE                           ! case 2 
     298                  zmaxv =  -ze3wv / e3w_n(ji,jj,ikv) 
     299                  ! interpolated values of tracers 
     300                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) ) 
     301                  ! gradient of tracers 
     302                  pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     303               ENDIF 
     304 
     305            END DO 
    151306         END DO 
    152307         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     
    156311      ! horizontal derivative of density anomalies (rd) 
    157312      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level 
    158 # if defined key_vectopt_loop 
    159          jj = 1 
    160          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled) 
    161 # else 
    162          DO jj = 1, jpjm1 
    163             DO ji = 1, jpim1 
    164 # endif 
     313         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ;  
     314         ! 
     315         DO jj = 1, jpjm1 
     316            DO ji = 1, jpim1 
     317 
    165318               iku = mbku(ji,jj) 
    166319               ikv = mbkv(ji,jj) 
    167                ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    168                ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
    169                IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1 
    170                ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2 
    171                ENDIF 
    172                IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1 
    173                ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2 
    174                ENDIF 
    175 # if ! defined key_vectopt_loop 
    176             END DO 
    177 # endif 
     320               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku) 
     321               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv) 
     322               ! 
     323               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)    ! i-direction: case 1 
     324               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)    ! -     -      case 2 
     325               ENDIF 
     326               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)    ! j-direction: case 1 
     327               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)    ! -     -      case 2 
     328               ENDIF 
     329 
     330            END DO 
    178331         END DO 
    179332 
    180333         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial 
    181334         ! step and store it in  zri, zrj for each  case 
    182          CALL eos( zti, zhi, zri )   
     335         CALL eos( zti, zhi, zri ) 
    183336         CALL eos( ztj, zhj, zrj ) 
    184337 
    185          ! Gradient of density at the last level  
    186 # if defined key_vectopt_loop 
    187          jj = 1 
    188          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled) 
    189 # else 
    190          DO jj = 1, jpjm1 
    191             DO ji = 1, jpim1 
    192 # endif 
     338         DO jj = 1, jpjm1                 ! Gradient of density at the last level  
     339            DO ji = 1, jpim1 
    193340               iku = mbku(ji,jj) 
    194341               ikv = mbkv(ji,jj) 
    195                ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku) 
    196                ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv) 
    197                IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1 
    198                ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2 
    199                ENDIF 
    200                IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1 
    201                ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2 
    202                ENDIF 
    203 # if ! defined key_vectopt_loop 
    204             END DO 
    205 # endif 
    206          END DO 
     342               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku) 
     343               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv) 
     344 
     345               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1 
     346               ELSE                        ;   pgru(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2 
     347               ENDIF 
     348               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1 
     349               ELSE                        ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2 
     350               ENDIF 
     351 
     352            END DO 
     353         END DO 
     354 
    207355         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions 
    208356         ! 
    209357      END IF 
    210358      ! 
    211       CALL wrk_dealloc( jpi, jpj,       zri, zrj, zhi, zhj )  
    212       CALL wrk_dealloc( jpi, jpj, kjpt, zti, ztj           )  
    213       ! 
    214       IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde') 
    215       ! 
    216    END SUBROUTINE zps_hde 
    217  
     359      !     !==  (ISH)  compute grui and gruvi  ==! 
     360      ! 
     361      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            ! 
     362         DO jj = 1, jpjm1 
     363            DO ji = 1, jpim1 
     364               iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1 
     365               ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1 
     366               ! 
     367               ! (ISF) case partial step top and bottom in adjacent cell in vertical 
     368               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom 
     369               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj 
     370               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0 
     371               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku) 
     372               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv)  
     373 
     374               ! i- direction 
     375               IF( ze3wu >= 0._wp ) THEN      ! case 1 
     376                  zmaxu = ze3wu / e3w_n(ji+1,jj,ikup1) 
     377                  ! interpolated values of tracers 
     378                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) ) 
     379                  ! gradient of tracers 
     380                  pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) ) 
     381               ELSE                           ! case 2 
     382                  zmaxu = - ze3wu / e3w_n(ji,jj,ikup1) 
     383                  ! interpolated values of tracers 
     384                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) ) 
     385                  ! gradient of  tracers 
     386                  pgtui(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) ) 
     387               ENDIF 
     388               ! 
     389               ! j- direction 
     390               IF( ze3wv >= 0._wp ) THEN      ! case 1 
     391                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikvp1) 
     392                  ! interpolated values of tracers 
     393                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) ) 
     394                  ! gradient of tracers 
     395                  pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) ) 
     396               ELSE                           ! case 2 
     397                  zmaxv =  - ze3wv / e3w_n(ji,jj,ikvp1) 
     398                  ! interpolated values of tracers 
     399                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) ) 
     400                  ! gradient of tracers 
     401                  pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) ) 
     402               ENDIF 
     403 
     404            END DO 
     405         END DO 
     406         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ); CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond. 
     407         ! 
     408      END DO 
     409 
     410      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part) 
     411         ! 
     412         pgrui(:,:)  =0.0_wp; pgrvi(:,:)  =0.0_wp; 
     413         DO jj = 1, jpjm1 
     414            DO ji = 1, jpim1 
     415 
     416               iku = miku(ji,jj) 
     417               ikv = mikv(ji,jj) 
     418               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku) 
     419               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv)  
     420               ! 
     421               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)    ! i-direction: case 1 
     422               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)    ! -     -      case 2 
     423               ENDIF 
     424 
     425               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)    ! j-direction: case 1 
     426               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)    ! -     -      case 2 
     427               ENDIF 
     428 
     429            END DO 
     430         END DO 
     431         ! 
     432         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj  
     433         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj  
     434         ! 
     435         DO jj = 1, jpjm1                 ! Gradient of density at the last level  
     436            DO ji = 1, jpim1 
     437               iku = miku(ji,jj)  
     438               ikv = mikv(ji,jj)  
     439               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku) 
     440               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv)  
     441 
     442               IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) ) ! i: 1 
     443               ELSE                      ; pgrui(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj  ,iku) - zri(ji,jj    ) ) ! i: 2 
     444               ENDIF 
     445               IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( zrj(ji  ,jj      ) - prd(ji,jj,ikv) ) ! j: 1 
     446               ELSE                      ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( prd(ji  ,jj+1,ikv) - zrj(ji,jj    ) ) ! j: 2 
     447               ENDIF 
     448 
     449            END DO 
     450         END DO 
     451         CALL lbc_lnk( pgrui   , 'U', -1. ); CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions 
     452         ! 
     453      END IF   
     454      ! 
     455      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde_isf') 
     456      ! 
     457   END SUBROUTINE zps_hde_isf 
    218458   !!====================================================================== 
    219459END MODULE zpshde 
Note: See TracChangeset for help on using the changeset viewer.