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

source: trunk/NEMO/OPA_SRC/TRA/zpshde.F90 @ 87

Last change on this file since 87 was 87, checked in by opalod, 20 years ago

CT : UPDATE057 : # General syntax, alignement, comments corrections

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