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/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/zpshde.F90 @ 12340

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

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

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