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

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
  • Property svn:keywords set to Id
File size: 32.6 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 "domzgr_substitute.h90"
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   &
44      &                          prd, pgru, pgrv    )
45      !!----------------------------------------------------------------------
46      !!                     ***  ROUTINE zps_hde  ***
47      !!                   
48      !! ** Purpose :   Compute the horizontal derivative of T, S and rho
49      !!      at u- and v-points with a linear interpolation for z-coordinate
50      !!      with partial steps.
51      !!
52      !! ** Method  :   In z-coord with partial steps, scale factors on last
53      !!      levels are different for each grid point, so that T, S and rd
54      !!      points are not at the same depth as in z-coord. To have horizontal
55      !!      gradients again, we interpolate T and S at the good depth :
56      !!      Linear interpolation of T, S   
57      !!         Computation of di(tb) and dj(tb) by vertical interpolation:
58      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~
59      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~
60      !!         This formulation computes the two cases:
61      !!                 CASE 1                   CASE 2 
62      !!         k-1  ___ ___________   k-1   ___ ___________
63      !!                    Ti  T~                  T~  Ti+1
64      !!                  _____                        _____
65      !!         k        |   |Ti+1     k           Ti |   |
66      !!                  |   |____                ____|   |
67      !!              ___ |   |   |           ___  |   |   |
68      !!                 
69      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
70      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
71      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  )
72      !!          or
73      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
74      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
75      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
76      !!          Idem for di(s) and dj(s)         
77      !!
78      !!      For rho, we call eos which will compute rd~(t~,s~) at the right
79      !!      depth zh from interpolated T and S for the different formulations
80      !!      of the equation of state (eos).
81      !!      Gradient formulation for rho :
82      !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~
83      !!
84      !! ** Action  : compute for top interfaces
85      !!              - pgtu, pgtv: horizontal gradient of tracer at u- & v-points
86      !!              - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points
87      !!----------------------------------------------------------------------
88      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step 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( nn_timing == 1 )   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 jj = 1, jpjm1
110            DO ji = 1, jpim1
111               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
112               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
113               ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
114               ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
115               !
116               ! i- direction
117               IF( ze3wu >= 0._wp ) THEN      ! case 1
118                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
119                  ! interpolated values of tracers
120                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
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 / fse3w(ji,jj,iku)
125                  ! interpolated values of tracers
126                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
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 / fse3w(ji,jj+1,ikv)
134                  ! interpolated values of tracers
135                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
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 / fse3w(ji,jj,ikv)
140                  ! interpolated values of tracers
141                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
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 DO
146         END DO
147         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
148         !
149      END DO
150      !               
151      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
152         pgru(:,:) = 0._wp
153         pgrv(:,:) = 0._wp                ! depth of the partial step level
154         DO jj = 1, jpjm1
155            DO ji = 1, jpim1
156               iku = mbku(ji,jj)
157               ikv = mbkv(ji,jj)
158               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
159               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
160               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1
161               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2
162               ENDIF
163               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1
164               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2
165               ENDIF
166            END DO
167         END DO
168         !
169         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj
170         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj
171         !
172         DO jj = 1, jpjm1                 ! Gradient of density at the last level
173            DO ji = 1, jpim1
174               iku = mbku(ji,jj)
175               ikv = mbkv(ji,jj)
176               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
177               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
178               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj    ) - prd(ji,jj,iku) )   ! i: 1
179               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj    ) )   ! i: 2
180               ENDIF
181               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj      ) - prd(ji,jj,ikv) )   ! j: 1
182               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj    ) )   ! j: 2
183               ENDIF
184            END DO
185         END DO
186         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions
187         !
188      END IF
189      !
190      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde')
191      !
192   END SUBROUTINE zps_hde
193
194
195   SUBROUTINE zps_hde_isf( kt, kjpt, pta, pgtu , pgtv , pgtui, pgtvi,                                   &
196      &                              prd, pgru , pgrv , pmru , pmrv , pgzu , pgzv , pge3ru , pge3rv ,   &
197      &                                   pgrui, pgrvi, pmrui, pmrvi, pgzui, pgzvi, pge3rui, pge3rvi )
198      !!----------------------------------------------------------------------
199      !!                     ***  ROUTINE zps_hde  ***
200      !!                   
201      !! ** Purpose :   Compute the horizontal derivative of T, S and rho
202      !!      at u- and v-points with a linear interpolation for z-coordinate
203      !!      with partial steps.
204      !!
205      !! ** Method  :   In z-coord with partial steps, scale factors on last
206      !!      levels are different for each grid point, so that T, S and rd
207      !!      points are not at the same depth as in z-coord. To have horizontal
208      !!      gradients again, we interpolate T and S at the good depth :
209      !!      Linear interpolation of T, S   
210      !!         Computation of di(tb) and dj(tb) by vertical interpolation:
211      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~
212      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~
213      !!         This formulation computes the two cases:
214      !!                 CASE 1                   CASE 2 
215      !!         k-1  ___ ___________   k-1   ___ ___________
216      !!                    Ti  T~                  T~  Ti+1
217      !!                  _____                        _____
218      !!         k        |   |Ti+1     k           Ti |   |
219      !!                  |   |____                ____|   |
220      !!              ___ |   |   |           ___  |   |   |
221      !!                 
222      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
223      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
224      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  )
225      !!          or
226      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
227      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
228      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
229      !!          Idem for di(s) and dj(s)         
230      !!
231      !!      For rho, we call eos which will compute rd~(t~,s~) at the right
232      !!      depth zh from interpolated T and S for the different formulations
233      !!      of the equation of state (eos).
234      !!      Gradient formulation for rho :
235      !!          di(rho) = rd~ - rd(i,j,k)   or   rd(i+1,j,k) - rd~
236      !!
237      !! ** Action  : compute for top and bottom interfaces
238      !!              - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points
239      !!              - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points
240      !!              - pmru, pmrv, pmrui, pmrvi: horizontal sum of rho at u- & v- point (used in dynhpg with vvl)
241      !!              - pgzu, pgzv, pgzui, pgzvi: horizontal gradient of z at u- and v- point (used in dynhpg with vvl)
242      !!              - pge3ru, pge3rv, pge3rui, pge3rvi: horizontal gradient of rho weighted by local e3w at u- & v-points
243      !!----------------------------------------------------------------------
244      INTEGER                              , INTENT(in   )           ::  kt                ! ocean time-step index
245      INTEGER                              , INTENT(in   )           ::  kjpt              ! number of tracers
246      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta               ! 4D tracers fields
247      !                                                              !!  u-point ! v-point !
248      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu    , pgtv    ! bottom GRADh( ptra ) 
249      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtui   , pgtvi   ! top    GRADh( ptra )
250      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd               ! 3D density anomaly fields
251      !                                                              !!  u-point ! v-point !
252      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru    , pgrv    ! bottom GRADh( prd  )
253      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmru    , pmrv    ! bottom SUM  ( prd  )
254      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzu    , pgzv    ! bottom GRADh( z    )
255      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3ru  , pge3rv  ! bottom GRADh( prd  ) weighted by e3w
256      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgrui   , pgrvi   ! top    GRADh( prd  )
257      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pmrui   , pmrvi   ! top    SUM  ( prd  )
258      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgzui   , pgzvi   ! top    GRADh( z    )
259      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pge3rui , pge3rvi ! top    GRADh( prd  ) weighted by e3w
260      !
261      INTEGER  ::   ji, jj, jn      ! Dummy loop indices
262      INTEGER  ::   iku, ikv, ikum1, ikvm1,ikup1, ikvp1   ! partial step level (ocean bottom level) at u- and v-points
263      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv, zdzwu, zdzwv, zdzwuip1, zdzwvjp1  ! temporary scalars
264      REAL(wp), DIMENSION(jpi,jpj)      ::  zri, zrj, zhi, zhj   ! NB: 3rd dim=1 to use eos
265      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::  zti, ztj             !
266      !!----------------------------------------------------------------------
267      !
268      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_isf')
269      !
270      pgtu (:,:,:) = 0._wp   ;   pgtv (:,:,:) =0._wp
271      pgtui(:,:,:) = 0._wp   ;   pgtvi(:,:,:) =0._wp
272      zti  (:,:,:) = 0._wp   ;   ztj  (:,:,:) =0._wp
273      zhi  (:,:  ) = 0._wp   ;   zhj  (:,:  ) =0._wp
274      !
275      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
276         !
277         DO jj = 1, jpjm1
278            DO ji = 1, jpim1
279               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
280               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
281               ! (ISF) case partial step top and bottom in adjacent cell in vertical
282               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
283               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
284               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
285               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
286               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
287               !
288               ! i- direction
289               IF( ze3wu >= 0._wp ) THEN      ! case 1
290                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
291                  ! interpolated values of tracers
292                  zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
293                  ! gradient of  tracers
294                  pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
295               ELSE                           ! case 2
296                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
297                  ! interpolated values of tracers
298                  zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
299                  ! gradient of tracers
300                  pgtu(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
301               ENDIF
302               !
303               ! j- direction
304               IF( ze3wv >= 0._wp ) THEN      ! case 1
305                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
306                  ! interpolated values of tracers
307                  ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
308                  ! gradient of tracers
309                  pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
310               ELSE                           ! case 2
311                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
312                  ! interpolated values of tracers
313                  ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
314                  ! gradient of tracers
315                  pgtv(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
316               ENDIF
317            END DO
318         END DO
319         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
320         !
321      END DO
322
323      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
324         !
325         pgru  (:,:)=0._wp   ;   pgrv  (:,:) = 0._wp
326         pgzu  (:,:)=0._wp   ;   pgzv  (:,:) = 0._wp 
327         pmru  (:,:)=0._wp   ;   pmru  (:,:) = 0._wp 
328         pge3ru(:,:)=0._wp   ;   pge3rv(:,:) = 0._wp 
329         !
330         DO jj = 1, jpjm1                 ! depth of the partial step level
331            DO ji = 1, jpim1
332               iku = mbku(ji,jj)
333               ikv = mbkv(ji,jj)
334               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
335               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
336               !
337               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) - ze3wu     ! i-direction: case 1
338               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) + ze3wu    ! -     -      case 2
339               ENDIF
340               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) - ze3wv    ! j-direction: case 1
341               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) + ze3wv    ! -     -      case 2
342               ENDIF
343            END DO
344         END DO
345         !
346         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj
347         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj
348
349         DO jj = 1, jpjm1                 ! Gradient of density at the last level
350            DO ji = 1, jpim1
351               iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
352               ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 )    ! last and before last ocean level at u- & v-points
353               ze3wu  = (gdept_0(ji+1,jj,iku) - gdepw_0(ji+1,jj,iku)) - (gdept_0(ji,jj,iku) - gdepw_0(ji,jj,iku))
354               ze3wv  = (gdept_0(ji,jj+1,ikv) - gdepw_0(ji,jj+1,ikv)) - (gdept_0(ji,jj,ikv) - gdepw_0(ji,jj,ikv))
355               IF( ze3wu >= 0._wp ) THEN
356                  pgzu(ji,jj) = (fsde3w(ji+1,jj,iku) - ze3wu) - fsde3w(ji,jj,iku)
357                  pgru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1
358                  pmru(ji,jj) = umask(ji,jj,iku) * ( zri(ji  ,jj) + prd(ji,jj,iku) )   ! i: 1
359                  pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  &
360                                * ( (fse3w(ji+1,jj,iku) - ze3wu )* ( zri(ji  ,jj    ) + prd(ji+1,jj,ikum1) + 2._wp) &
361                                   - fse3w(ji  ,jj,iku)          * ( prd(ji  ,jj,iku) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2
362               ELSE 
363                  pgzu(ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) + ze3wu)
364                  pgru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2
365                  pmru(ji,jj) = umask(ji,jj,iku) * ( prd(ji+1,jj,iku) + zri(ji,jj) )   ! i: 2
366                  pge3ru(ji,jj) = umask(ji,jj,iku)                                                                  &
367                                * (  fse3w(ji+1,jj,iku)          * ( prd(ji+1,jj,iku) + prd(ji+1,jj,ikum1) + 2._wp) &
368                                   -(fse3w(ji  ,jj,iku) + ze3wu) * ( zri(ji  ,jj    ) + prd(ji  ,jj,ikum1) + 2._wp) )  ! j: 2
369               ENDIF
370               IF( ze3wv >= 0._wp ) THEN
371                  pgzv(ji,jj) = (fsde3w(ji,jj+1,ikv) - ze3wv) - fsde3w(ji,jj,ikv) 
372                  pgrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1
373                  pmrv(ji,jj) = vmask(ji,jj,ikv) * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )   ! j: 1
374                  pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  &
375                                * ( (fse3w(ji,jj+1,ikv) - ze3wv )* ( zrj(ji,jj      ) + prd(ji,jj+1,ikvm1) + 2._wp) &
376                                   - fse3w(ji,jj  ,ikv)          * ( prd(ji,jj  ,ikv) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2
377               ELSE
378                  pgzv(ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) + ze3wv)
379                  pgrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2
380                  pmrv(ji,jj) = vmask(ji,jj,ikv) * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )   ! j: 2
381                  pge3rv(ji,jj) = vmask(ji,jj,ikv)                                                                  &
382                                * (  fse3w(ji,jj+1,ikv)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikvm1) + 2._wp) &
383                                   -(fse3w(ji,jj  ,ikv) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikvm1) + 2._wp) )  ! j: 2
384               ENDIF
385            END DO
386         END DO
387         CALL lbc_lnk( pgru   , 'U', -1. )   ;   CALL lbc_lnk( pgrv   , 'V', -1. )   ! Lateral boundary conditions
388         CALL lbc_lnk( pmru   , 'U',  1. )   ;   CALL lbc_lnk( pmrv   , 'V',  1. )   ! Lateral boundary conditions
389         CALL lbc_lnk( pgzu   , 'U', -1. )   ;   CALL lbc_lnk( pgzv   , 'V', -1. )   ! Lateral boundary conditions
390         CALL lbc_lnk( pge3ru , 'U', -1. )   ;   CALL lbc_lnk( pge3rv , 'V', -1. )   ! Lateral boundary conditions
391         !
392      END IF
393      !
394      !     !==  (ISH)  compute grui and gruvi  ==!
395      !
396      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!            !
397         DO jj = 1, jpjm1
398            DO ji = 1, jpim1
399               iku = miku(ji,jj)   ;  ikup1 = miku(ji,jj) + 1
400               ikv = mikv(ji,jj)   ;  ikvp1 = mikv(ji,jj) + 1
401               !
402               ! (ISF) case partial step top and bottom in adjacent cell in vertical
403               ! cannot used e3w because if 2 cell water column, we have ps at top and bottom
404               ! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
405               ! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
406               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)) 
407               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))
408               ! i- direction
409               IF( ze3wu >= 0._wp ) THEN      ! case 1
410                  zmaxu = ze3wu / fse3w(ji+1,jj,iku+1)
411                  ! interpolated values of tracers
412                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku+1,jn) - pta(ji+1,jj,iku,jn) )
413                  ! gradient of tracers
414                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
415               ELSE                           ! case 2
416                  zmaxu = - ze3wu / fse3w(ji,jj,iku+1)
417                  ! interpolated values of tracers
418                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku+1,jn) - pta(ji,jj,iku,jn) )
419                  ! gradient of  tracers
420                  pgtui(ji,jj,jn) = umask(ji,jj,iku) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
421               ENDIF
422               !
423               ! j- direction
424               IF( ze3wv >= 0._wp ) THEN      ! case 1
425                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv+1)
426                  ! interpolated values of tracers
427                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv+1,jn) - pta(ji,jj+1,ikv,jn) )
428                  ! gradient of tracers
429                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
430               ELSE                           ! case 2
431                  zmaxv =  - ze3wv / fse3w(ji,jj,ikv+1)
432                  ! interpolated values of tracers
433                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv+1,jn) - pta(ji,jj,ikv,jn) )
434                  ! gradient of tracers
435                  pgtvi(ji,jj,jn) = vmask(ji,jj,ikv) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
436               ENDIF
437            END DO!!
438         END DO!!
439         CALL lbc_lnk( pgtui(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtvi(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
440         !
441      END DO
442
443      IF( PRESENT( prd ) ) THEN    !==  horizontal derivative of density anomalies (rd)  ==!    (optional part)
444         !
445         pgrui(:,:)  =0.0_wp ; pgrvi(:,:)  =0.0_wp ;
446         pgzui(:,:)  =0.0_wp ; pgzvi(:,:)  =0.0_wp ;
447         pmrui(:,:)  =0.0_wp ; pmrui(:,:)  =0.0_wp ;
448         pge3rui(:,:)=0.0_wp ; pge3rvi(:,:)=0.0_wp ;
449         !
450         DO jj = 1, jpjm1        ! depth of the partial step level
451            DO ji = 1, jpim1
452               iku = miku(ji,jj)
453               ikv = mikv(ji,jj)
454               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))
455               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))
456               !
457               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji+1,jj,iku) + ze3wu    ! i-direction: case 1
458               ELSE                        ;   zhi(ji,jj) = fsdept(ji  ,jj,iku) - ze3wu    ! -     -      case 2
459               ENDIF
460               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv) + ze3wv    ! j-direction: case 1
461               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv) - ze3wv    ! -     -      case 2
462               ENDIF
463            END DO
464         END DO
465         !
466         CALL eos( zti, zhi, zri )        ! interpolated density from zti, ztj
467         CALL eos( ztj, zhj, zrj )        ! at the partial step depth output in  zri, zrj
468         !
469         DO jj = 1, jpjm1                 ! Gradient of density at the last level
470            DO ji = 1, jpim1
471               iku = miku(ji,jj) ; ikup1 = miku(ji,jj) + 1
472               ikv = mikv(ji,jj) ; ikvp1 = mikv(ji,jj) + 1
473               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))
474               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))
475               IF( ze3wu >= 0._wp ) THEN
476                 pgzui  (ji,jj) = (fsde3w(ji+1,jj,iku) + ze3wu) - fsde3w(ji,jj,iku)
477                 pgrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) - prd(ji,jj,iku) )          ! i: 1
478                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( zri(ji,jj) + prd(ji,jj,iku) )          ! i: 1
479                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                  &
480                    &           * ( (fse3w(ji+1,jj,iku+1) - ze3wu) * (zri(ji,jj    ) + prd(ji+1,jj,iku+1) + 2._wp)   &
481                    &              - fse3w(ji  ,jj,iku+1)          * (prd(ji,jj,iku) + prd(ji  ,jj,iku+1) + 2._wp)   ) ! i: 1
482               ELSE
483                 pgzui  (ji,jj) = fsde3w(ji+1,jj,iku) - (fsde3w(ji,jj,iku) - ze3wu)
484                 pgrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) - zri(ji,jj) )      ! i: 2
485                 pmrui  (ji,jj) = umask(ji,jj,iku)   * ( prd(ji+1,jj,iku) + zri(ji,jj) )      ! i: 2
486                 pge3rui(ji,jj) = umask(ji,jj,iku+1)                                                                   &
487                    &           * (  fse3w(ji+1,jj,iku+1)          * (prd(ji+1,jj,iku) + prd(ji+1,jj,iku+1) + 2._wp)  &
488                    &              -(fse3w(ji  ,jj,iku+1) + ze3wu) * (zri(ji,jj      ) + prd(ji  ,jj,iku+1) + 2._wp)  )     ! i: 2
489               ENDIF
490               IF( ze3wv >= 0._wp ) THEN
491                 pgzvi  (ji,jj) = (fsde3w(ji,jj+1,ikv) + ze3wv) - fsde3w(ji,jj,ikv) 
492                 pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )        ! j: 1
493                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( zrj(ji,jj  ) + prd(ji,jj,ikv) )        ! j: 1
494                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                  & 
495                     &           * ( (fse3w(ji,jj+1,ikv+1) - ze3wv) * ( zrj(ji,jj    ) + prd(ji,jj+1,ikv+1) + 2._wp)  &
496                     &                - fse3w(ji,jj  ,ikv+1)          * ( prd(ji,jj,ikv) + prd(ji,jj  ,ikv+1) + 2._wp)  ) ! j: 1
497                                  ! + 2 due to the formulation in density and not in anomalie in hpg sco
498               ELSE
499                 pgzvi  (ji,jj) = fsde3w(ji,jj+1,ikv) - (fsde3w(ji,jj,ikv) - ze3wv)
500                 pgrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )     ! j: 2
501                 pmrvi  (ji,jj) = vmask(ji,jj,ikv)   * ( prd(ji,jj+1,ikv) + zrj(ji,jj) )     ! j: 2
502                 pge3rvi(ji,jj) = vmask(ji,jj,ikv+1)                                                                   &
503                    &           * (  fse3w(ji,jj+1,ikv+1)          * ( prd(ji,jj+1,ikv) + prd(ji,jj+1,ikv+1) + 2._wp) &
504                    &              -(fse3w(ji,jj  ,ikv+1) + ze3wv) * ( zrj(ji,jj      ) + prd(ji,jj  ,ikv+1) + 2._wp) )  ! j: 2
505               ENDIF
506            END DO
507         END DO
508         CALL lbc_lnk( pgrui   , 'U', -1. )   ;   CALL lbc_lnk( pgrvi   , 'V', -1. )   ! Lateral boundary conditions
509         CALL lbc_lnk( pmrui   , 'U',  1. )   ;   CALL lbc_lnk( pmrvi   , 'V',  1. )   ! Lateral boundary conditions
510         CALL lbc_lnk( pgzui   , 'U', -1. )   ;   CALL lbc_lnk( pgzvi   , 'V', -1. )   ! Lateral boundary conditions
511         CALL lbc_lnk( pge3rui , 'U', -1. )   ;   CALL lbc_lnk( pge3rvi , 'V', -1. )   ! Lateral boundary conditions
512         !
513      END IF 
514      !
515      IF( nn_timing == 1 )   CALL timing_stop( 'zps_hde_isf')
516      !
517   END SUBROUTINE zps_hde_isf
518   !!======================================================================
519END MODULE zpshde
Note: See TracBrowser for help on using the repository browser.