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_tam.F90 in branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/TRA – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/TRA/zpshde_tam.F90 @ 3659

Last change on this file since 3659 was 3659, checked in by pabouttier, 12 years ago

Use ikum1 and ikvm1 in zpshde_tam - See ticket #1023
  • Property svn:executable set to *
File size: 36.2 KB
Line 
1MODULE zpshde_tam
2#ifdef key_tam
3   !!==============================================================================
4   !!                       ***  MODULE zpshde_tam   ***
5   !! z-coordinate - partial step : Horizontal Derivative
6   !!                               Tangent and Adjoint Module
7   !!==============================================================================
8
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 par_kind
15   USE par_oce
16   USE oce
17   USE oce_tam
18   USE dom_oce
19   USE in_out_manager
20   USE eosbn2_tam
21   USE lbclnk
22   USE lbclnk_tam
23   USE gridrandom
24   USE dotprodfld
25   USE tstool_tam
26   USE lib_mpp
27   USE lib_mpp_tam
28   USE wrk_nemo
29   USE timing
30
31   IMPLICIT NONE
32   PRIVATE
33
34   !! * Routine accessibility
35   PUBLIC zps_hde_tan      ! routine called by step_tam.F90
36   PUBLIC zps_hde_adj      ! routine called by step_tam.F90
37   !PUBLIC zps_hde_adj_tst  ! routine called by tst.F90
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42
43CONTAINS
44
45   SUBROUTINE zps_hde_tan ( kt, kjpt, pta,            &
46      &                      pta_tl,  prd_tl,    &
47      &                      pgtu_tl, pgru_tl,   &
48      &                      pgtv_tl, pgrv_tl     )
49      !!----------------------------------------------------------------------
50      !!                     ***  ROUTINE zps_hde_tan  ***
51      !!
52      !! ** Purpose of the direct routine:
53      !!      Compute the horizontal derivative of T, S and rd
54      !!      at u- and v-points with a linear interpolation for z-coordinate
55      !!      with partial steps.
56      !!
57      !! ** Method of the direct routine:
58      !!      In z-coord with partial steps, scale factors on last
59      !!      levels are different for each grid point, so that T, S and rd
60      !!      points are not at the same depth as in z-coord. To have horizontal
61      !!      gradients again, we interpolate T and S at the good depth :
62      !!      Linear interpolation of T, S
63      !!         Computation of di(tb) and dj(tb) by vertical interpolation:
64      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~
65      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~
66      !!         This formulation computes the two cases:
67      !!                 CASE 1                   CASE 2
68      !!         k-1  ___ ___________   k-1   ___ ___________
69      !!                    Ti  T~                  T~  Ti+1
70      !!                  _____                        _____
71      !!         k        |   |Ti+1     k           Ti |   |
72      !!                  |   |____                ____|   |
73      !!              ___ |   |   |           ___  |   |   |
74      !!
75      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
76      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
77      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  )
78      !!          or
79      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
80      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
81      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
82      !!          Idem for di(s) and dj(s)
83      !!
84      !!      For rho, we call eos_insitu_2d which will compute rd~(t~,s~) at
85      !!      the good depth zh from interpolated T and S for the different
86      !!      formulation of the equation of state (eos).
87      !!      Gradient formulation for rho :
88      !!          di(rho) = rd~ - rd(i,j,k) or rd (i+1,j,k) - rd~
89      !!
90      !! ** Action  : - pgtu, pgsu, pgru: horizontal gradient of T, S
91      !!                and rd at U-points
92      !!              - pgtv, pgsv, pgrv: horizontal gradient of T, S
93      !!                and rd at V-points
94      !!
95      !! History of the direct routine:
96      !!   8.5  !  02-04  (A. Bozec)  Original code
97      !!   8.5  !  02-08  (G. Madec E. Durand)  Optimization and Free form
98      !! History of the TAM routine:
99      !!   9.0  !  08-06 (A. Vidard) Skeleton
100      !!        !  08-06 (A. Vidard) tangent of the 02-08 version
101      !!----------------------------------------------------------------------
102      !! * Arguments
103      INTEGER, INTENT( in ) ::   kt, kjpt          ! ocean time-step index
104      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT( in ) ::   &
105         pta, pta_tl            ! 3D T, S and rd direct fields
106      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ), OPTIONAL ::   &
107         prd_tl
108      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ), OPTIONAL ::   &
109         pgtu_tl, pgtv_tl                                    ! 3D T, S and rd tangent fields
110      REAL(wp), DIMENSION(jpi,jpj), INTENT( out ), OPTIONAL ::   &
111         pgru_tl,                               &  ! horizontal grad. of T, S and rd at u-
112         pgrv_tl                                   ! and v-points of the partial step level
113      !! * Local declarations
114      INTEGER ::   ji, jj,jk, jn,            &        ! Dummy loop indices
115               &    iku,ikv, ikum1, ikvm1          ! partial step level at u- and v-points
116      REAL(wp), POINTER, DIMENSION(:,:) ::   &
117         zri, zrj,               &  ! and rd
118         zritl, zrjtl,           &  ! and rdtl
119         zhi, zhj                 ! depth of interpolation for eos2d
120      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zti, ztj, ztitl, ztjtl    ! interpolated value of tracer
121      REAL(wp) ::   &
122         ze3wu, ze3wv,           &  ! temporary scalars
123         zmaxu, zmaxu2,         &  !    "         "
124         zmaxv, zmaxv2             !    "         "
125      !!---------------------------------------------------------------------
126      !
127      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_tan')
128      !
129      CALL wrk_alloc( jpi, jpj,       zri, zrj, zhi, zhj, zritl, zrjtl )
130      CALL wrk_alloc( jpi, jpj, kjpt, zti, ztj, ztitl, ztjtl           )
131      !
132      DO jn = 1, kjpt
133         ! Interpolation of T and S at the last ocean level
134# if defined key_vectopt_loop
135            jj = 1
136            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
137# else
138         DO jj = 1, jpjm1
139            DO ji = 1, jpim1
140# endif
141               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
142               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
143               ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
144               ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
145               !
146               ! i- direction
147               IF( ze3wu >= 0. ) THEN      ! case 1
148                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
149                  ! interpolated values of T and S
150                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn)            &
151                     &       + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
152                  ztitl(ji,jj,jn) = pta_tl(ji+1,jj,iku,jn)       &
153                     &         + zmaxu * ( pta_tl(ji+1,jj,ikum1,jn) - pta_tl(ji+1,jj,iku,jn) )
154                  ! gradient of T and S
155                  pgtu_tl(ji,jj,jn) = umask(ji,jj,1) * ( ztitl(ji,jj,jn) - pta_tl(ji,jj,iku,jn) )
156               ELSE                        ! case 2
157                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
158                  ! interpolated values of T and S
159                  zti(ji,jj,jn) = pta(ji,jj,iku,jn)              &
160                     &       + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
161                  ! interpolated values of T and S
162                  ztitl(ji,jj,jn) = pta_tl(ji,jj,iku,jn)         &
163                     &         + zmaxu * ( pta_tl(ji,jj,ikum1,jn) - pta_tl(ji,jj,iku,jn) )
164                  ! gradient of T and S
165                  pgtu_tl(ji,jj,jn) = umask(ji,jj,1) * ( pta_tl(ji+1,jj,iku,jn) - ztitl (ji,jj,jn) )
166               ENDIF
167               ! j- direction
168               IF( ze3wv >= 0. ) THEN      ! case 1
169                  ! interpolated values of direct T and S
170                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
171                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn)            &
172                     &       + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
173                  ! interpolated values of tangent T and S
174                  ztjtl(ji,jj,jn) = pta_tl(ji,jj+1,ikv,jn)       &
175                     &         + zmaxv * ( pta_tl(ji,jj+1,ikvm1,jn) - pta_tl(ji,jj+1,ikv,jn) )
176                  ! gradient of T and S
177                  pgtv_tl(ji,jj,jn) = vmask(ji,jj,1) * ( ztjtl(ji,jj,jn) - pta_tl(ji,jj,ikv,jn) )
178               ELSE                        ! case 2
179                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
180                  ! interpolated values of T and S
181                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn)       &
182                     &       + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
183                  ! interpolated values of T and S
184                  ztjtl(ji,jj,jn) = pta_tl(ji,jj,ikv,jn)  &
185                     &         + zmaxv * ( pta_tl(ji,jj,ikvm1,jn) - pta_tl(ji,jj,ikv,jn) )
186                  ! gradient of T and S
187                  pgtv_tl(ji,jj,jn) = vmask(ji,jj,1) * ( pta_tl(ji,jj+1,ikv,jn) - ztjtl(ji,jj,jn) )
188               ENDIF
189# if ! defined key_vectopt_loop
190            END DO
191# endif
192         END DO
193         CALL lbc_lnk( pgtu_tl(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv_tl(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
194         !
195      END DO
196      !
197      ! horizontal derivative of density anomalies (rd)
198      IF( PRESENT( prd_tl ) ) THEN         ! depth of the partial step 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 = mbku(ji,jj)
207               ikv = mbkv(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._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1
211               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2
212               ENDIF
213               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1
214               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2
215               ENDIF
216# if ! defined key_vectopt_loop
217            END DO
218# endif
219         END DO
220         ! Compute interpolated rd from zti, zsi, ztj, zsj for the 2 cases at the depth of the partial
221         ! step and store it in  zri, zrj for each  case
222         CALL eos_tan( zti, zhi, ztitl, zritl )
223         CALL eos_tan( ztj, zhj, ztjtl, zrjtl )
224
225
226         ! Gradient of density at the last level
227# if defined key_vectopt_loop
228            jj = 1
229            DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
230# else
231         DO jj = 1, jpjm1
232            DO ji = 1, jpim1
233# endif
234               iku = mbku(ji,jj)
235               ikv = mbkv(ji,jj)
236               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
237               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
238               IF( ze3wu >= 0. ) THEN    ! i-direction: case 1
239                  pgru_tl(ji,jj) = umask(ji,jj,1) * ( zritl(ji,jj) - prd_tl(ji,jj,iku) )
240               ELSE                      ! i-direction: case 2
241                  pgru_tl(ji,jj) = umask(ji,jj,1) * ( prd_tl(ji+1,jj,iku) - zritl(ji,jj) )
242               ENDIF
243               IF( ze3wv >= 0. ) THEN    ! j-direction: case 1
244                  pgrv_tl(ji,jj) = vmask(ji,jj,1) * ( zrjtl(ji,jj) - prd_tl(ji,jj,ikv) )
245               ELSE                      ! j-direction: case 2
246                  pgrv_tl(ji,jj) = vmask(ji,jj,1) * ( prd_tl(ji,jj+1,ikv) - zrjtl(ji,jj) )
247               ENDIF
248# if ! defined key_vectopt_loop
249            END DO
250# endif
251         END DO
252         ! Lateral boundary conditions on each gradient
253         CALL lbc_lnk( pgru_tl , 'U', -1.0_wp )  ;  CALL lbc_lnk( pgrv_tl , 'V', -1.0_wp )
254      END IF
255      !
256      CALL wrk_dealloc( jpi, jpj,       zri, zrj, zhi, zhj, zritl, zrjtl )
257      CALL wrk_dealloc( jpi, jpj, kjpt, zti, ztj, ztitl, ztjtl           )
258      !
259      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde_tan')
260      !
261   END SUBROUTINE zps_hde_tan
262
263   SUBROUTINE zps_hde_adj ( kt, kjpt, pta, pgtu, pgtv,       &
264      &                      pta_ad,  prd_ad,    &
265      &                      pgtu_ad, pgru_ad,   &
266      &                      pgtv_ad, pgrv_ad     )
267      !!----------------------------------------------------------------------
268      !!                     ***  ROUTINE zps_hde_adj  ***
269      !!
270      !! ** Purpose of the direct routine:
271      !!      Compute the horizontal derivative of T, S and rd
272      !!      at u- and v-points with a linear interpolation for z-coordinate
273      !!      with partial steps.
274      !!
275      !! ** Method of the direct routine:
276      !!      In z-coord with partial steps, scale factors on last
277      !!      levels are different for each grid point, so that T, S and rd
278      !!      points are not at the same depth as in z-coord. To have horizontal
279      !!      gradients again, we interpolate T and S at the good depth :
280      !!      Linear interpolation of T, S
281      !!         Computation of di(tb) and dj(tb) by vertical interpolation:
282      !!          di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~
283      !!          dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~
284      !!         This formulation computes the two cases:
285      !!                 CASE 1                   CASE 2
286      !!         k-1  ___ ___________   k-1   ___ ___________
287      !!                    Ti  T~                  T~  Ti+1
288      !!                  _____                        _____
289      !!         k        |   |Ti+1     k           Ti |   |
290      !!                  |   |____                ____|   |
291      !!              ___ |   |   |           ___  |   |   |
292      !!
293      !!      case 1->   e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
294      !!          t~ = t(i+1,j  ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
295      !!        ( t~ = t(i  ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1)  )
296      !!          or
297      !!      case 2->   e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
298      !!          t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
299      !!        ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
300      !!          Idem for di(s) and dj(s)
301      !!
302      !!      For rho, we call eos_insitu_2d which will compute rd~(t~,s~) at
303      !!      the good depth zh from interpolated T and S for the different
304      !!      formulation of the equation of state (eos).
305      !!      Gradient formulation for rho :
306      !!          di(rho) = rd~ - rd(i,j,k) or rd (i+1,j,k) - rd~
307      !!
308      !! ** Action  : - pgtu, pgsu, pgru: horizontal gradient of T, S
309      !!                and rd at U-points
310      !!              - pgtv, pgsv, pgrv: horizontal gradient of T, S
311      !!                and rd at V-points
312      !!
313      !! History of the direct routine:
314      !!   8.5  !  02-04  (A. Bozec)  Original code
315      !!   8.5  !  02-08  (G. Madec E. Durand)  Optimization and Free form
316      !! History of the TAM routine:
317      !!   9.0  !  08-06 (A. Vidard) Skeleton
318      !!        !  08-08 (A. Vidard) adjoint of the 02-08 version
319      !!----------------------------------------------------------------------
320      !! * Arguments
321      INTEGER, INTENT( in ) ::   kt, kjpt          ! ocean time-step index
322      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT( inout ) ::   &
323         pta, pta_ad            ! 3D T, S and rd direct fields
324      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ), OPTIONAL ::   &
325         prd_ad                                              ! 3D T, S and rd tangent fields
326      REAL(wp), DIMENSION(jpi,jpj,kjpt), INTENT( inout ), OPTIONAL ::   &
327         pgtu, pgtv, pgtu_ad, pgtv_ad                                     ! 3D T, S and rd tangent fields
328      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ), OPTIONAL ::   &
329         pgru_ad,                               &  ! horizontal grad. of T, S and rd at u-
330         pgrv_ad                                   ! and v-points of the partial step level
331      !! * Local declarations
332      INTEGER ::   ji, jj,jk, jn,  &                    ! Dummy loop indices
333                &   iku,ikv, ikum1, ikvm1          ! partial step level at u- and v-points
334      REAL(wp), POINTER, DIMENSION(:,:) ::   &
335         zri, zrj,               &  ! and rd
336         zriad, zrjad,           &  ! and rdtl
337         zhi, zhj                 ! depth of interpolation for eos2d
338      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zti, ztj, ztiad, ztjad    ! interpolated value of tracer
339      REAL(wp) ::   &
340         ze3wu, ze3wv,           &  ! temporary scalars
341         zmaxu, zmaxu2,         &  !    "         "
342         zmaxv, zmaxv2             !    "         "
343      !!---------------------------------------------------------------------
344      !
345      IF( nn_timing == 1 )  CALL timing_start( 'zps_hde_adj')
346      !
347      CALL wrk_alloc( jpi, jpj,       zri, zrj, zhi, zhj, zriad, zrjad )
348      CALL wrk_alloc( jpi, jpj, kjpt, zti, ztj, ztiad, ztjad           )
349      !
350      ! 1. Direct model recomputation
351      DO jn = 1, kjpt      !==   Interpolation of tracers at the last ocean level   ==!
352         !
353# if defined key_vectopt_loop
354         jj = 1
355         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
356# else
357         DO jj = 1, jpjm1
358            DO ji = 1, jpim1
359# endif
360               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
361               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
362               ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
363               ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
364               !
365               ! i- direction
366               IF( ze3wu >= 0._wp ) THEN      ! case 1
367                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
368                  ! interpolated values of tracers
369                  zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
370                  ! gradient of  tracers
371                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
372               ELSE                           ! case 2
373                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
374                  ! interpolated values of tracers
375                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
376                  ! gradient of tracers
377                  pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
378               ENDIF
379               !
380               ! j- direction
381               IF( ze3wv >= 0._wp ) THEN      ! case 1
382                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
383                  ! interpolated values of tracers
384                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
385                  ! gradient of tracers
386                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
387               ELSE                           ! case 2
388                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
389                  ! interpolated values of tracers
390                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
391                  ! gradient of tracers
392                  pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
393               ENDIF
394# if ! defined key_vectopt_loop
395            END DO
396# endif
397         END DO
398         CALL lbc_lnk( pgtu(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk( pgtv(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
399         !
400      END DO
401
402      !2. Adjoint model counterpart
403      ztiad = 0.0_wp ; ztjad = 0.0_wp
404      zriad = 0.0_wp ; zrjad = 0.0_wp
405! horizontal derivative of density anomalies (rd)
406      IF( PRESENT( prd_ad ) ) THEN         ! depth of the partial step level
407         ! Lateral boundary conditions on each gradient
408         CALL lbc_lnk_adj( pgru_ad , 'U', -1.0_wp )
409         CALL lbc_lnk_adj( pgrv_ad , 'V', -1.0_wp )
410# if defined key_vectopt_loop
411         jj = 1
412         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
413# else
414         DO jj = 1, jpjm1
415            DO ji = 1, jpim1
416# endif
417               iku = mbku(ji,jj)
418               ikv = mbkv(ji,jj)
419               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
420               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
421               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1
422               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2
423               ENDIF
424               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1
425               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2
426               ENDIF
427# if ! defined key_vectopt_loop
428            END DO
429# endif
430         END DO
431         ! Gradient of density at the last level
432# if defined key_vectopt_loop
433            jj = 1
434            DO ji = jpij-jpi, -1    ! vector opt. (forced unrolled)
435# else
436         DO jj = jpjm1, 1, -1
437            DO ji = jpim1, 1, -1
438# endif
439               iku = mbku(ji,jj)
440               ikv = mbkv(ji,jj)
441               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
442               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
443               IF( ze3wv >= 0. ) THEN    ! j-direction: case 1
444                  zrjad(ji,jj)        = zrjad(ji,jj)                     &
445                     &                + pgrv_ad(ji,jj) * vmask(ji,jj,1)
446                  prd_ad(ji,jj,ikv)   = prd_ad(ji,jj,ikv)                &
447                     &                - pgrv_ad(ji,jj) * vmask(ji,jj,1)
448                  pgrv_ad(ji,jj)      = 0.0_wp
449               ELSE                      ! j-direction: case 2
450                  prd_ad(ji,jj+1,ikv) = prd_ad(ji,jj+1,ikv)            &
451                     &                + pgrv_ad(ji,jj) * vmask(ji,jj,1)
452                  zrjad(ji,jj)        = zrjad(ji,jj)                     &
453                     &                - pgrv_ad(ji,jj) * vmask(ji,jj,1)
454                  pgrv_ad(ji,jj)      = 0.0_wp
455               ENDIF
456               IF( ze3wu >= 0. ) THEN    ! i-direction: case 1
457                  zriad(ji,jj)        = zriad(ji,jj)                     &
458                     &                + pgru_ad(ji,jj) * umask(ji,jj,1)
459                  prd_ad(ji,jj,iku)   = prd_ad(ji,jj,iku)                &
460                     &                - pgru_ad(ji,jj) * umask(ji,jj,1)
461                  pgru_ad(ji,jj)      = 0.0_wp
462               ELSE                      ! i-direction: case 2
463                  prd_ad(ji+1,jj,iku) = prd_ad(ji+1,jj,iku)            &
464                     &                + pgru_ad(ji,jj) * umask(ji,jj,1)
465                  zriad(ji,jj)        = zriad(ji,jj)                   &
466                     &                - pgru_ad(ji,jj) * umask(ji,jj,1)
467                  pgru_ad(ji,jj)      = 0.0_wp
468               ENDIF
469
470# if ! defined key_vectopt_loop
471            END DO
472# endif
473         END DO
474
475         ! Compute interpolated rd from zti, zsi, ztj, zsj for the 2 cases at the depth of the partial
476         ! step and store it in  zri, zrj for each  case
477         CALL eos_adj( ztj, zhj, ztjad, zrjad )
478         CALL eos_adj( zti, zhi, ztiad, zriad )
479      END IF
480
481      DO jn = 1, kjpt
482         CALL lbc_lnk_adj( pgtu_ad(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk_adj( pgtv_ad(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
483# if defined key_vectopt_loop
484            jj = 1
485            DO ji = jpij-jpi, 1, -1   ! vector opt. (forced unrolled)
486# else
487         DO jj = jpjm1, 1, -1
488            DO ji = jpim1, 1, -1
489# endif
490               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
491               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
492               ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
493               ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
494               !
495               ! j- direction
496               IF( ze3wv >= 0. ) THEN      ! case 1
497                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
498                  ! gradient of T and S
499                  ztjad(ji,jj,jn)           = ztjad(ji,jj,jn)                &
500                     &                   + pgtv_ad(ji,jj,jn) * vmask(ji,jj,1)
501                  pta_ad(ji,jj,ikv,jn)     = pta_ad(ji,jj,ikv,jn)          &
502                     &                   - pgtv_ad(ji,jj,jn) * vmask(ji,jj,1)
503                  pgtv_ad(ji,jj,jn)         = 0.0_wp
504                  ! interpolated values of T and S
505                  pta_ad(ji,jj+1,ikv,jn)   = pta_ad(ji,jj+1,ikv,jn)        &
506                     &                   + ztjad(ji,jj,jn) * (1 - zmaxv)
507                  pta_ad(ji,jj+1,ikv-1,jn) = pta_ad(ji,jj+1,ikv-1,jn)      &
508                     &                   + ztjad(ji,jj,jn)* zmaxv
509                  ztjad(ji,jj,jn)           = 0.0_wp
510               ELSE                        ! case 2
511                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
512                  ! gradient of T and S
513                  pta_ad(ji,jj+1,ikv,jn)   = pta_ad(ji,jj+1,ikv,jn)        &
514                     &                   + pgtv_ad(ji,jj,jn) * vmask(ji,jj,1)
515                  ztjad(ji,jj,jn)           = ztjad(ji,jj,jn)                &
516                     &                   - pgtv_ad(ji,jj,jn) * vmask(ji,jj,1)
517                  pgtv_ad(ji,jj,jn)         = 0.0_wp
518
519                  ! interpolated values of T and S
520                  pta_ad(ji,jj,ikv,jn)     = pta_ad(ji,jj,ikv,jn)          &
521                     &                   + ztjad(ji,jj,jn) * (1 - zmaxv)
522                  pta_ad(ji,jj,ikv-1,jn)   = pta_ad(ji,jj,ikv-1,jn)        &
523                     &                   + ztjad(ji,jj,jn) * zmaxv
524                  ztjad(ji,jj,jn)           = 0.0_wp
525               ENDIF
526               ! i- direction
527               IF( ze3wu >= 0. ) THEN      ! case 1
528                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
529                  ! gradient of T and S
530                  ztiad(ji,jj,jn)       = ztiad(ji,jj,jn)                    &
531                     &               + pgtu_ad(ji,jj,jn) * umask(ji,jj,1)
532                  pta_ad(ji,jj,iku,jn) = pta_ad(ji,jj,iku,jn)              &
533                     &               - pgtu_ad(ji,jj,jn) * umask(ji,jj,1)
534                  pgtu_ad(ji,jj,jn)     = 0.0_wp
535                  ! interpolated values of T and S
536                  pta_ad(ji+1,jj,iku,jn)   = pta_ad(ji+1,jj,iku,jn)        &
537                     &                   + ztiad(ji,jj,jn) * (1 - zmaxu)
538                  pta_ad(ji+1,jj,iku-1,jn) = pta_ad(ji+1,jj,iku-1,jn)      &
539                     &                   + ztiad(ji,jj,jn) * zmaxu
540                  ztiad(ji,jj,jn)           = 0.0_wp
541               ELSE                        ! case 2
542                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
543                  ! gradient of T and S
544                  pta_ad(ji+1,jj,iku,jn)   = pta_ad(ji+1,jj,iku,jn)        &
545                     &                   + pgtu_ad(ji,jj,jn) * umask(ji,jj,1)
546                  ztiad (ji,jj,jn)          = ztiad (ji,jj,jn)               &
547                     &                   - pgtu_ad(ji,jj,jn) * umask(ji,jj,1)
548                  pgtu_ad(ji,jj,jn)         = 0.0_wp
549                  ! interpolated values of T and S
550                  pta_ad(ji,jj,iku,jn)     = pta_ad(ji,jj,iku,jn)          &
551                     &                   + ztiad(ji,jj,jn) * (1 - zmaxu)
552                  pta_ad(ji,jj,iku-1,jn)   = pta_ad(ji,jj,iku-1,jn)        &
553                     &                   + ztiad(ji,jj,jn) * zmaxu
554                  ztiad(ji,jj,jn)           = 0.0_wp
555               ENDIF
556# if ! defined key_vectopt_loop
557            END DO
558# endif
559         END DO
560      END DO
561      !
562      CALL wrk_dealloc( jpi, jpj,       zri, zrj, zhi, zhj, zriad, zrjad )
563      CALL wrk_dealloc( jpi, jpj, kjpt, zti, ztj, ztiad, ztjad           )
564      !
565      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde_adj')
566      !
567   END SUBROUTINE zps_hde_adj
568   !SUBROUTINE zps_hde_adj_tst( kumadt )
569      !!!-----------------------------------------------------------------------
570      !!!
571      !!!                  ***  ROUTINE zps_hde_adj_tst ***
572      !!!
573      !!! ** Purpose : Test the adjoint routine.
574      !!!
575      !!! ** Method  : Verify the scalar product
576      !!!
577      !!!                 ( L dx )^T W dy  =  dx^T L^T W dy
578      !!!
579      !!!              where  L   = tangent routine
580      !!!                     L^T = adjoint routine
581      !!!                     W   = diagonal matrix of scale factors
582      !!!                     dx  = input perturbation (random field)
583      !!!                     dy  = L dx
584      !!!
585      !!!
586      !!! History :
587      !!!        ! 08-08 (A. Vidard)
588      !!!-----------------------------------------------------------------------
589      !!! * Modules used
590
591      !!! * Arguments
592      !INTEGER, INTENT(IN) :: &
593         !& kumadt             ! Output unit
594
595      !INTEGER :: &
596         !& ji,    &        ! dummy loop indices
597         !& jj,    &
598         !& jk,    &
599         !& kt,    &
600         !& jn
601
602      !!! * Local declarations
603      !REAL(KIND=wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
604         !& zts,         & ! Direct field : temperature
605         !& zts_tlin,    & ! Tangent input: temperature
606         !& zts_adout,   & ! Adjoint output: temperature
607         !& zats           ! 3D random field for t
608
609      !REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
610         !& zgtu_tlout,   & ! Tangent output: horizontal gradient
611         !& zgtv_tlout,   & ! Tangent output: horizontal gradient
612         !& zrd_adout,    & ! Adjoint output:
613         !& zar,          & ! 3D random field for rd
614         !& zrd_tlin,     & ! Tangent input:
615         !& zgtu_adin,    & ! Adjoint input : horizontal gradient
616         !& zgtv_adin       ! Adjoint input : horizontal gradient
617
618      !REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
619         !& zgru_tlout,   & ! Tangent output: horizontal gradient
620         !& zgrv_tlout,   & ! Tangent output: horizontal gradient
621         !& zgru_adin,    & ! Adjoint input : horizontal gradient
622         !& zgrv_adin       ! Adjoint input : horizontal gradient
623
624      !REAL(KIND=wp) :: &
625                           !! random field standard deviation for:
626         !& zsp1,         & ! scalar product involving the tangent routine
627         !& zsp1_1,       & !   scalar product components
628         !& zsp1_2,       &
629         !& zsp1_3,       & !   scalar product components
630         !& zsp1_4,       &
631         !& zsp1_5,       & !   scalar product components
632         !& zsp1_6,       &
633         !& zsp2,         & ! scalar product involving the adjoint routine
634         !& zsp2_1,       & !   scalar product components
635         !& zsp2_2,       &
636         !& zsp2_3
637      !CHARACTER (LEN=14) :: &
638         !& cl_name
639
640      !kt = nit000
641      !! Allocate memory
642      !ALLOCATE( &
643         !& zts(jpi,jpj,jpk,jpts),         &
644         !& zts_tlin(jpi,jpj,jpk,jpts),    &
645         !& zrd_tlin(jpi,jpj,jpk),     &
646         !& zts_adout(jpi,jpj,jpk,jpts),   &
647         !& zrd_adout(jpi,jpj,jpk),    &
648         !& zar(jpi,jpj,jpk),          &
649         !& zats(jpi,jpj,jpk,jpts),          &
650         !& zgtu_tlout(jpi,jpj,jpts),       &
651         !& zgtv_tlout(jpi,jpj,jpts),       &
652         !& zgru_tlout(jpi,jpj),       &
653         !& zgrv_tlout(jpi,jpj),       &
654         !& zgtu_adin(jpi,jpj,jpts),        &
655         !& zgtv_adin(jpi,jpj,jpts),        &
656         !& zgru_adin(jpi,jpj),        &
657         !& zgrv_adin(jpi,jpj)         &
658         !& )
659      !! Initialize random field standard deviationsthe reference state
660      !zts = tsn(:,:,:,:)
661
662      !!=============================================================
663      !! 1) dx = ( T ) and dy = ( T )
664      !!=============================================================
665
666      !!--------------------------------------------------------------------
667      !! Reset the tangent and adjoint variables
668      !!--------------------------------------------------------------------
669      !zts_tlin(:,:,:,:)  = 0.0_wp
670      !zrd_tlin(:,:,:)   = 0.0_wp
671      !zts_adout(:,:,:,:) = 0.0_wp
672      !zrd_adout(:,:,:)  = 0.0_wp
673      !zgtu_tlout(:,:,:)   = 0.0_wp
674      !zgtv_tlout(:,:,:)   = 0.0_wp
675      !zgru_tlout(:,:)   = 0.0_wp
676      !zgrv_tlout(:,:)   = 0.0_wp
677      !zgtu_adin(:,:,:)    = 0.0_wp
678      !zgtv_adin(:,:,:)    = 0.0_wp
679      !zgru_adin(:,:)    = 0.0_wp
680      !zgrv_adin(:,:)    = 0.0_wp
681
682      !!--------------------------------------------------------------------
683      !! Initialize the tangent input with random noise: dx
684      !!--------------------------------------------------------------------
685      !DO jn = 1, jpts
686         !CALL grid_random(  zats(:,:,:,jn), 'T', 0.0_wp, stdt )
687      !END DO
688      !CALL grid_random(  zar, 'T', 0.0_wp, stdr )
689
690      !zts_tlin(:,:,:,:) = zats(:,:,:,:)
691      !zrd_tlin(:,:,:) = zar(:,:,:)
692      !CALL zps_hde_tan ( nit000, jpts, zts,           &
693         !&                   zts_tlin , zrd_tlin  ,   &
694         !&                   zgtu_tlout, zgru_tlout,  &
695         !&                   zgtv_tlout, zgrv_tlout   )
696      !DO jn = 1, jpts
697         !DO jj = nldj, nlej
698            !DO ji = nldi, nlei
699               !jk = mbku(ji,jj)
700               !zgtu_adin(ji,jj,jn) = zgtu_tlout(ji,jj,jn) &
701                  !&             * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
702                  !&             * umask(ji,jj,jk)
703               !jk = mbkv(ji,jj)
704               !zgtv_adin(ji,jj,jn) = zgtv_tlout(ji,jj,jn) &
705                  !&             * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
706                  !&             * vmask(ji,jj,jk)
707            !END DO
708         !END DO
709      !END DO
710      !DO jj = nldj, nlej
711         !DO ji = nldi, nlei
712            !jk = mbku(ji,jj)
713            !zgru_adin(ji,jj) = zgru_tlout(ji,jj) &
714               !&             * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
715               !&             * umask(ji,jj,jk)
716            !jk = mbkv(ji,jj)
717            !zgrv_adin(ji,jj) = zgrv_tlout(ji,jj) &
718               !&             * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
719               !&             * vmask(ji,jj,jk)
720         !END DO
721      !END DO
722      !!--------------------------------------------------------------------
723      !! Compute the scalar product: ( L dx )^T W dy
724      !!--------------------------------------------------------------------
725
726      !zsp1_1 = DOT_PRODUCT( zgtu_adin(:,:,jp_tem), zgtu_tlout(:,:,jp_tem) )
727      !zsp1_2 = DOT_PRODUCT( zgtu_adin(:,:,jp_sal), zgtu_tlout(:,:,jp_sal) )
728      !zsp1_3 = DOT_PRODUCT( zgru_adin, zgru_tlout )
729      !zsp1_4 = DOT_PRODUCT( zgtv_adin(:,:,jp_tem), zgtv_tlout(:,:,jp_tem) )
730      !zsp1_5 = DOT_PRODUCT( zgtv_adin(:,:,jp_sal), zgtv_tlout(:,:,jp_sal) )
731      !zsp1_6 = DOT_PRODUCT( zgrv_adin, zgrv_tlout )
732      !zsp1 = zsp1_1 + zsp1_2 + zsp1_3 + zsp1_4 + zsp1_5 + zsp1_6
733
734
735      !!--------------------------------------------------------------------
736      !! Call the adjoint routine: dx^* = L^T dy^*
737      !!--------------------------------------------------------------------
738      !CALL zps_hde_adj ( kt, jpts, zts, gtsu, gtsv ,   &
739         !&                   tsa_ad,   ,   &
740         !&                   zgtu_adin , zgru_adin ,   &
741         !&                   zgtv_adin , zgrv_adin     )
742
743      !zsp2_1 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem) )
744      !zsp2_2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal) )
745      !zsp2_3 = DOT_PRODUCT( zrd_tlin , zrd_adout  )
746      !zsp2   = zsp2_1 + zsp2_2 + zsp2_3
747
748      !! Compare the scalar products
749
750      !cl_name = 'zps_hde_adj   '
751      !CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
752
753      !! Deallocate memory
754      !DEALLOCATE( &
755         !& zts,         &
756         !& zts_tlin,    &
757         !& zrd_tlin,     &
758         !& zts_adout,   &
759         !& zrd_adout,    &
760         !& zar,          &
761         !& zats,          &
762         !& zgtu_tlout,   &
763         !& zgtv_tlout,   &
764         !& zgru_tlout,   &
765         !& zgrv_tlout,   &
766         !& zgtu_adin,    &
767         !& zgtv_adin,    &
768         !& zgru_adin,    &
769         !& zgrv_adin     &
770         !& )
771
772   !END SUBROUTINE zps_hde_adj_tst
773#endif
774END MODULE zpshde_tam
Note: See TracBrowser for help on using the repository browser.