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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 12.0 KB
Line 
1MODULE zpshde
2   !!======================================================================
3   !!                       ***  MODULE zpshde   ***
4   !! z-coordinate + partial step : Horizontal Derivative at ocean bottom level
5   !!======================================================================
6   !! History :  OPA  !  2002-04  (A. Bozec)  Original code
7   !!   NEMO     1.0  !  2002-08  (G. Madec E. Durand)  Optimization and Free form
8   !!             -   !  2004-03  (C. Ethe)  adapted for passive tracers
9   !!            3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
10   !!======================================================================
11   
12   !!----------------------------------------------------------------------
13   !!   zps_hde      :  Horizontal DErivative of T, S and rd at the last
14   !!                   ocean level (Z-coord. with Partial Steps)
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean: dynamics and tracers variables
17   USE dom_oce         ! domain: ocean variables
18   USE phycst          ! physical constants
19   USE eosbn2          ! ocean equation of state
20   USE in_out_manager  ! I/O manager
21   USE lbclnk          ! lateral boundary conditions (or mpp link)
22   USE lib_mpp         ! MPP library
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC   zps_hde    ! routine called by step.F90
28
29   !! * Control permutation of array indices
30#  include "oce_ftrans.h90"
31#  include "dom_oce_ftrans.h90"
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_insitu_2d which will compute rd~(t~,s~) at
79      !!      the good depth zh from interpolated T and S for the different
80      !!      formulation 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  : - 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      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
88      USE wrk_nemo, ONLY:   zri => wrk_2d_1 , zrj => wrk_2d_2   ! interpolated value of rd
89      USE wrk_nemo, ONLY:   zhi => wrk_2d_3 , zhj => wrk_2d_4   ! depth of interpolation for eos2d
90      !
91      INTEGER                              , INTENT(in   )           ::  kt          ! ocean time-step index
92      INTEGER                              , INTENT(in   )           ::  kjpt        ! number of tracers
93
94      !! DCSE_NEMO: This style defeats ftrans
95!     REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   )           ::  pta         ! 4D tracers fields
96!FTRANS pta :I :I :z :
97      REAL(wp), INTENT(in   )                      ::  pta(jpi,jpj,jpk,kjpt)         ! 4D tracers fields
98
99      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(  out)           ::  pgtu, pgtv  ! hor. grad. of ptra at u- & v-pts
100
101!     REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ), OPTIONAL ::  prd         ! 3D density anomaly fields
102!FTRANS prd :I :I :z
103      REAL(wp), INTENT(in   ), OPTIONAL            ::  prd(jpi,jpj,jpk)              ! 3D density anomaly fields
104
105      REAL(wp), DIMENSION(jpi,jpj         ), INTENT(  out), OPTIONAL ::  pgru, pgrv  ! hor. grad. of prd at u- & v-pts
106      !
107      INTEGER  ::   ji, jj, jn      ! Dummy loop indices
108      INTEGER  ::   iku, ikv, ikum1, ikvm1   ! partial step level (ocean bottom level) at u- and v-points
109      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars
110      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zti, ztj    ! interpolated value of tracer
111      !!----------------------------------------------------------------------
112
113      IF( wrk_in_use(2, 1,2,3,4) ) THEN
114         CALL ctl_stop('zps_hde: requested workspace arrays unavailable')  ;  RETURN
115      END IF
116
117      ! Allocate workspaces whose dimension is > jpk
118      ALLOCATE( zti(jpi,jpj,kjpt) )
119      ALLOCATE( ztj(jpi,jpj,kjpt) )
120
121      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
122         !
123# if defined key_vectopt_loop
124         jj = 1
125         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
126# else
127         DO jj = 1, jpjm1
128            DO ji = 1, jpim1
129# endif
130               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
131               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
132               ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
133               ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
134               !
135               ! i- direction
136               IF( ze3wu >= 0._wp ) THEN      ! case 1
137                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
138                  ! interpolated values of tracers
139                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
140                  ! gradient of  tracers
141#if defined key_z_first
142                  pgtu(ji,jj,jn) = umask_1(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
143#else
144                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
145#endif
146               ELSE                           ! case 2
147                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
148                  ! interpolated values of tracers
149                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
150                  ! gradient of tracers
151#if defined key_z_first
152                  pgtu(ji,jj,jn) = umask_1(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
153#else
154                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
155#endif
156               ENDIF
157               !
158               ! j- direction
159               IF( ze3wv >= 0._wp ) THEN      ! case 1
160                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
161                  ! interpolated values of tracers
162                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
163                  ! gradient of tracers
164#if defined key_z_first
165                  pgtv(ji,jj,jn) = vmask_1(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
166#else
167                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
168#endif
169               ELSE                           ! case 2
170                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
171                  ! interpolated values of tracers
172                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
173                  ! gradient of tracers
174#if defined key_z_first
175                  pgtv(ji,jj,jn) = vmask_1(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
176#else
177                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
178#endif
179               ENDIF
180# if ! defined key_vectopt_loop
181            END DO
182# endif
183         END DO
184         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
185         !
186      END DO
187
188      ! horizontal derivative of density anomalies (rd)
189      IF( PRESENT( prd ) ) THEN         ! depth of the partial step level
190# if defined key_vectopt_loop
191         jj = 1
192         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
193# else
194         DO jj = 1, jpjm1
195            DO ji = 1, jpim1
196# endif
197               iku = mbku(ji,jj)
198               ikv = mbkv(ji,jj)
199               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
200               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
201               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1
202               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2
203               ENDIF
204               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1
205               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2
206               ENDIF
207# if ! defined key_vectopt_loop
208            END DO
209# endif
210         END DO
211
212         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
213         ! step and store it in  zri, zrj for each  case
214         CALL eos( zti, zhi, zri )   ;   CALL eos( ztj, zhj, zrj )
215
216         ! Gradient of density at the last level
217# if defined key_vectopt_loop
218         jj = 1
219         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
220# else
221         DO jj = 1, jpjm1
222            DO ji = 1, jpim1
223# endif
224               iku = mbku(ji,jj)
225               ikv = mbkv(ji,jj)
226               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
227               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
228               IF( ze3wu >= 0._wp ) THEN   ;   pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji  ,jj) - prd(ji,jj,iku) )   ! i: 1
229               ELSE                        ;   pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) )   ! i: 2
230               ENDIF
231               IF( ze3wv >= 0._wp ) THEN   ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj  ) - prd(ji,jj,ikv) )   ! j: 1
232               ELSE                        ;   pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )   ! j: 2
233               ENDIF
234# if ! defined key_vectopt_loop
235            END DO
236# endif
237         END DO
238         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )   ! Lateral boundary conditions
239         !
240      END IF
241      !
242      IF( wrk_not_released(2, 1,2,3,4) )   CALL ctl_stop('zps_hde: failed to release workspace arrays')
243      !
244      DEALLOCATE( zti )
245      DEALLOCATE( ztj )
246      !
247   END SUBROUTINE zps_hde
248
249   !!======================================================================
250END MODULE zpshde
Note: See TracBrowser for help on using the repository browser.