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 trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

update trunk with OpenMP parallelization
  • Property svn:keywords set to Id
File size: 25.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   !!            3.6  !  2014-11  (P. Mathiot) Add zps_hde_isf (needed to open a cavity)
11   !!======================================================================
12   
13   !!----------------------------------------------------------------------
14   !!   zps_hde      :  Horizontal DErivative of T, S and rd at the last
15   !!                   ocean level (Z-coord. with Partial Steps)
16   !!----------------------------------------------------------------------
17   USE oce             ! ocean: dynamics and tracers variables
18   USE dom_oce         ! domain: ocean variables
19   USE phycst          ! physical constants
20   USE eosbn2          ! ocean equation of state
21   USE in_out_manager  ! I/O manager
22   USE lbclnk          ! lateral boundary conditions (or mpp link)
23   USE lib_mpp         ! MPP library
24   USE wrk_nemo        ! Memory allocation
25   USE timing          ! Timing
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC   zps_hde     ! routine called by step.F90
31   PUBLIC   zps_hde_isf ! routine called by step.F90
32
33   !! * Substitutions
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
37   !! $Id$
38   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
39   !!----------------------------------------------------------------------
40CONTAINS
41
42   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   &
43      &                          prd, pgru, pgrv    )
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 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).
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 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      !!----------------------------------------------------------------------
87      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index
88      INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers
89      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields
90      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts
91      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields
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      DO jn = 1, kjpt
104!$OMP PARALLEL DO schedule(static) private(jj,ji)
105         DO jj = 1, jpjm1
106            DO ji = 1, jpim1
107               pgtu(ji,jj,jn)=0._wp   ;   zti (ji,jj,jn)=0._wp
108               pgtv(ji,jj,jn)=0._wp   ;   ztj (ji,jj,jn)=0._wp
109            END DO
110         END DO
111      END DO
112!$OMP PARALLEL DO schedule(static) private(jj,ji)
113      DO jj = 1, jpjm1
114         DO ji = 1, jpim1
115            zhi (ji,jj  )=0._wp
116            zhj (ji,jj  )=0._wp
117         END DO
118       END DO
119      !
120      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
121         !
122!$OMP PARALLEL DO schedule(static) private(jj,ji,iku,ikv,ze3wu,ze3wv,zmaxu,zmaxv)
123         DO jj = 1, jpjm1
124            DO ji = 1, jpim1
125               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
126               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
127!!gm BUG ? when applied to before fields, e3w_b should be used....
128               ze3wu = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku)
129               ze3wv = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv)
130               !
131               ! i- direction
132               IF( ze3wu >= 0._wp ) THEN      ! case 1
133                  zmaxu =  ze3wu / e3w_n(ji+1,jj,iku)
134                  ! interpolated values of tracers
135                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
136                  ! gradient of  tracers
137                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
138               ELSE                           ! case 2
139                  zmaxu = -ze3wu / e3w_n(ji,jj,iku)
140                  ! interpolated values of tracers
141                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
142                  ! gradient of tracers
143                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
144               ENDIF
145               !
146               ! j- direction
147               IF( ze3wv >= 0._wp ) THEN      ! case 1
148                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikv)
149                  ! interpolated values of tracers
150                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
151                  ! gradient of tracers
152                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
153               ELSE                           ! case 2
154                  zmaxv =  -ze3wv / e3w_n(ji,jj,ikv)
155                  ! interpolated values of tracers
156                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
157                  ! gradient of tracers
158                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
159               ENDIF
160            END DO
161         END DO
162         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
163         !
164      END DO
165      !               
166      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
167!$OMP PARALLEL
168!$OMP DO schedule(static) private(jj,ji)
169         DO jj = 1, jpjm1
170            DO ji = 1, jpim1
171               pgru(ji,jj) = 0._wp
172               pgrv(ji,jj) = 0._wp                ! depth of the partial step level
173            END DO
174         END DO
175!$OMP END DO NOWAIT
176!$OMP DO schedule(static) private(jj,ji,iku,ikv,ze3wu,ze3wv)
177         DO jj = 1, jpjm1
178            DO ji = 1, jpim1
179               iku = mbku(ji,jj)
180               ikv = mbkv(ji,jj)
181               ze3wu  = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku)
182               ze3wv  = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv)
183               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)     ! i-direction: case 1
184               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)     ! -     -      case 2
185               ENDIF
186               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)     ! j-direction: case 1
187               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)     ! -     -      case 2
188               ENDIF
189            END DO
190         END DO
191!$OMP END DO NOWAIT
192!$OMP END PARALLEL
193         !
194         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj
195         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj
196         !
197!$OMP PARALLEL DO schedule(static) private(jj,ji,iku,ikv,ze3wu,ze3wv)
198         DO jj = 1, jpjm1                 ! Gradient of density at the last level
199            DO ji = 1, jpim1
200               iku = mbku(ji,jj)
201               ikv = mbkv(ji,jj)
202               ze3wu  = e3w_n(ji+1,jj  ,iku) - e3w_n(ji,jj,iku)
203               ze3wv  = e3w_n(ji  ,jj+1,ikv) - e3w_n(ji,jj,ikv)
204               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
205               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2
206               ENDIF
207               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1
208               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2
209               ENDIF
210            END DO
211         END DO
212         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions
213         !
214      END IF
215      !
216      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde')
217      !
218   END SUBROUTINE zps_hde
219   !
220   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  &
221      &                          prd, pgru, pgrv, pgrui, pgrvi )
222      !!----------------------------------------------------------------------
223      !!                     ***  ROUTINE zps_hde_isf  ***
224      !!                   
225      !! ** Purpose :   Compute the horizontal derivative of T, S and rho
226      !!      at u- and v-points with a linear interpolation for z-coordinate
227      !!      with partial steps for top (ice shelf) and bottom.
228      !!
229      !! ** Method  :   In z-coord with partial steps, scale factors on last
230      !!      levels are different for each grid point, so that T, S and rd
231      !!      points are not at the same depth as in z-coord. To have horizontal
232      !!      gradients again, we interpolate T and S at the good depth :
233      !!      For the bottom case:
234      !!      Linear interpolation of T, S   
235      !!         Computation of di(tb) and dj(tb) by vertical interpolation:
236      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~
237      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~
238      !!         This formulation computes the two cases:
239      !!                 CASE 1                   CASE 2 
240      !!         k-1  ___ ___________   k-1   ___ ___________
241      !!                    Ti  T~                  T~  Ti+1
242      !!                  _____                        _____
243      !!         k        |   |Ti+1     k           Ti |   |
244      !!                  |   |____                ____|   |
245      !!              ___ |   |   |           ___  |   |   |
246      !!                 
247      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
248      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
249      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  )
250      !!          or
251      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
252      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
253      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
254      !!          Idem for di(s) and dj(s)         
255      !!
256      !!      For rho, we call eos which will compute rd~(t~,s~) at the right
257      !!      depth zh from interpolated T and S for the different formulations
258      !!      of the equation of state (eos).
259      !!      Gradient formulation for rho :
260      !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~
261      !!
262      !!      For the top case (ice shelf): As for the bottom case but upside down
263      !!
264      !! ** Action  : compute for top and bottom interfaces
265      !!              - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points
266      !!              - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points
267      !!----------------------------------------------------------------------
268      INTEGER                              , INTENT(in   )           ::  kt           ! ocean time-step index
269      INTEGER                              , INTENT(in   )           ::  kjpt         ! number of tracers
270      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta          ! 4D tracers fields
271      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv   ! hor. grad. of ptra at u- & v-pts
272      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF)
273      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd          ! 3D density anomaly fields
274      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv   ! hor. grad of prd at u- & v-pts (bottom)
275      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top)
276      !
277      INTEGER  ::   ji, jj, jn      ! Dummy loop indices
278      INTEGER  ::   iku, ikv, ikum1, ikvm1,ikup1, ikvp1   ! partial step level (ocean bottom level) at u- and v-points
279      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv             ! temporary scalars
280      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos
281      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !
282      !!----------------------------------------------------------------------
283      !
284      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf')
285      !
286      pgtu (:,:,:) = 0._wp   ;   pgtv (:,:,:) =0._wp
287      pgtui(:,:,:) = 0._wp   ;   pgtvi(:,:,:) =0._wp
288      zti  (:,:,:) = 0._wp   ;   ztj  (:,:,:) =0._wp
289      zhi  (:,:  ) = 0._wp   ;   zhj  (:,:  ) =0._wp
290      !
291      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
292         !
293         DO jj = 1, jpjm1
294            DO ji = 1, jpim1
295
296               iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
297               ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
298               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku)
299               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv)
300               !
301               ! i- direction
302               IF( ze3wu >= 0._wp ) THEN      ! case 1
303                  zmaxu =  ze3wu / e3w_n(ji+1,jj,iku)
304                  ! interpolated values of tracers
305                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
306                  ! gradient of  tracers
307                  pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
308               ELSE                           ! case 2
309                  zmaxu = -ze3wu / e3w_n(ji,jj,iku)
310                  ! interpolated values of tracers
311                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
312                  ! gradient of tracers
313                  pgtu(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
314               ENDIF
315               !
316               ! j- direction
317               IF( ze3wv >= 0._wp ) THEN      ! case 1
318                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikv)
319                  ! interpolated values of tracers
320                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
321                  ! gradient of tracers
322                  pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
323               ELSE                           ! case 2
324                  zmaxv =  -ze3wv / e3w_n(ji,jj,ikv)
325                  ! interpolated values of tracers
326                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
327                  ! gradient of tracers
328                  pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
329               ENDIF
330
331            END DO
332         END DO
333         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
334         !
335      END DO
336
337      ! horizontal derivative of density anomalies (rd)
338      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
339         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ; 
340         !
341         DO jj = 1, jpjm1
342            DO ji = 1, jpim1
343
344               iku = mbku(ji,jj)
345               ikv = mbkv(ji,jj)
346               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku)
347               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv)
348               !
349               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)    ! i-direction: case 1
350               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)    ! -     -      case 2
351               ENDIF
352               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)    ! j-direction: case 1
353               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)    ! -     -      case 2
354               ENDIF
355
356            END DO
357         END DO
358
359         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
360         ! step and store it in  zri, zrj for each  case
361         CALL eos( zti, zhi, zri )
362         CALL eos( ztj, zhj, zrj )
363
364         DO jj = 1, jpjm1                 ! Gradient of density at the last level
365            DO ji = 1, jpim1
366               iku = mbku(ji,jj)
367               ikv = mbkv(ji,jj)
368               ze3wu = gdept_n(ji+1,jj,iku) - gdept_n(ji,jj,iku)
369               ze3wv = gdept_n(ji,jj+1,ikv) - gdept_n(ji,jj,ikv)
370
371               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
372               ELSE                        ;   pgru(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2
373               ENDIF
374               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1
375               ELSE                        ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2
376               ENDIF
377
378            END DO
379         END DO
380
381         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions
382         !
383      END IF
384      !
385      !     !==  (ISH)  compute grui and gruvi  ==!
386      !
387      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            !
388         DO jj = 1, jpjm1
389            DO ji = 1, jpim1
390               iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1
391               ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1
392               !
393               ! (ISF) case partial step top and bottom in adjacent cell in vertical
394               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
395               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
396               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
397               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku)
398               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 
399
400               ! i- direction
401               IF( ze3wu >= 0._wp ) THEN      ! case 1
402                  zmaxu = ze3wu / e3w_n(ji+1,jj,ikup1)
403                  ! interpolated values of tracers
404                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) )
405                  ! gradient of tracers
406                  pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
407               ELSE                           ! case 2
408                  zmaxu = - ze3wu / e3w_n(ji,jj,ikup1)
409                  ! interpolated values of tracers
410                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) )
411                  ! gradient of  tracers
412                  pgtui(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
413               ENDIF
414               !
415               ! j- direction
416               IF( ze3wv >= 0._wp ) THEN      ! case 1
417                  zmaxv =  ze3wv / e3w_n(ji,jj+1,ikvp1)
418                  ! interpolated values of tracers
419                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) )
420                  ! gradient of tracers
421                  pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
422               ELSE                           ! case 2
423                  zmaxv =  - ze3wv / e3w_n(ji,jj,ikvp1)
424                  ! interpolated values of tracers
425                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) )
426                  ! gradient of tracers
427                  pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
428               ENDIF
429
430            END DO
431         END DO
432         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. ); CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
433         !
434      END DO
435
436      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
437         !
438         pgrui(:,:)  =0.0_wp; pgrvi(:,:)  =0.0_wp;
439         DO jj = 1, jpjm1
440            DO ji = 1, jpim1
441
442               iku = miku(ji,jj)
443               ikv = mikv(ji,jj)
444               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku)
445               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 
446               !
447               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept_n(ji  ,jj,iku)    ! i-direction: case 1
448               ELSE                        ;   zhi(ji,jj) = gdept_n(ji+1,jj,iku)    ! -     -      case 2
449               ENDIF
450
451               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept_n(ji,jj  ,ikv)    ! j-direction: case 1
452               ELSE                        ;   zhj(ji,jj) = gdept_n(ji,jj+1,ikv)    ! -     -      case 2
453               ENDIF
454
455            END DO
456         END DO
457         !
458         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj
459         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj
460         !
461         DO jj = 1, jpjm1                 ! Gradient of density at the last level
462            DO ji = 1, jpim1
463               iku = miku(ji,jj) 
464               ikv = mikv(ji,jj) 
465               ze3wu  =  gdept_n(ji,jj,iku) - gdept_n(ji+1,jj,iku)
466               ze3wv  =  gdept_n(ji,jj,ikv) - gdept_n(ji,jj+1,ikv) 
467
468               IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) ) ! i: 1
469               ELSE                      ; pgrui(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj  ,iku) - zri(ji,jj    ) ) ! i: 2
470               ENDIF
471               IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( zrj(ji  ,jj      ) - prd(ji,jj,ikv) ) ! j: 1
472               ELSE                      ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( prd(ji  ,jj+1,ikv) - zrj(ji,jj    ) ) ! j: 2
473               ENDIF
474
475            END DO
476         END DO
477         CALL lbc_lnk( pgrui   , 'U', -1. ); CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions
478         !
479      END IF 
480      !
481      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde_isf')
482      !
483   END SUBROUTINE zps_hde_isf
484   !!======================================================================
485END MODULE zpshde
Note: See TracBrowser for help on using the repository browser.