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_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/TRA – NEMO

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/TRA/zpshde_tam.F90 @ 4594

Last change on this file since 4594 was 4594, checked in by pabouttier, 10 years ago

Remove pgtu and pgtv args in zps_hde_adj definition, see Ticket #1283
  • Property svn:executable set to *
File size: 35.0 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( in ), 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,       &
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_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               ELSE                           ! case 2
371                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
372                  ! interpolated values of tracers
373                  zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
374               ENDIF
375               !
376               ! j- direction
377               IF( ze3wv >= 0._wp ) THEN      ! case 1
378                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
379                  ! interpolated values of tracers
380                  ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
381               ELSE                           ! case 2
382                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
383                  ! interpolated values of tracers
384                  ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
385               ENDIF
386# if ! defined key_vectopt_loop
387            END DO
388# endif
389         END DO
390         !
391      END DO
392
393      !2. Adjoint model counterpart
394      ztiad = 0.0_wp ; ztjad = 0.0_wp
395      zriad = 0.0_wp ; zrjad = 0.0_wp
396! horizontal derivative of density anomalies (rd)
397      IF( PRESENT( prd_ad ) ) THEN         ! depth of the partial step level
398         ! Lateral boundary conditions on each gradient
399         CALL lbc_lnk_adj( pgru_ad , 'U', -1.0_wp )
400         CALL lbc_lnk_adj( pgrv_ad , 'V', -1.0_wp )
401# if defined key_vectopt_loop
402         jj = 1
403         DO ji = 1, jpij-jpi   ! vector opt. (forced unrolled)
404# else
405         DO jj = 1, jpjm1
406            DO ji = 1, jpim1
407# endif
408               iku = mbku(ji,jj)
409               ikv = mbkv(ji,jj)
410               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
411               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
412               IF( ze3wu >= 0._wp ) THEN   ;   zhi(ji,jj) = fsdept(ji  ,jj,iku)     ! i-direction: case 1
413               ELSE                        ;   zhi(ji,jj) = fsdept(ji+1,jj,iku)     ! -     -      case 2
414               ENDIF
415               IF( ze3wv >= 0._wp ) THEN   ;   zhj(ji,jj) = fsdept(ji,jj  ,ikv)     ! j-direction: case 1
416               ELSE                        ;   zhj(ji,jj) = fsdept(ji,jj+1,ikv)     ! -     -      case 2
417               ENDIF
418# if ! defined key_vectopt_loop
419            END DO
420# endif
421         END DO
422         ! Gradient of density at the last level
423# if defined key_vectopt_loop
424            jj = 1
425            DO ji = jpij-jpi, -1    ! vector opt. (forced unrolled)
426# else
427         DO jj = jpjm1, 1, -1
428            DO ji = jpim1, 1, -1
429# endif
430               iku = mbku(ji,jj)
431               ikv = mbkv(ji,jj)
432               ze3wu  = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
433               ze3wv  = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
434               IF( ze3wv >= 0. ) THEN    ! j-direction: case 1
435                  zrjad(ji,jj)        = zrjad(ji,jj)                     &
436                     &                + pgrv_ad(ji,jj) * vmask(ji,jj,1)
437                  prd_ad(ji,jj,ikv)   = prd_ad(ji,jj,ikv)                &
438                     &                - pgrv_ad(ji,jj) * vmask(ji,jj,1)
439                  pgrv_ad(ji,jj)      = 0.0_wp
440               ELSE                      ! j-direction: case 2
441                  prd_ad(ji,jj+1,ikv) = prd_ad(ji,jj+1,ikv)            &
442                     &                + pgrv_ad(ji,jj) * vmask(ji,jj,1)
443                  zrjad(ji,jj)        = zrjad(ji,jj)                     &
444                     &                - pgrv_ad(ji,jj) * vmask(ji,jj,1)
445                  pgrv_ad(ji,jj)      = 0.0_wp
446               ENDIF
447               IF( ze3wu >= 0. ) THEN    ! i-direction: case 1
448                  zriad(ji,jj)        = zriad(ji,jj)                     &
449                     &                + pgru_ad(ji,jj) * umask(ji,jj,1)
450                  prd_ad(ji,jj,iku)   = prd_ad(ji,jj,iku)                &
451                     &                - pgru_ad(ji,jj) * umask(ji,jj,1)
452                  pgru_ad(ji,jj)      = 0.0_wp
453               ELSE                      ! i-direction: case 2
454                  prd_ad(ji+1,jj,iku) = prd_ad(ji+1,jj,iku)            &
455                     &                + pgru_ad(ji,jj) * umask(ji,jj,1)
456                  zriad(ji,jj)        = zriad(ji,jj)                   &
457                     &                - pgru_ad(ji,jj) * umask(ji,jj,1)
458                  pgru_ad(ji,jj)      = 0.0_wp
459               ENDIF
460
461# if ! defined key_vectopt_loop
462            END DO
463# endif
464         END DO
465
466         ! Compute interpolated rd from zti, zsi, ztj, zsj for the 2 cases at the depth of the partial
467         ! step and store it in  zri, zrj for each  case
468         CALL eos_adj( ztj, zhj, ztjad, zrjad )
469         CALL eos_adj( zti, zhi, ztiad, zriad )
470      END IF
471
472      DO jn = 1, kjpt
473         CALL lbc_lnk_adj( pgtu_ad(:,:,jn), 'U', -1. )   ;   CALL lbc_lnk_adj( pgtv_ad(:,:,jn), 'V', -1. )   ! Lateral boundary cond.
474# if defined key_vectopt_loop
475            jj = 1
476            DO ji = jpij-jpi, 1, -1   ! vector opt. (forced unrolled)
477# else
478         DO jj = jpjm1, 1, -1
479            DO ji = jpim1, 1, -1
480# endif
481               iku = mbku(ji,jj)   ;   ikum1 = MAX( iku - 1 , 1 )    ! last and before last ocean level at u- & v-points
482               ikv = mbkv(ji,jj)   ;   ikvm1 = MAX( ikv - 1 , 1 )    ! if level first is a p-step, ik.m1=1
483               ze3wu = fse3w(ji+1,jj  ,iku) - fse3w(ji,jj,iku)
484               ze3wv = fse3w(ji  ,jj+1,ikv) - fse3w(ji,jj,ikv)
485               !
486               ! j- direction
487               IF( ze3wv >= 0. ) THEN      ! case 1
488                  zmaxv =  ze3wv / fse3w(ji,jj+1,ikv)
489                  ! gradient of T and S
490                  ztjad(ji,jj,jn)           = ztjad(ji,jj,jn)                &
491                     &                   + pgtv_ad(ji,jj,jn) * vmask(ji,jj,1)
492                  pta_ad(ji,jj,ikv,jn)     = pta_ad(ji,jj,ikv,jn)          &
493                     &                   - pgtv_ad(ji,jj,jn) * vmask(ji,jj,1)
494                  pgtv_ad(ji,jj,jn)         = 0.0_wp
495                  ! interpolated values of T and S
496                  pta_ad(ji,jj+1,ikv,jn)   = pta_ad(ji,jj+1,ikv,jn)        &
497                     &                   + ztjad(ji,jj,jn) * (1 - zmaxv)
498                  pta_ad(ji,jj+1,ikvm1,jn) = pta_ad(ji,jj+1,ikvm1,jn)      &
499                     &                   + ztjad(ji,jj,jn)* zmaxv
500                  ztjad(ji,jj,jn)           = 0.0_wp
501               ELSE                        ! case 2
502                  zmaxv =  -ze3wv / fse3w(ji,jj,ikv)
503                  ! gradient of T and S
504                  pta_ad(ji,jj+1,ikv,jn)   = pta_ad(ji,jj+1,ikv,jn)        &
505                     &                   + pgtv_ad(ji,jj,jn) * vmask(ji,jj,1)
506                  ztjad(ji,jj,jn)           = ztjad(ji,jj,jn)                &
507                     &                   - pgtv_ad(ji,jj,jn) * vmask(ji,jj,1)
508                  pgtv_ad(ji,jj,jn)         = 0.0_wp
509
510                  ! interpolated values of T and S
511                  pta_ad(ji,jj,ikv,jn)     = pta_ad(ji,jj,ikv,jn)          &
512                     &                   + ztjad(ji,jj,jn) * (1 - zmaxv)
513                  pta_ad(ji,jj,ikvm1,jn)   = pta_ad(ji,jj,ikvm1,jn)        &
514                     &                   + ztjad(ji,jj,jn) * zmaxv
515                  ztjad(ji,jj,jn)           = 0.0_wp
516               ENDIF
517               ! i- direction
518               IF( ze3wu >= 0. ) THEN      ! case 1
519                  zmaxu =  ze3wu / fse3w(ji+1,jj,iku)
520                  ! gradient of T and S
521                  ztiad(ji,jj,jn)       = ztiad(ji,jj,jn)                    &
522                     &               + pgtu_ad(ji,jj,jn) * umask(ji,jj,1)
523                  pta_ad(ji,jj,iku,jn) = pta_ad(ji,jj,iku,jn)              &
524                     &               - pgtu_ad(ji,jj,jn) * umask(ji,jj,1)
525                  pgtu_ad(ji,jj,jn)     = 0.0_wp
526                  ! interpolated values of T and S
527                  pta_ad(ji+1,jj,iku,jn)   = pta_ad(ji+1,jj,iku,jn)        &
528                     &                   + ztiad(ji,jj,jn) * (1 - zmaxu)
529                  pta_ad(ji+1,jj,ikum1,jn) = pta_ad(ji+1,jj,ikum1,jn)      &
530                     &                   + ztiad(ji,jj,jn) * zmaxu
531                  ztiad(ji,jj,jn)           = 0.0_wp
532               ELSE                        ! case 2
533                  zmaxu = -ze3wu / fse3w(ji,jj,iku)
534                  ! gradient of T and S
535                  pta_ad(ji+1,jj,iku,jn)   = pta_ad(ji+1,jj,iku,jn)        &
536                     &                   + pgtu_ad(ji,jj,jn) * umask(ji,jj,1)
537                  ztiad (ji,jj,jn)          = ztiad (ji,jj,jn)               &
538                     &                   - pgtu_ad(ji,jj,jn) * umask(ji,jj,1)
539                  pgtu_ad(ji,jj,jn)         = 0.0_wp
540                  ! interpolated values of T and S
541                  pta_ad(ji,jj,iku,jn)     = pta_ad(ji,jj,iku,jn)          &
542                     &                   + ztiad(ji,jj,jn) * (1 - zmaxu)
543                  pta_ad(ji,jj,ikum1,jn)   = pta_ad(ji,jj,ikum1,jn)        &
544                     &                   + ztiad(ji,jj,jn) * zmaxu
545                  ztiad(ji,jj,jn)           = 0.0_wp
546               ENDIF
547# if ! defined key_vectopt_loop
548            END DO
549# endif
550         END DO
551      END DO
552      !
553      CALL wrk_dealloc( jpi, jpj,       zri, zrj, zhi, zhj, zriad, zrjad )
554      CALL wrk_dealloc( jpi, jpj, kjpt, zti, ztj, ztiad, ztjad           )
555      !
556      IF( nn_timing == 1 )  CALL timing_stop( 'zps_hde_adj')
557      !
558   END SUBROUTINE zps_hde_adj
559      SUBROUTINE zps_hde_adj_tst( kumadt )
560      !!-----------------------------------------------------------------------
561      !!
562      !!                  ***  ROUTINE zps_hde_adj_tst ***
563      !!
564      !! ** Purpose : Test the adjoint routine.
565      !!
566      !! ** Method  : Verify the scalar product
567      !!
568      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
569      !!
570      !!              where  L   = tangent routine
571      !!                     L^T = adjoint routine
572      !!                     W   = diagonal matrix of scale factors
573      !!                     dx  = input perturbation (random field)
574      !!                     dy  = L dx
575      !!
576      !!
577      !! History :
578      !!        ! 08-08 (A. Vidard)
579      !!-----------------------------------------------------------------------
580      !! * Modules used
581      !! * Arguments
582      INTEGER, INTENT(IN) :: &
583         & kumadt             ! Output unit
584
585      INTEGER :: &
586         & ji,    &        ! dummy loop indices
587         & jj,    &
588         & jk,    &
589         & kt,    &
590         & jn
591
592      !! * Local declarations
593      REAL(KIND=wp), DIMENSION(:,:,:,:), POINTER :: &
594         & zts,         & ! Direct field : temperature
595         & zts_tlin,    & ! Tangent input: temperature
596         & zts_adout,   & ! Adjoint output: temperature
597         & zats           ! 3D random field for t
598
599      REAL(KIND=wp), DIMENSION(:,:,:), POINTER :: &
600         & zgtu_tlout,   & ! Tangent output: horizontal gradient
601         & zgtv_tlout,   & ! Tangent output: horizontal gradient
602         & zrd_adout,    & ! Adjoint output:
603         & zar,          & ! 3D random field for rd
604         & zrd_tlin,     & ! Tangent input:
605         & zgtu_adin,    & ! Adjoint input : horizontal gradient
606         & zgtv_adin       ! Adjoint input : horizontal gradient
607
608      REAL(KIND=wp), DIMENSION(:,:), POINTER :: &
609         & zgru_tlout,   & ! Tangent output: horizontal gradient
610         & zgrv_tlout,   & ! Tangent output: horizontal gradient
611         & zgru_adin,    & ! Adjoint input : horizontal gradient
612         & zgrv_adin       ! Adjoint input : horizontal gradient
613
614      REAL(KIND=wp) :: &
615                           ! random field standard deviation for:
616         & zsp1,         & ! scalar product involving the tangent routine
617         & zsp1_1,       & !   scalar product components
618         & zsp1_2,       &
619         & zsp1_3,       & !   scalar product components
620         & zsp1_4,       &
621         & zsp1_5,       & !   scalar product components
622         & zsp1_6,       &
623         & zsp2,         & ! scalar product involving the adjoint routine
624         & zsp2_1,       & !   scalar product components
625         & zsp2_2,       &
626         & zsp2_3
627      CHARACTER (LEN=14) :: &
628         & cl_name
629
630      kt = nit000
631      ! Allocate memory
632      CALL wrk_alloc(jpi, jpj, jpk, jpts, zts, zts_tlin,zts_adout, zats )
633
634      CALL wrk_alloc(jpi, jpj, jpts, zgtu_tlout, zgtv_tlout, zgtu_adin, zgtv_adin )
635
636      CALL wrk_alloc(jpi, jpj, jpk , zrd_adout, zrd_tlin, zar )
637
638      CALL wrk_alloc(jpi, jpj, zgru_tlout, zgrv_tlout, zgru_adin, zgrv_adin )
639
640      ! Initialize random field standard deviationsthe reference state
641      zts = tsn(:,:,:,:)
642
643      !=============================================================
644      ! 1) dx = ( T ) and dy = ( T )
645      !=============================================================
646
647      !--------------------------------------------------------------------
648      ! Reset the tangent and adjoint variables
649      !--------------------------------------------------------------------
650      zts_tlin(:,:,:,:)  = 0.0_wp
651      zrd_tlin(:,:,:)   = 0.0_wp
652      zts_adout(:,:,:,:) = 0.0_wp
653      zrd_adout(:,:,:)  = 0.0_wp
654      zgtu_tlout(:,:,:)   = 0.0_wp
655      zgtv_tlout(:,:,:)   = 0.0_wp
656      zgru_tlout(:,:)   = 0.0_wp
657      zgrv_tlout(:,:)   = 0.0_wp
658      zgtu_adin(:,:,:)    = 0.0_wp
659      zgtv_adin(:,:,:)    = 0.0_wp
660      zgru_adin(:,:)    = 0.0_wp
661      zgrv_adin(:,:)    = 0.0_wp
662
663      !--------------------------------------------------------------------
664      ! Initialize the tangent input with random noise: dx
665      !--------------------------------------------------------------------
666      DO jn = 1, jpts
667         CALL grid_random(  zats(:,:,:,jn), 'T', 0.0_wp, stdt )
668      END DO
669      CALL grid_random(  zar, 'T', 0.0_wp, stdr )
670
671      DO jj = nldj, nlej
672         DO ji = nldi, nlei
673            zts_tlin(ji,jj,:,:) = zats(ji,jj,:,:)
674            zrd_tlin(ji,jj,:)   = zar (ji,jj,:)
675         END DO
676      END DO
677
678      CALL zps_hde_tan ( nit000, jpts, zts,           &
679         &                   zts_tlin , zrd_tlin  ,   &
680         &                   zgtu_tlout, zgru_tlout,  &
681         &                   zgtv_tlout, zgrv_tlout   )
682
683      DO jn = 1, jpts
684         DO jj = nldj, nlej
685            DO ji = nldi, nlei
686               jk = mbku(ji,jj)
687               zgtu_adin(ji,jj,jn) = zgtu_tlout(ji,jj,jn) &
688                  &             * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
689                  &             * umask(ji,jj,jk)
690               jk = mbkv(ji,jj)
691               zgtv_adin(ji,jj,jn) = zgtv_tlout(ji,jj,jn) &
692                  &             * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
693                  &             * vmask(ji,jj,jk)
694            END DO
695         END DO
696      END DO
697      DO jj = nldj, nlej
698         DO ji = nldi, nlei
699            jk = mbku(ji,jj)
700            zgru_adin(ji,jj) = zgru_tlout(ji,jj) &
701               &             * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
702               &             * umask(ji,jj,jk)
703            jk = mbkv(ji,jj)
704            zgrv_adin(ji,jj) = zgrv_tlout(ji,jj) &
705               &             * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
706               &             * vmask(ji,jj,jk)
707         END DO
708      END DO
709      !!--------------------------------------------------------------------
710      !! Compute the scalar product: ( L dx )^T W dy
711      !!--------------------------------------------------------------------
712
713      zsp1_1 = DOT_PRODUCT( zgtu_adin(:,:,jp_tem), zgtu_tlout(:,:,jp_tem) )
714      zsp1_2 = DOT_PRODUCT( zgtu_adin(:,:,jp_sal), zgtu_tlout(:,:,jp_sal) )
715      zsp1_3 = DOT_PRODUCT( zgru_adin, zgru_tlout )
716      zsp1_4 = DOT_PRODUCT( zgtv_adin(:,:,jp_tem), zgtv_tlout(:,:,jp_tem) )
717      zsp1_5 = DOT_PRODUCT( zgtv_adin(:,:,jp_sal), zgtv_tlout(:,:,jp_sal) )
718      zsp1_6 = DOT_PRODUCT( zgrv_adin, zgrv_tlout )
719      zsp1 = zsp1_1 + zsp1_2 + zsp1_3 + zsp1_4 + zsp1_5 + zsp1_6
720
721
722      !--------------------------------------------------------------------
723      ! Call the adjoint routine: dx^* = L^T dy^*
724      !--------------------------------------------------------------------
725      CALL zps_hde_adj ( kt, jpts, zts ,   &
726         &                   zts_adout , zrd_adout ,   &
727         &                   zgtu_adin , zgru_adin ,   &
728         &                   zgtv_adin , zgrv_adin     )
729
730      zsp2_1 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem) )
731      zsp2_2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal) )
732      zsp2_3 = DOT_PRODUCT( zrd_tlin , zrd_adout  )
733      zsp2   = zsp2_1 + zsp2_2 + zsp2_3
734
735      ! Compare the scalar products
736
737      cl_name = 'zps_hde_adj   '
738      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
739
740      ! Deallocate memory
741      CALL wrk_dealloc(jpi, jpj, jpk, jpts, zts, zts_tlin,zts_adout, zats )
742
743      CALL wrk_dealloc(jpi, jpj, jpts, zgtu_tlout, zgtv_tlout, zgtu_adin, zgtv_adin )
744
745      CALL wrk_dealloc(jpi, jpj, jpk , zrd_adout, zrd_tlin, zar )
746
747      CALL wrk_dealloc(jpi, jpj, zgru_tlout, zgrv_tlout, zgru_adin, zgrv_adin )
748     
749
750   END SUBROUTINE zps_hde_adj_tst
751#endif
752END MODULE zpshde_tam
Note: See TracBrowser for help on using the repository browser.