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

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/zpshde.F90 @ 2082

Last change on this file since 2082 was 2082, checked in by cetlod, 14 years ago

Improve the merge of TRA-TRC, see ticket #717

  • Property svn:eol-style set to native
  • 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
5   !!==============================================================================
6   !! History :
7   !!       OPA   8.5  !  2002-04  (A. Bozec)  Original code
8   !!             8.5  !  2002-08  (G. Madec E. Durand)  Optimization and Free form
9   !!             9.0  !  2004-03  (C. Ethe)  adapted for passive tracers
10   !!      NEMO   3.3  !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
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   !! * Modules used
18   USE dom_oce         ! ocean space domain variables
19   USE oce             ! ocean dynamics and tracers variables
20   USE phycst          ! physical constants
21   USE in_out_manager  ! I/O manager
22   USE eosbn2          ! ocean equation of state
23   USE lbclnk          ! lateral boundary conditions (or mpp link)
24
25   IMPLICIT NONE
26   PRIVATE
27
28   !! * Routine accessibility
29   PUBLIC zps_hde          ! routine called by step.F90
30   PUBLIC zps_hde_init     ! routine called by opa.F90
31
32   !! * module variables
33   INTEGER, DIMENSION(jpi,jpj) ::   &
34      mbatu, mbatv      ! bottom ocean level index at U- and V-points
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !!----------------------------------------------------------------------
41   !!  OPA 9.0 , LOCEAN-IPSL (2005)
42   !! $Id$
43   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE zps_hde( kt, kjpt, pta, pgtu, pgtv,   &
48                                 prd, pgru, pgrv    )
49      !!----------------------------------------------------------------------
50      !!                     ***  ROUTINE zps_hde  ***
51      !!                   
52      !! ** Purpose :   Compute the horizontal derivative of T, S and rd
53      !!      at u- and v-points with a linear interpolation for z-coordinate
54      !!      with partial steps.
55      !!
56      !! ** Method  :   In z-coord with partial steps, scale factors on last
57      !!      levels are different for each grid point, so that T, S and rd
58      !!      points are not at the same depth as in z-coord. To have horizontal
59      !!      gradients again, we interpolate T and S at the good depth :
60      !!      Linear interpolation of T, S   
61      !!         Computation of di(tb) and dj(tb) by vertical interpolation:
62      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~
63      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~
64      !!         This formulation computes the two cases:
65      !!                 CASE 1                   CASE 2 
66      !!         k-1  ___ ___________   k-1   ___ ___________
67      !!                    Ti  T~                  T~  Ti+1
68      !!                  _____                        _____
69      !!         k        |   |Ti+1     k           Ti |   |
70      !!                  |   |____                ____|   |
71      !!              ___ |   |   |           ___  |   |   |
72      !!                 
73      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
74      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
75      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  )
76      !!          or
77      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
78      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
79      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
80      !!          Idem for di(s) and dj(s)         
81      !!
82      !!      For rho, we call eos_insitu_2d which will compute rd~(t~,s~) at
83      !!      the good depth zh from interpolated T and S for the different
84      !!      formulation of the equation of state (eos).
85      !!      Gradient formulation for rho :
86      !!          di(rho) = rd~ - rd(i,j,k) or rd (i+1,j,k) - rd~
87      !!
88      !! ** Action  : - pgtu, pgtv: horizontal gradient of tracer at U/V-points
89      !!              - pgru, pgrv: horizontal gradient of rd if present at U/V-points
90      !!                and rd at V-points
91      !!----------------------------------------------------------------------
92      !! * Arguments
93      INTEGER                              , INTENT( in )           ::  kt    ! ocean time-step index
94      INTEGER                              , INTENT( in )           ::  kjpt  ! number of tracers
95      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT( in )           ::  pta   ! 4D active or passive tracers fields
96      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT( out)           ::  pgtu, pgtv  ! horizontal grad. of ptra u- and v-points
97      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT( in ), OPTIONAL ::  prd   ! 3D rd fields
98      REAL(wp), DIMENSION(jpi,jpj         ), INTENT( out), OPTIONAL ::  pgru, pgrv  ! horizontal grad. of prd u- and v-points
99      !! * Local declarations
100      INTEGER  ::   ji, jj, jn      ! Dummy loop indices
101      INTEGER  ::   iku, ikv        ! partial step level at u- and v-points
102      REAL(wp), DIMENSION(jpi,jpj,kjpt) ::   zti, ztj     ! interpolated value of tracer
103      REAL(wp), DIMENSION(jpi,jpj)      ::   zri, zrj     ! interpolated value of rd
104      REAL(wp), DIMENSION(jpi,jpj)      ::   zhi, zhj     ! depth of interpolation for eos2d
105      REAL(wp) ::  ze3wu, ze3wv, zmaxu, zmaxv  ! temporary scalars
106      !!----------------------------------------------------------------------
107
108
109      ! Interpolation of tracers at the last ocean level
110      DO jn = 1, kjpt
111# if defined key_vectopt_loop
112         jj = 1
113         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
114# else
115         DO jj = 1, jpjm1
116            DO ji = 1, jpim1
117# endif
118               ! last level
119               iku = mbatu(ji,jj)
120               ikv = mbatv(ji,jj)
121               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
122               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
123
124               ! i- direction
125               IF( ze3wu >= 0. ) THEN      ! case 1
126                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
127                  ! interpolated values of tracers
128                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,iku-1,jn) - pta(ji+1,jj,iku,jn) )
129                  ! gradient of  tracers
130                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
131               ELSE                        ! case 2
132                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
133                  ! interpolated values of tracers
134                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,iku-1,jn) - pta(ji,jj,iku,jn) )
135                  ! gradient of tracers
136                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
137               ENDIF
138
139               ! j- direction
140               IF( ze3wv >= 0. ) THEN      ! case 1
141                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
142                  ! interpolated values of tracers
143                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikv-1,jn) - pta(ji,jj+1,ikv,jn) )
144                  ! gradient of tracers
145                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
146               ELSE                        ! case 2
147                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
148                  ! interpolated values of tracers
149                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikv-1,jn) - pta(ji,jj,ikv,jn) )
150                  ! gradient of tracers
151                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
152               ENDIF
153# if ! defined key_vectopt_loop
154            END DO
155# endif
156         END DO
157
158         ! Lateral boundary conditions on each gradient
159         CALL lbc_lnk( pgtu(:,:,jn) , 'U', -1. )
160         CALL lbc_lnk( pgtv(:,:,jn) , 'V', -1. )
161
162      END DO
163
164      ! horizontal derivative of rd
165      IF( PRESENT( prd ) ) THEN
166         ! depth of the partial step level
167# if defined key_vectopt_loop
168         jj = 1
169         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
170# else
171         DO jj = 1, jpjm1
172            DO ji = 1, jpim1
173# endif
174               iku = mbatu(ji,jj)
175               ikv = mbatv(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. ) THEN   
179                  zhi(ji,jj) = fsdept(ji  ,jj,iku)
180               ELSE                   
181                  zhi(ji,jj) = fsdept(ji+1,jj,iku)
182               ENDIF
183               IF( ze3wv >= 0. ) THEN 
184                  zhj(ji,jj) = fsdept(ji,jj  ,ikv)
185               ELSE                   
186                  zhj(ji,jj) = fsdept(ji,jj+1,ikv)
187               ENDIF
188# if ! defined key_vectopt_loop
189            END DO
190# endif
191         END DO
192
193         ! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
194         ! step and store it in  zri, zrj for each  case
195         CALL eos( zti, zhi, zri )
196         CALL eos( ztj, zhj, zrj )
197
198         ! Gradient of density at the last level
199# if defined key_vectopt_loop
200         jj = 1
201         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
202# else
203         DO jj = 1, jpjm1
204            DO ji = 1, jpim1
205# endif
206               iku = mbatu(ji,jj)
207               ikv = mbatv(ji,jj)
208               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
209               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
210               IF( ze3wu >= 0. ) THEN    ! i-direction: case 1
211                  pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji,jj) - prd(ji,jj,iku) )
212               ELSE                      ! i-direction: case 2
213                  pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj) )
214               ENDIF
215               IF( ze3wv >= 0. ) THEN    ! j-direction: case 1
216                  pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj) - prd(ji,jj,ikv) )
217               ELSE                      ! j-direction: case 2
218                  pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj) )
219               ENDIF
220# if ! defined key_vectopt_loop
221            END DO
222# endif
223         END DO
224
225         ! Lateral boundary conditions on each gradient
226         CALL lbc_lnk( pgru , 'U', -1. )   ;   CALL lbc_lnk( pgrv , 'V', -1. )
227         !
228      END IF
229      !
230   END SUBROUTINE zps_hde
231
232   SUBROUTINE zps_hde_init
233      !!----------------------------------------------------------------------
234      !!                     ***  ROUTINE zps_hde_init  ***
235      !!
236      !! ** Purpose : Computation of bottom ocean level index at U- and V-points
237      !!                   
238      !!----------------------------------------------------------------------
239      !! * Local declarations
240      INTEGER ::   ji, jj           ! Dummy loop indices
241      REAL(wp), DIMENSION(jpi,jpj) :: zti, ztj     !  temporary arrays
242      !!----------------------------------------------------------------------
243
244      mbatu(:,:) = 0
245      mbatv(:,:) = 0
246      DO jj = 1, jpjm1
247         DO ji = 1, fs_jpim1   ! vector opt.
248            mbatu(ji,jj) = MAX( MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1, 2 )
249            mbatv(ji,jj) = MAX( MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1, 2 )
250         END DO
251      END DO
252      zti(:,:) = FLOAT( mbatu(:,:) )
253      ztj(:,:) = FLOAT( mbatv(:,:) )
254      ! lateral boundary conditions: T-point, sign unchanged
255      CALL lbc_lnk( zti , 'U', 1. )
256      CALL lbc_lnk( ztj , 'V', 1. )
257      mbatu(:,:) = MAX( INT( zti(:,:) ), 2 )
258      mbatv(:,:) = MAX( INT( ztj(:,:) ), 2 )
259
260   END SUBROUTINE zps_hde_init
261   !!======================================================================
262END MODULE zpshde
Note: See TracBrowser for help on using the repository browser.