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.
zpshde.F90 in branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90 @ 4747

Last change on this file since 4747 was 4747, checked in by mathiot, 10 years ago

ISF branch: change to deal with non mask bathymetry (land processor definition, building bathy and ice shelf draft variable), update of hpg (definition of ze3wu in case of zps and vvl) and bfr (in case of 2 cell water column thickness, each cell feels top and bottom friction).

  • Property svn:keywords set to Id
File size: 24.4 KB
Line 
1MODULE zpshde
2   !!======================================================================
3   !!                       ***  MODULE zpshde   ***
4   !! z-coordinate + partial step : Horizontal Derivative at ocean bottom level
5   !!======================================================================
6   !! History :  OPA  !  2002-04  (A. Bozec)  Original code
7   !!   NEMO     1.0  !  2002-08  (G. Madec E. Durand)  Optimization and Free form
8   !!             -   !  2004-03  (C. Ethe)  adapted for passive tracers
9   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
10   !!======================================================================
11   
12   !!----------------------------------------------------------------------
13   !!   zps_hde      :  Horizontal DErivative of T, S and rd at the last
14   !!                   ocean level (Z-coord. with Partial Steps)
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean: dynamics and tracers variables
17   USE dom_oce         ! domain: ocean variables
18   USE phycst          ! physical constants
19   USE eosbn2          ! ocean equation of state
20   USE in_out_manager  ! I/O manager
21   USE lbclnk          ! lateral boundary conditions (or mpp link)
22   USE lib_mpp         ! MPP library
23   USE wrk_nemo        ! Memory allocation
24   USE timing          ! Timing
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   zps_hde    ! routine called by step.F90
30
31   !! * Substitutions
32#  include "domzgr_substitute.h90"
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
36   !! $Id$
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   &
42      &                          prd, pgru, pgrv, pmru, pmrv, pgzu, pgzv, pge3ru, pge3rv,  &
43      &                   sgtu, sgtv, sgru, sgrv, smru, smrv, sgzu, sgzv, sge3ru, sge3rv )
44      !!----------------------------------------------------------------------
45      !!                     ***  ROUTINE zps_hde  ***
46      !!                   
47      !! ** Purpose :   Compute the horizontal derivative of T, S and rho
48      !!      at u- and v-points with a linear interpolation for z-coordinate
49      !!      with partial steps.
50      !!
51      !! ** Method  :   In z-coord with partial steps, scale factors on last
52      !!      levels are different for each grid point, so that T, S and rd
53      !!      points are not at the same depth as in z-coord. To have horizontal
54      !!      gradients again, we interpolate T and S at the good depth :
55      !!      Linear interpolation of T, S   
56      !!         Computation of di(tb) and dj(tb) by vertical interpolation:
57      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~
58      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~
59      !!         This formulation computes the two cases:
60      !!                 CASE 1                   CASE 2 
61      !!         k-1  ___ ___________   k-1   ___ ___________
62      !!                    Ti  T~                  T~  Ti+1
63      !!                  _____                        _____
64      !!         k        |   |Ti+1     k           Ti |   |
65      !!                  |   |____                ____|   |
66      !!              ___ |   |   |           ___  |   |   |
67      !!                 
68      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
69      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
70      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  )
71      !!          or
72      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
73      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
74      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
75      !!          Idem for di(s) and dj(s)         
76      !!
77      !!      For rho, we call eos_insitu_2d which will compute rd~(t~,s~) at
78      !!      the good depth zh from interpolated T and S for the different
79      !!      formulation of the equation of state (eos).
80      !!      Gradient formulation for rho :
81      !!          di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~
82      !!
83      !! ** Action  : compute for top and bottom interfaces
84      !!              - pgtu, pgtv, sgtu, sgtv: horizontal gradient of tracer at u- & v-points
85      !!              - pgru, pgrv, sgru, sgtv: horizontal gradient of rho (if present) at u- & v-points
86      !!              - pmru, pmrv, smru, smrv: horizontal sum of rho at u- & v- point (used in dynhpg with vvl)
87      !!              - pgzu, pgzv, sgzu, sgzv: horizontal gradient of z at u- and v- point (used in dynhpg with vvl)
88      !!              - pge3ru, pge3rv, sge3ru, sge3rv: horizontal gradient of rho weighted by local e3w at u- & v-points
89      !!----------------------------------------------------------------------
90      !
91      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index
92      INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers
93      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields
94      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts
95      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  sgtu, sgtv  ! hor. grad. of stra at u- & v-pts (ISF)
96      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields
97      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv      ! hor. grad of prd at u- & v-pts (bottom)
98      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru, pmrv      ! hor. sum  of prd at u- & v-pts (bottom)
99      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu, pgzv      ! hor. grad of z   at u- & v-pts (bottom)
100      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru, pge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (bottom)
101      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sgru, sgrv      ! hor. grad of prd at u- & v-pts (top)
102      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  smru, smrv      ! hor. sum  of prd at u- & v-pts (top)
103      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sgzu, sgzv      ! hor. grad of z   at u- & v-pts (top)
104      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  sge3ru, sge3rv  ! hor. grad of prd weighted by local e3w at u- & v-pts (top)
105      !
106      INTEGER  ::   ji, jj, jn      ! Dummy loop indices
107      INTEGER  ::   iku, ikv, ikum1, ikvm1,ikup1, ikvp1   ! partial step level (ocean bottom level) at u- and v-points
108      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv, zdzwu, zdzwv, zdzwuip1, zdzwvjp1  ! temporary scalars
109      REAL(wp), POINTER, DIMENSION(:,:  ) ::  zri, zrj, zhi, zhj
110      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zti, ztj    ! interpolated value of tracer
111      !!----------------------------------------------------------------------
112      !
113      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde')
114      !
115      CALL wrk_alloc( jpi, jpj,       zri, zrj, zhi, zhj ) 
116      CALL wrk_alloc( jpi, jpj, kjpt, zti, ztj           ) 
117      !
118      pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ; pgtu(:,:,:)=0.0_wp ; pgtv(:,:,:)=0.0_wp ;
119      !
120      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
121         !
122# if defined key_vectopt_loop
123         jj = 1
124         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
125# else
126         DO jj = 1, jpjm1
127            DO ji = 1, jpim1
128# endif
129               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
130               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
131               ! (ISF) case partial step top and bottom in adjacent cell in vertical
132               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
133               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
134               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
135               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
136               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
137               !
138               ! i- direction
139               IF (iku .GT. 1) THEN
140               IF( ze3wu >= 0._wp ) THEN      ! case 1
141                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
142                  ! interpolated values of tracers
143                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
144                  ! gradient of  tracers
145                  pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
146                  pgzu(ji,jj)    = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku)
147               ELSE                           ! case 2
148                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
149                  ! interpolated values of tracers
150                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
151                  ! gradient of tracers
152                  pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
153                  pgzu(ji,jj)    = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu)
154               ENDIF
155               ENDIF
156               !
157               ! j- direction
158               IF (ikv .GT. 1) THEN
159               IF( ze3wv >= 0._wp ) THEN      ! case 1
160                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
161                  ! interpolated values of tracers
162                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
163                  ! gradient of tracers
164                  pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
165                  pgzv(ji,jj)    = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv) 
166               ELSE                           ! case 2
167                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
168                  ! interpolated values of tracers
169                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
170                  ! gradient of tracers
171                  pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
172                  pgzv(ji,jj)    = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv)
173               ENDIF
174               ENDIF
175# if ! defined key_vectopt_loop
176            END DO
177# endif
178         END DO
179         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
180         !
181      END DO
182
183      ! horizontal derivative of density anomalies (rd)
184      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
185# if defined key_vectopt_loop
186         jj = 1
187         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
188# else
189         DO jj = 1, jpjm1
190            DO ji = 1, jpim1
191# endif
192               iku = mbku(ji,jj)
193               ikv = mbkv(ji,jj)
194               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
195               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
196
197               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1
198               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2
199               ENDIF
200               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv    ! j-direction: case 1
201               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) + ze3wv    ! -     -      case 2
202               ENDIF
203# if ! defined key_vectopt_loop
204            END DO
205# endif
206         END DO
207         
208         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
209         ! step and store it in  zri, zrj for each  case
210         CALL eos( zti, zhi, zri ) 
211         CALL eos( ztj, zhj, zrj )
212
213         ! Gradient of density at the last level
214# if defined key_vectopt_loop
215         jj = 1
216         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
217# else
218         DO jj = 1, jpjm1
219            DO ji = 1, jpim1
220# endif
221               iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
222               ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 )    ! last and before last ocean level at u- & v-points
223               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
224               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
225               IF( ze3wu >= 0._wp ) THEN
226                  pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1
227                  pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i: 1
228                  pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  &
229                                * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji  ,jj    ) + prd(ji+1,jj,ikum1) + 2._wp) &
230                                   - fse3w(ji  ,jj,iku)          * ( prd(ji  ,jj,iku) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2
231               ELSE 
232                  pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2
233                  pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) )   ! i: 2
234                  pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  &
235                                * (  fse3w(ji+1,jj,iku)          * ( prd(ji+1,jj,iku) + prd(ji+1,jj,ikum1) + 2._wp) &
236                                   -(fse3w(ji  ,jj,iku) + ze3wu) * ( zri(ji  ,jj    ) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2
237               ENDIF
238               IF( ze3wv >= 0._wp ) THEN
239                  pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1
240                  pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )   ! j: 1
241                  pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  &
242                                * ( (fse3w(ji,jj+1,ikv) - ze3wv )* ( zrj(ji,jj      ) + prd(ji,jj+1,ikvm1) + 2._wp) &
243                                   - fse3w(ji,jj  ,ikv)          * ( prd(ji,jj  ,ikv) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2
244               ELSE
245                  pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2
246                  pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )   ! j: 2
247                  pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  &
248                                * (  fse3w(ji,jj+1,ikv)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikvm1) + 2._wp) &
249                                   -(fse3w(ji,jj  ,ikv) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2
250               ENDIF
251# if ! defined key_vectopt_loop
252            END DO
253# endif
254         END DO
255         CALL lbc_lnk( pgru   , 'U', -1. )   ;   CALL lbc_lnk( pgrv   , 'V', -1. )   ! Lateral boundary conditions
256         CALL lbc_lnk( pmru   , 'U',  1. )   ;   CALL lbc_lnk( pmrv   , 'V',  1. )   ! Lateral boundary conditions
257         CALL lbc_lnk( pgzu   , 'U', -1. )   ;   CALL lbc_lnk( pgzv   , 'V', -1. )   ! Lateral boundary conditions
258         CALL lbc_lnk( pge3ru , 'U', -1. )   ;   CALL lbc_lnk( pge3rv , 'V', -1. )   ! Lateral boundary conditions
259         !
260      END IF
261         ! (ISH)  compute grui and gruvi
262      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            !
263# if defined key_vectopt_loop
264         jj = 1
265         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
266# else
267         DO jj = 1, jpjm1
268            DO ji = 1, jpim1
269# endif
270               iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1
271               ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1
272               !
273               ! (ISF) case partial step top and bottom in adjacent cell in vertical
274               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
275               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
276               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
277               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 
278               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
279               ! i- direction
280               IF( ze3wu >= 0._wp ) THEN      ! case 1
281                  zmaxu = ze3wu / fse3w(ji+1,jj,iku+1)
282                  ! interpolated values of tracers
283                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) )
284                  ! gradient of tracers
285                  sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
286                  sgzu(ji,jj)    = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku)
287               ELSE                           ! case 2
288                  zmaxu = - ze3wu / fse3w(ji,jj,iku+1)
289                  ! interpolated values of tracers
290                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) )
291                  sgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
292                  ! gradient of  tracers
293                  sgzu(ji,jj)    = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu)
294               ENDIF
295               !
296               ! j- direction
297               IF( ze3wv >= 0._wp ) THEN      ! case 1
298                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv+1)
299                  ! interpolated values of tracers
300                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) )
301                  ! gradient of tracers
302                  sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
303                  sgzv(ji,jj)    = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv) 
304               ELSE                           ! case 2
305                  zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1)
306                  ! interpolated values of tracers
307                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) )
308                  ! gradient of tracers
309                  sgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
310                  sgzv(ji,jj)    = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv)
311               ENDIF
312# if ! defined key_vectopt_loop
313            END DO
314# endif
315         END DO!!
316         CALL lbc_lnk( sgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( sgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
317         !
318      END DO
319
320      ! horizontal derivative of density anomalies (rd)
321      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
322# if defined key_vectopt_loop
323         jj = 1
324         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
325# else
326         DO jj = 1, jpjm1
327            DO ji = 1, jpim1
328# endif
329               iku = miku(ji,jj)
330               ikv = mikv(ji,jj)
331               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
332               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
333
334               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1
335               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2
336               ENDIF
337               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1
338               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2
339               ENDIF
340# if ! defined key_vectopt_loop
341            END DO
342# endif
343         END DO
344
345         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
346         ! step and store it in  zri, zrj for each  case
347         CALL eos( zti, zhi, zri ) 
348         CALL eos( ztj, zhj, zrj )
349
350         ! Gradient of density at the last level
351# if defined key_vectopt_loop
352         jj = 1
353         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
354# else
355         DO jj = 1, jpjm1
356            DO ji = 1, jpim1
357# endif
358               iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1
359               ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1
360               ze3wu  = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
361               ze3wv  = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
362               IF( ze3wu >= 0._wp ) THEN
363                 sgru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1
364                 smru  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1
365                 sge3ru(ji,jj) = umask(ji,jj,iku+1)                                                                  &
366                                * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   &
367                                   - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1
368               ELSE
369                 sgru  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2
370                 smru  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2
371                 sge3ru(ji,jj) = umask(ji,jj,iku+1)                                                                   &
372                                * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  &
373                                   -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2
374               ENDIF
375               IF( ze3wv >= 0._wp ) THEN
376                 sgrv  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1
377                 smrv  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1
378                 sge3rv(ji,jj) = vmask(ji,jj,ikv+1)                                                                  & 
379                                * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  &
380                                   - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1
381                                  ! + 2 due to the formulation in density and not in anomalie in hpg sco
382               ELSE
383                 sgrv  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2
384                 smrv  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2
385                 sge3rv(ji,jj) = vmask(ji,jj,ikv+1)                                                                   &
386                                * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) &
387                                   -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2
388               ENDIF
389# if ! defined key_vectopt_loop
390            END DO
391# endif
392         END DO
393         CALL lbc_lnk( sgru   , 'U', -1. )   ;   CALL lbc_lnk( sgrv   , 'V', -1. )   ! Lateral boundary conditions
394         CALL lbc_lnk( smru   , 'U',  1. )   ;   CALL lbc_lnk( smrv   , 'V',  1. )   ! Lateral boundary conditions
395         CALL lbc_lnk( sgzu   , 'U', -1. )   ;   CALL lbc_lnk( sgzv   , 'V', -1. )   ! Lateral boundary conditions
396         CALL lbc_lnk( sge3ru , 'U', -1. )   ;   CALL lbc_lnk( sge3rv , 'V', -1. )   ! Lateral boundary conditions
397         !
398      END IF 
399      !
400      CALL wrk_dealloc( jpi, jpj,       zri, zrj, zhi, zhj) 
401      CALL wrk_dealloc( jpi, jpj, kjpt, zti, ztj           ) 
402      !
403      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde')
404      !
405   END SUBROUTINE zps_hde
406
407   !!======================================================================
408END MODULE zpshde
Note: See TracBrowser for help on using the repository browser.