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 NEMO/trunk/src/OCE/TRA – NEMO

source: NEMO/trunk/src/OCE/TRA/zpshde.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
​svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.
  • Property svn:keywords set to Id
File size: 23.8 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 timing          ! Timing
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   zps_hde     ! routine called by step.F90
30   PUBLIC   zps_hde_isf ! routine called by step.F90
31
32   !! * Substitutions
33#  include "do_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
36   !! $Id$
37   !! Software governed by the CeCILL license (see ./LICENSE)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE zps_hde( kt, Kmm, kjpt, pta, pgtu, pgtv,   &
42      &                          prd, pgru, pgrv    )
43      !!----------------------------------------------------------------------
44      !!                     ***  ROUTINE zps_hde  ***
45      !!                   
46      !! ** Purpose :   Compute the horizontal derivative of T, S and rho
47      !!      at u- and v-points with a linear interpolation for z-coordinate
48      !!      with partial steps.
49      !!
50      !! ** Method  :   In z-coord with partial steps, scale factors on last
51      !!      levels are different for each grid point, so that T, S and rd
52      !!      points are not at the same depth as in z-coord. To have horizontal
53      !!      gradients again, we interpolate T and S at the good depth :
54      !!      Linear interpolation of T, S   
55      !!         Computation of di(tb) and dj(tb) by vertical interpolation:
56      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~
57      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~
58      !!         This formulation computes the two cases:
59      !!                 CASE 1                   CASE 2 
60      !!         k-1  ___ ___________   k-1   ___ ___________
61      !!                    Ti  T~                  T~  Ti+1
62      !!                  _____                        _____
63      !!         k        |   |Ti+1     k           Ti |   |
64      !!                  |   |____                ____|   |
65      !!              ___ |   |   |           ___  |   |   |
66      !!                 
67      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
68      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
69      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  )
70      !!          or
71      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
72      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
73      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
74      !!          Idem for di(s) and dj(s)         
75      !!
76      !!      For rho, we call eos which will compute rd~(t~,s~) at the right
77      !!      depth zh from interpolated T and S for the different formulations
78      !!      of the equation of state (eos).
79      !!      Gradient formulation for rho :
80      !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~
81      !!
82      !! ** Action  : compute for top interfaces
83      !!              - pgtu, pgtv: horizontal gradient of tracer at u- & v-points
84      !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points
85      !!----------------------------------------------------------------------
86      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index
87      INTEGER                              , INTENT(in   )           ::  Kmm         ! ocean time level 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( ln_timing )   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
105      !
106      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
107         !
108         DO_2D_10_10
109            iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
110            ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
111!!gm BUG ? when applied to before fields, e3w(:,:,:,Kbb) should be used....
112            ze3wu = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
113            ze3wv = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
114            !
115            ! i- direction
116            IF( ze3wu >= 0._wp ) THEN      ! case 1
117               zmaxu =  ze3wu / e3w(ji+1,jj,iku,Kmm)
118               ! interpolated values of tracers
119               zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
120               ! gradient of  tracers
121               pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
122            ELSE                           ! case 2
123               zmaxu = -ze3wu / e3w(ji,jj,iku,Kmm)
124               ! interpolated values of tracers
125               zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
126               ! gradient of tracers
127               pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
128            ENDIF
129            !
130            ! j- direction
131            IF( ze3wv >= 0._wp ) THEN      ! case 1
132               zmaxv =  ze3wv / e3w(ji,jj+1,ikv,Kmm)
133               ! interpolated values of tracers
134               ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
135               ! gradient of tracers
136               pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
137            ELSE                           ! case 2
138               zmaxv =  -ze3wv / e3w(ji,jj,ikv,Kmm)
139               ! interpolated values of tracers
140               ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
141               ! gradient of tracers
142               pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
143            ENDIF
144         END_2D
145      END DO
146      !
147      CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1. , pgtv(:,:,:), 'V', -1. )   ! Lateral boundary cond.
148      !               
149      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
150         pgru(:,:) = 0._wp
151         pgrv(:,:) = 0._wp                ! depth of the partial step level
152         DO_2D_10_10
153            iku = mbku(ji,jj)
154            ikv = mbkv(ji,jj)
155            ze3wu  = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
156            ze3wv  = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
157            IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku,Kmm)     ! i-direction: case 1
158            ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm)     ! -     -      case 2
159            ENDIF
160            IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv,Kmm)     ! j-direction: case 1
161            ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm)     ! -     -      case 2
162            ENDIF
163         END_2D
164         !
165         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj
166         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj
167         !
168         DO_2D_10_10
169            iku = mbku(ji,jj)
170            ikv = mbkv(ji,jj)
171            ze3wu  = e3w(ji+1,jj  ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
172            ze3wv  = e3w(ji  ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
173            IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
174            ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2
175            ENDIF
176            IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1
177            ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2
178            ENDIF
179         END_2D
180         CALL lbc_lnk_multi( 'zpshde', pgru , 'U', -1. , pgrv , 'V', -1. )   ! Lateral boundary conditions
181         !
182      END IF
183      !
184      IF( ln_timing )   CALL timing_stop( 'zps_hde')
185      !
186   END SUBROUTINE zps_hde
187
188
189   SUBROUTINE zps_hde_isf( kt, Kmm, kjpt, pta, pgtu, pgtv, pgtui, pgtvi,  &
190      &                          prd, pgru, pgrv, pgrui, pgrvi )
191      !!----------------------------------------------------------------------
192      !!                     ***  ROUTINE zps_hde_isf  ***
193      !!                   
194      !! ** Purpose :   Compute the horizontal derivative of T, S and rho
195      !!      at u- and v-points with a linear interpolation for z-coordinate
196      !!      with partial steps for top (ice shelf) and bottom.
197      !!
198      !! ** Method  :   In z-coord with partial steps, scale factors on last
199      !!      levels are different for each grid point, so that T, S and rd
200      !!      points are not at the same depth as in z-coord. To have horizontal
201      !!      gradients again, we interpolate T and S at the good depth :
202      !!      For the bottom case:
203      !!      Linear interpolation of T, S   
204      !!         Computation of di(tb) and dj(tb) by vertical interpolation:
205      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~
206      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~
207      !!         This formulation computes the two cases:
208      !!                 CASE 1                   CASE 2 
209      !!         k-1  ___ ___________   k-1   ___ ___________
210      !!                    Ti  T~                  T~  Ti+1
211      !!                  _____                        _____
212      !!         k        |   |Ti+1     k           Ti |   |
213      !!                  |   |____                ____|   |
214      !!              ___ |   |   |           ___  |   |   |
215      !!                 
216      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
217      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
218      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  )
219      !!          or
220      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
221      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
222      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
223      !!          Idem for di(s) and dj(s)         
224      !!
225      !!      For rho, we call eos which will compute rd~(t~,s~) at the right
226      !!      depth zh from interpolated T and S for the different formulations
227      !!      of the equation of state (eos).
228      !!      Gradient formulation for rho :
229      !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~
230      !!
231      !!      For the top case (ice shelf): As for the bottom case but upside down
232      !!
233      !! ** Action  : compute for top and bottom interfaces
234      !!              - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points
235      !!              - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points
236      !!----------------------------------------------------------------------
237      INTEGER                              , INTENT(in   )           ::  kt           ! ocean time-step index
238      INTEGER                              , INTENT(in   )           ::  Kmm          ! ocean time level index
239      INTEGER                              , INTENT(in   )           ::  kjpt         ! number of tracers
240      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta          ! 4D tracers fields
241      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv   ! hor. grad. of ptra at u- & v-pts
242      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF)
243      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd          ! 3D density anomaly fields
244      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv   ! hor. grad of prd at u- & v-pts (bottom)
245      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top)
246      !
247      INTEGER  ::   ji, jj, jn      ! Dummy loop indices
248      INTEGER  ::   iku, ikv, ikum1, ikvm1,ikup1, ikvp1   ! partial step level (ocean bottom level) at u- and v-points
249      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv             ! temporary scalars
250      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos
251      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !
252      !!----------------------------------------------------------------------
253      !
254      IF( ln_timing )   CALL timing_start( 'zps_hde_isf')
255      !
256      pgtu (:,:,:) = 0._wp   ;   pgtv (:,:,:) =0._wp
257      pgtui(:,:,:) = 0._wp   ;   pgtvi(:,:,:) =0._wp
258      zti  (:,:,:) = 0._wp   ;   ztj  (:,:,:) =0._wp
259      zhi  (:,:  ) = 0._wp   ;   zhj  (:,:  ) =0._wp
260      !
261      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
262         !
263         DO_2D_10_10
264
265            iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
266            ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
267            ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
268            ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
269            !
270            ! i- direction
271            IF( ze3wu >= 0._wp ) THEN      ! case 1
272               zmaxu =  ze3wu / e3w(ji+1,jj,iku,Kmm)
273               ! interpolated values of tracers
274               zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
275               ! gradient of  tracers
276               pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
277            ELSE                           ! case 2
278               zmaxu = -ze3wu / e3w(ji,jj,iku,Kmm)
279               ! interpolated values of tracers
280               zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
281               ! gradient of tracers
282               pgtu(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
283            ENDIF
284            !
285            ! j- direction
286            IF( ze3wv >= 0._wp ) THEN      ! case 1
287               zmaxv =  ze3wv / e3w(ji,jj+1,ikv,Kmm)
288               ! interpolated values of tracers
289               ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
290               ! gradient of tracers
291               pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
292            ELSE                           ! case 2
293               zmaxv =  -ze3wv / e3w(ji,jj,ikv,Kmm)
294               ! interpolated values of tracers
295               ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
296               ! gradient of tracers
297               pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
298            ENDIF
299
300         END_2D
301      END DO
302      !
303      CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1. , pgtv(:,:,:), 'V', -1. )   ! Lateral boundary cond.
304
305      ! horizontal derivative of density anomalies (rd)
306      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
307         pgru(:,:)=0.0_wp   ; pgrv(:,:)=0.0_wp ; 
308         !
309         DO_2D_10_10
310
311            iku = mbku(ji,jj)
312            ikv = mbkv(ji,jj)
313            ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
314            ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
315            !
316            IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku,Kmm)    ! i-direction: case 1
317            ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm)    ! -     -      case 2
318            ENDIF
319            IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv,Kmm)    ! j-direction: case 1
320            ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm)    ! -     -      case 2
321            ENDIF
322
323         END_2D
324
325         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
326         ! step and store it in  zri, zrj for each  case
327         CALL eos( zti, zhi, zri )
328         CALL eos( ztj, zhj, zrj )
329
330         DO_2D_10_10
331            iku = mbku(ji,jj)
332            ikv = mbkv(ji,jj)
333            ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
334            ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
335
336            IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
337            ELSE                        ;   pgru(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2
338            ENDIF
339            IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1
340            ELSE                        ;   pgrv(ji,jj) = ssvmask(ji,jj) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2
341            ENDIF
342
343         END_2D
344
345         CALL lbc_lnk_multi( 'zpshde', pgru , 'U', -1. , pgrv , 'V', -1. )   ! Lateral boundary conditions
346         !
347      END IF
348      !
349      !     !==  (ISH)  compute grui and gruvi  ==!
350      !
351      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            !
352         DO_2D_10_10
353            iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1
354            ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1
355            !
356            ! (ISF) case partial step top and bottom in adjacent cell in vertical
357            ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
358            ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
359            ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
360            ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
361            ze3wv  =  gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 
362
363            ! i- direction
364            IF( ze3wu >= 0._wp ) THEN      ! case 1
365               zmaxu = ze3wu / e3w(ji+1,jj,ikup1,Kmm)
366               ! interpolated values of tracers
367               zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) )
368               ! gradient of tracers
369               pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
370            ELSE                           ! case 2
371               zmaxu = - ze3wu / e3w(ji,jj,ikup1,Kmm)
372               ! interpolated values of tracers
373               zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) )
374               ! gradient of  tracers
375               pgtui(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
376            ENDIF
377            !
378            ! j- direction
379            IF( ze3wv >= 0._wp ) THEN      ! case 1
380               zmaxv =  ze3wv / e3w(ji,jj+1,ikvp1,Kmm)
381               ! interpolated values of tracers
382               ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) )
383               ! gradient of tracers
384               pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
385            ELSE                           ! case 2
386               zmaxv =  - ze3wv / e3w(ji,jj,ikvp1,Kmm)
387               ! interpolated values of tracers
388               ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) )
389               ! gradient of tracers
390               pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
391            ENDIF
392
393         END_2D
394         !
395      END DO
396      CALL lbc_lnk_multi( 'zpshde', pgtui(:,:,:), 'U', -1. , pgtvi(:,:,:), 'V', -1. )   ! Lateral boundary cond.
397
398      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
399         !
400         pgrui(:,:)  =0.0_wp; pgrvi(:,:)  =0.0_wp;
401         DO_2D_10_10
402
403            iku = miku(ji,jj)
404            ikv = mikv(ji,jj)
405            ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
406            ze3wv  =  gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 
407            !
408            IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = gdept(ji  ,jj,iku,Kmm)    ! i-direction: case 1
409            ELSE                        ;   zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm)    ! -     -      case 2
410            ENDIF
411
412            IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = gdept(ji,jj  ,ikv,Kmm)    ! j-direction: case 1
413            ELSE                        ;   zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm)    ! -     -      case 2
414            ENDIF
415
416         END_2D
417         !
418         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj
419         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj
420         !
421         DO_2D_10_10
422            iku = miku(ji,jj) 
423            ikv = mikv(ji,jj) 
424            ze3wu  =  gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
425            ze3wv  =  gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm) 
426
427            IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji  ,jj      ) - prd(ji,jj,iku) ) ! i: 1
428            ELSE                      ; pgrui(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj  ,iku) - zri(ji,jj    ) ) ! i: 2
429            ENDIF
430            IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( zrj(ji  ,jj      ) - prd(ji,jj,ikv) ) ! j: 1
431            ELSE                      ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( prd(ji  ,jj+1,ikv) - zrj(ji,jj    ) ) ! j: 2
432            ENDIF
433
434         END_2D
435         CALL lbc_lnk_multi( 'zpshde', pgrui, 'U', -1. , pgrvi, 'V', -1. )   ! Lateral boundary conditions
436         !
437      END IF 
438      !
439      IF( ln_timing )   CALL timing_stop( 'zps_hde_isf')
440      !
441   END SUBROUTINE zps_hde_isf
442
443   !!======================================================================
444END MODULE zpshde
Note: See TracBrowser for help on using the repository browser.