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.
traldf_bilap_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/traldf_bilap_tam.F90 @ 3627

Last change on this file since 3627 was 3627, checked in by rblod, 11 years ago

Correct long lines in TAM 4.3 see ticket #1010
  • Property svn:executable set to *
File size: 28.5 KB
Line 
1MODULE traldf_bilap_tam
2#ifdef key_tam
3   !!===========================================================================
4   !!                       ***  MODULE  dynldf_bilap_tam  ***
5   !! Ocean dynamics:  lateral viscosity trend
6   !!                  tangent and Adjoint Module
7   !!===========================================================================
8
9   !!---------------------------------------------------------------------------
10   !!   dyn_ldf_bilap_tan : update the momentum trend with the lateral diffusion
11   !!                       using an iso-level bilaplacian operator (tangent)
12   !!   dyn_ldf_bilap_adj : update the momentum trend with the lateral diffusion
13   !!                       using an iso-level bilaplacian operator (adjoint)
14   !!---------------------------------------------------------------------------
15
16   !! * Modules used
17   USE lbclnk
18   USE lbclnk_tam
19   USE par_oce
20   USE oce_tam
21   USE dom_oce
22   USE ldftra_oce
23   USE in_out_manager
24   USE gridrandom
25   USE dotprodfld
26   USE tstool_tam
27   USE paresp
28   USE trc_oce
29   USE lib_mpp
30   USE wrk_nemo
31   USE timing
32
33IMPLICIT NONE
34   PRIVATE
35
36   !! * Routine accessibility
37   PUBLIC tra_ldf_bilap_tan  ! called by dynldf_tam.F90
38   PUBLIC tra_ldf_bilap_adj  ! called by dynldf_tam.F90
39   PUBLIC tra_ldf_bilap_adj_tst  ! routine called by tradldf_tam.F90
40   !! * Substitutions
41#  include "domzgr_substitute.h90"
42#  include "ldftra_substitute.h90"
43#  include "ldfeiv_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46
47CONTAINS
48   SUBROUTINE tra_ldf_bilap_tan( kt, kit000, cdtype, pgu_tl, pgv_tl,      &
49      &                                  ptb_tl, pta_tl, kjpt  )
50      !!----------------------------------------------------------------------
51      !!                  ***  ROUTINE tra_ldf_bilap  ***
52      !!
53      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
54      !!      trend and add it to the general trend of tracer equation.
55      !!
56      !! ** Method  :   4th order diffusive operator along model level surfaces
57      !!      evaluated using before fields (forward time scheme). The hor.
58      !!      diffusive trends of temperature (idem for salinity) is given by:
59      !!      Laplacian of tb:
60      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(tb) ]
61      !!                                  + dj-1[ e1v*e3v/e2v dj(tb) ]  }
62      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
63      !!        zlt   = ahtt * zlt
64      !!        call to lbc_lnk
65      !!      Bilaplacian (laplacian of zlt):
66      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
67      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
68      !!      Note: if key_zco defined, e3t=e3u=e3v, they are simplified.
69      !!
70      !!      Add this trend to the general trend (ta,sa):
71      !!         (ta,sa) = (ta,sa) + ( difft , diffs )
72      !!
73      !! ** Action : - Update (ta,sa) arrays with the before iso-level
74      !!               biharmonic mixing trend.
75      !!
76      !! History :
77      !!        !  91-11  (G. Madec)  Original code
78      !!        !  93-03  (M. Guyon)  symetrical conditions
79      !!        !  95-11  (G. Madec)  suppress volumetric scale factors
80      !!        !  96-01  (G. Madec)  statement function for e3
81      !!        !  96-01  (M. Imbard)  mpp exchange
82      !!        !  97-07  (G. Madec)  optimization, and ahtt
83      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
84      !!   9.0  !  04-08  (C. talandier) New trends organization
85      !!        !  05-11  (G. Madec)  zps or sco as default option
86      !!----------------------------------------------------------------------
87      !! History of the tangent routine
88      !!   9.0  !  10-05 (P.A. Bouttier) tangent of 9.0
89      !!----------------------------------------------------------------------
90      !! * Modules used
91      !! * Arguments
92      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
93      INTEGER                              , INTENT(in   ) ::   kit000           ! first time step index
94      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype           ! =TRA or TRC (tracer indicator)
95      INTEGER                              , INTENT(in   ) ::   kjpt             ! number of tracers
96      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(in   ) ::   pgu_tl, pgv_tl   ! tracer gradient at pstep levels
97      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb_tl           ! before and now tracer fields
98      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta_tl           ! tracer trend
99      !! * Local declarations
100      INTEGER ::   ji, jj, jk, jn      ! dummy loop indices
101      INTEGER ::   iku, ikv               ! temporary integers
102      REAL(wp) ::   ztatl, zsatl, zbtr              ! temporary scalars
103      REAL(wp), POINTER, DIMENSION(:,:) ::   &
104         & zeeu, zeev,               & ! 2D workspace
105         & zlttl
106      REAL(wp), POINTER, DIMENSION(:,:,:) ::   &
107         & ztutl, ztvtl       ! 3D workspace
108      !!----------------------------------------------------------------------
109      !
110      IF( nn_timing == 1 )  CALL timing_start( 'tra_ldf_bilap_tan')
111      !
112      CALL wrk_alloc( jpi, jpj, zeeu, zeev, zlttl )
113      CALL wrk_alloc( jpi, jpj, jpk, ztutl, ztvtl )
114      !
115      IF( kt == kit000 ) THEN
116         IF(lwp) WRITE(numout,*)
117         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap_tan : iso-level biharmonic operator on ', cdtype
118         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
119      ENDIF
120
121      DO jn = 1, kjpt
122      !                                                ! ===============
123         DO jk = 1, jpkm1                                 ! Horizontal slab
124            !                                             ! ===============
125
126            ! 0. Initialization of metric arrays (for z- or s-coordinates)
127            ! ----------------------------------
128
129            DO jj = 1, jpjm1
130               DO ji = 1, fs_jpim1   ! vector opt.
131                  zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
132                  zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
133               END DO
134            END DO
135
136            ! 1. Laplacian
137            ! ------------
138
139            ! First derivative (gradient)
140            DO jj = 1, jpjm1
141               DO ji = 1, fs_jpim1   ! vector opt.
142                  ztutl(ji,jj,jk) = zeeu(ji,jj) * ( ptb_tl(ji+1,jj  ,jk,jn) - ptb_tl(ji,jj,jk,jn) )
143                  ztvtl(ji,jj,jk) = zeev(ji,jj) * ( ptb_tl(ji  ,jj+1,jk,jn) - ptb_tl(ji,jj,jk,jn) )
144               END DO
145            END DO
146            IF( ln_zps ) THEN      ! set gradient at partial step level
147               DO jj = 1, jpjm1
148                  DO ji = 1, jpim1
149                     IF( mbku(ji,jj) == jk )  ztutl(ji,jj,jk) = zeeu(ji,jj) * pgu_tl(ji,jj,jn)
150                     IF( mbkv(ji,jj) == jk )  ztvtl(ji,jj,jk) = zeev(ji,jj) * pgv_tl(ji,jj,jn)
151                  END DO
152               END DO
153            ENDIF
154
155            ! Second derivative (divergence)
156            DO jj = 2, jpjm1
157               DO ji = fs_2, fs_jpim1   ! vector opt.
158                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
159                  zlttl(ji,jj) = fsahtt(ji,jj,jk) * zbtr * (  ztutl(ji,jj,jk) - ztutl(ji-1,jj,jk)   &
160                               &                            + ztvtl(ji,jj,jk) - ztvtl(ji,jj-1,jk)  )
161               END DO
162            END DO
163
164            ! Lateral boundary conditions on the laplacian (zlttl,zlstl)   (unchanged sgn)
165            CALL lbc_lnk( zlttl, 'T', 1.0_wp )
166
167            ! 2. Bilaplacian
168            ! --------------
169
170            ! third derivative (gradient)
171            DO jj = 1, jpjm1
172               DO ji = 1, fs_jpim1   ! vector opt.
173                  ztutl(ji,jj,jk) = zeeu(ji,jj) * ( zlttl(ji+1,jj  ) - zlttl(ji,jj) )
174                  ztvtl(ji,jj,jk) = zeev(ji,jj) * ( zlttl(ji  ,jj+1) - zlttl(ji,jj) )
175               END DO
176            END DO
177
178            ! fourth derivative (divergence) and add to the general tracer trend
179            DO jj = 2, jpjm1
180               DO ji = fs_2, fs_jpim1   ! vector opt.
181                  ! horizontal diffusive trends
182                  zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
183                  ztatl = zbtr * (  ztutl(ji,jj,jk) - ztutl(ji-1,jj,jk) + ztvtl(ji,jj,jk) - ztvtl(ji,jj-1,jk)  )
184                  ! add it to the general tracer trends
185                  pta_tl(ji,jj,jk,jn) = pta_tl(ji,jj,jk,jn) + ztatl
186               END DO
187            END DO
188            !                                             ! ===============
189         END DO                                           ! Horizontal slab
190         !                                                ! ===============
191      END DO
192      IF( nn_timing == 1 )  CALL timing_stop( 'tra_ldf_bilap_tan')
193      !
194      CALL wrk_dealloc( jpi, jpj, zeeu, zeev, zlttl )
195      CALL wrk_dealloc( jpi, jpj, jpk, ztutl, ztvtl )
196      !
197   END SUBROUTINE tra_ldf_bilap_tan
198
199
200   SUBROUTINE tra_ldf_bilap_adj( kt, kit000, cdtype, pgu_ad, pgv_ad,      &
201      &                                  ptb_ad, pta_ad, kjpt   )
202      !!----------------------------------------------------------------------
203      !!                  ***  ROUTINE tra_ldf_bilap  ***
204      !!
205      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
206      !!      trend and add it to the general trend of tracer equation.
207      !!
208      !! ** Method  :   4th order diffusive operator along model level surfaces
209      !!      evaluated using before fields (forward time scheme). The hor.
210      !!      diffusive trends of temperature (idem for salinity) is given by:
211      !!      Laplacian of tb:
212      !!         zlt   = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(tb) ]
213      !!                                  + dj-1[ e1v*e3v/e2v dj(tb) ]  }
214      !!      Multiply by the eddy diffusivity coef. and insure lateral bc:
215      !!        zlt   = ahtt * zlt
216      !!        call to lbc_lnk
217      !!      Bilaplacian (laplacian of zlt):
218      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ e2u*e3u/e1u di(zlt) ]
219      !!                                  + dj-1[ e1v*e3v/e2v dj(zlt) ]  }
220      !!      Note: if key_zco defined, e3t=e3u=e3v, they are simplified.
221      !!
222      !!      Add this trend to the general trend (ta,sa):
223      !!         (ta,sa) = (ta,sa) + ( difft , diffs )
224      !!
225      !! ** Action : - Update (ta,sa) arrays with the before iso-level
226      !!               biharmonic mixing trend.
227      !!
228      !! History :
229      !!        !  91-11  (G. Madec)  Original code
230      !!        !  93-03  (M. Guyon)  symetrical conditions
231      !!        !  95-11  (G. Madec)  suppress volumetric scale factors
232      !!        !  96-01  (G. Madec)  statement function for e3
233      !!        !  96-01  (M. Imbard)  mpp exchange
234      !!        !  97-07  (G. Madec)  optimization, and ahtt
235      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
236      !!   9.0  !  04-08  (C. talandier) New trends organization
237      !!        !  05-11  (G. Madec)  zps or sco as default option
238      !!----------------------------------------------------------------------
239      !! History of the tangent routine
240      !!   9.0  !  10-05 (P.A. Bouttier) tangent of 9.0
241      !!----------------------------------------------------------------------
242      !! * Modules used
243      !! * Arguments
244      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
245      INTEGER                              , INTENT(in   ) ::   kit000           ! first time step index
246      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype           ! =TRA or TRC (tracer indicator)
247      INTEGER                              , INTENT(in   ) ::   kjpt             ! number of tracers
248      REAL(wp), DIMENSION(jpi,jpj,    kjpt), INTENT(inout) ::   pgu_ad, pgv_ad   ! tracer gradient at pstep levels
249      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   ptb_ad           ! before and now tracer fields
250      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta_ad           ! tracer trend
251
252      !! * Local declarations
253      INTEGER ::   ji, jj, jk, jn             ! dummy loop indices
254      INTEGER ::   iku, ikv               ! temporary integers
255      REAL(wp) ::  ztaad, zsaad, ztmp, zbtr     ! temporary scalars
256      REAL(wp), POINTER, DIMENSION(:,:) ::   &
257         zeeu, zeev,              & ! 2D workspace
258         zltad
259      REAL(wp), POINTER, DIMENSION(:,:,:) ::   &
260         ztuad, ztvad                         ! 3D workspace
261      !!----------------------------------------------------------------------
262      !
263      IF( nn_timing == 1 )  CALL timing_start( 'tra_ldf_bilap_adj')
264      !
265      CALL wrk_alloc( jpi, jpj, zeeu, zeev, zltad )
266      CALL wrk_alloc( jpi, jpj, jpk, ztuad, ztvad )
267      !
268      IF( kt == nit000 ) THEN
269         IF(lwp) WRITE(numout,*)
270         IF(lwp) WRITE(numout,*) 'tra_ldf_bilap_adj : iso-level biharmonic operator on ', cdtype
271         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
272      ENDIF
273      ztaad = 0.e0_wp
274      zsaad = 0.e0_wp
275      ztuad = 0.e0_wp
276      ztvad = 0.e0_wp
277
278      DO jn = 1, kjpt
279         DO jk = jpkm1, 1, -1
280               ! 0. Initialization of metric arrays (for z- or s-coordinates)
281               ! ----------------------------------
282                  DO jj = jpjm1, 1, -1
283                     DO ji = fs_jpim1, 1 ,-1   ! vector opt.
284                        zeeu(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
285                        zeev(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
286                     END DO
287                  END DO
288               !
289               ! 2. Bilaplacian
290               ! --------------
291               ! fourth derivative (divergence) and add to the general tracer trend
292               DO jj = jpjm1, 2, -1
293                  DO ji = fs_jpim1, fs_2, -1   ! vector opt.
294                     zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
295                     ! add it to the general tracer trends
296                     ztaad = pta_ad(ji,jj,jk,jn) * zbtr
297                     ! horizontal diffusive trends
298                     ztuad(ji  ,jj  ,jk) = ztuad(ji  ,jj  ,jk) + ztaad
299                     ztuad(ji-1,jj  ,jk) = ztuad(ji-1,jj  ,jk) - ztaad
300                     ztvad(ji  ,jj  ,jk) = ztvad(ji  ,jj  ,jk) + ztaad
301                     ztvad(ji  ,jj-1,jk) = ztvad(ji  ,jj-1,jk) - ztaad
302                  END DO
303               END DO
304
305               ! third derivative (gradient)
306               DO jj = jpjm1, 1, -1
307                  DO ji = fs_jpim1, 1 ,-1   ! vector opt.
308                     ztmp = zeev(ji,jj) * ztvad(ji,jj,jk)
309                     zltad(ji  ,jj+1) = zltad(ji  ,jj+1) + ztmp
310                     zltad(ji  ,jj  ) = zltad(ji  ,jj  ) - ztmp
311                     ztmp = zeeu(ji,jj) * ztuad(ji,jj,jk)
312                     zltad(ji+1,jj  ) = zltad(ji+1,jj  ) + ztmp
313                     zltad(ji  ,jj  ) = zltad(ji  ,jj  ) - ztmp
314                     ztuad(ji  ,jj  ,jk) = 0.0_wp
315                     ztvad(ji  ,jj  ,jk) = 0.0_wp
316                  END DO
317               END DO
318               ! Lateral boundary conditions on the laplacian (zltad,zlsad)   (unchanged sgn)
319               CALL lbc_lnk_adj( zltad, 'T', 1.0_wp )
320               ! Multiply by the eddy diffusivity coefficient
321               DO jj = jpjm1, 2, -1
322                  DO ji = fs_jpim1, fs_2, -1   ! vector opt.
323                     zltad(ji,jj) = fsahtt(ji,jj,jk) * zltad(ji,jj)
324                  END DO
325               END DO
326               ! Second derivative (divergence)
327               DO jj = jpjm1, 2, -1
328                  DO ji = fs_jpim1, fs_2, -1   ! vector opt.
329                     zbtr = 1.0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
330                     ztmp = zbtr * zltad(ji,jj)
331                     ztvad(ji  ,jj-1,jk) = ztvad(ji  ,jj-1,jk) - ztmp
332                     ztvad(ji  ,jj  ,jk) = ztvad(ji  ,jj  ,jk) + ztmp
333                     ztuad(ji-1,jj  ,jk) = ztuad(ji-1,jj  ,jk) - ztmp
334                     ztuad(ji  ,jj  ,jk) = ztuad(ji  ,jj  ,jk) + ztmp
335                     zltad(ji,jj) = 0.0_wp
336                  END DO
337               END DO
338               IF( ln_zps ) THEN      ! set gradient at partial step level
339               DO jj = jpjm1, 1, -1
340                  DO ji = fs_jpim1, 1 ,-1   ! vector opt.
341                        ! last level
342                        IF( mbku(ji,jj) == jk ) THEN
343                           pgu_ad(ji,jj,jn)   = pgu_ad(ji,jj,jn) + zeeu(ji,jj) * ztuad(ji,jj,jk)
344                           ztuad(ji,jj,jk) = 0.0_wp
345                        ENDIF
346                        IF( mbkv(ji,jj) == jk ) THEN
347                           pgv_ad(ji,jj,jn)   = pgv_ad(ji,jj,jn) + zeev(ji,jj) * ztvad(ji,jj,jk)
348                           ztvad(ji,jj,jk) = 0.0_wp
349                        ENDIF
350                     END DO
351                  END DO
352               ENDIF
353               ! 1. Laplacian
354               ! ------------
355
356               ! First derivative (gradient)
357               DO jj = jpjm1, 1, -1
358                  DO ji = fs_jpim1, 1 ,-1   ! vector opt.
359                     ztmp = zeev(ji,jj) * ztvad(ji,jj,jk)
360                     ptb_ad(ji  ,jj  ,jk,jn) = ptb_ad(ji  ,jj  ,jk,jn) - ztmp
361                     ptb_ad(ji  ,jj+1,jk,jn) = ptb_ad(ji  ,jj+1,jk,jn) + ztmp
362                     ztmp = zeeu(ji,jj) * ztuad(ji,jj,jk)
363                     ptb_ad(ji  ,jj  ,jk,jn) = ptb_ad(ji  ,jj  ,jk,jn) - ztmp
364                     ptb_ad(ji+1,jj  ,jk,jn) = ptb_ad(ji+1,jj  ,jk,jn) + ztmp
365                     ztuad(ji,jj,jk) = 0.0_wp
366                     ztvad(ji,jj,jk) = 0.0_wp
367                  END DO
368               END DO
369          END DO
370      END DO
371      !
372      CALL wrk_dealloc( jpi, jpj, zeeu, zeev, zltad )
373      CALL wrk_dealloc( jpi, jpj, jpk, ztuad, ztvad )
374      !
375      IF( nn_timing == 1 )  CALL timing_stop( 'tra_ldf_bilap_adj')
376      !
377   END SUBROUTINE tra_ldf_bilap_adj
378
379   SUBROUTINE tra_ldf_bilap_adj_tst ( kumadt )
380      !!-----------------------------------------------------------------------
381      !!
382      !!                  ***  ROUTINE example_adj_tst ***
383      !!
384      !! ** Purpose : Test the adjoint routine.
385      !!
386      !! ** Method  : Verify the scalar product
387      !!
388      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
389      !!
390      !!              where  L   = tangent routine
391      !!                     L^T = adjoint routine
392      !!                     W   = diagonal matrix of scale factors
393      !!                     dx  = input perturbation (random field)
394      !!                     dy  = L dx
395      !!
396      !! History :
397      !!        ! 08-08 (A. Vidard)
398      !!-----------------------------------------------------------------------
399      !! * Modules used
400
401      !! * Arguments
402      INTEGER, INTENT(IN) :: &
403         & kumadt             ! Output unit
404
405      !! * Local declarations
406      INTEGER ::  &
407         & ji,    &        ! dummy loop indices
408         & jj,    &
409         & jk
410      INTEGER, DIMENSION(jpi,jpj) :: &
411         & iseed_2d        ! 2D seed for the random number generator
412      REAL(KIND=wp) :: &
413         & zsp1,         & ! scalar product involving the tangent routine
414         & zsp1_T,       &
415         & zsp1_S,       &
416         & zsp2,         & ! scalar product involving the adjoint routine
417         & zsp2_1,       &
418         & zsp2_2,       &
419         & zsp2_3,       &
420         & zsp2_4,       &
421         & zsp2_5,       &
422         & zsp2_6,       &
423         & zsp2_7,       &
424         & zsp2_8,       &
425         & zsp2_T,       &
426         & zsp2_S
427      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
428         & ztb_tlin ,      & ! Tangent input
429         & zsb_tlin ,      & ! Tangent input
430         & zta_tlin ,      & ! Tangent input
431         & zsa_tlin ,      & ! Tangent input
432         & zta_tlout,      & ! Tangent output
433         & zsa_tlout,      & ! Tangent output
434         & zta_adin,       & ! Adjoint input
435         & zsa_adin,       & ! Adjoint input
436         & ztb_adout ,     & ! Adjoint output
437         & zsb_adout ,     & ! Adjoint output
438         & zta_adout ,     & ! Adjoint output
439         & zsa_adout ,     & ! Adjoint output
440         & z3r               ! 3D random field
441      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
442         & zgtu_tlin ,     & ! Tangent input
443         & zgsu_tlin ,     & ! Tangent input
444         & zgtv_tlin ,     & ! Tangent input
445         & zgsv_tlin ,     & ! Tangent input
446         & zgtu_adout ,    & ! Adjoint output
447         & zgsu_adout ,    & ! Adjoint output
448         & zgtv_adout ,    & ! Adjoint output
449         & zgsv_adout ,    & ! Adjoint output
450         & z2r               ! 2D random field
451      CHARACTER(LEN=14) :: cl_name
452      ! Allocate memory
453
454      ALLOCATE( &
455         & ztb_tlin(jpi,jpj,jpk),      &
456         & zsb_tlin(jpi,jpj,jpk),      &
457         & zta_tlin(jpi,jpj,jpk),      &
458         & zsa_tlin(jpi,jpj,jpk),      &
459         & zgtu_tlin(jpi,jpj),         &
460         & zgsu_tlin(jpi,jpj),         &
461         & zgtv_tlin(jpi,jpj),         &
462         & zgsv_tlin(jpi,jpj),         &
463         & zta_tlout(jpi,jpj,jpk),     &
464         & zsa_tlout(jpi,jpj,jpk),     &
465         & zta_adin(jpi,jpj,jpk),      &
466         & zsa_adin(jpi,jpj,jpk),      &
467         & ztb_adout(jpi,jpj,jpk),     &
468         & zsb_adout(jpi,jpj,jpk),     &
469         & zta_adout(jpi,jpj,jpk),     &
470         & zsa_adout(jpi,jpj,jpk),     &
471         & zgtu_adout(jpi,jpj),        &
472         & zgsu_adout(jpi,jpj),        &
473         & zgtv_adout(jpi,jpj),        &
474         & zgsv_adout(jpi,jpj),        &
475         & z3r(jpi,jpj,jpk),           &
476         & z2r(jpi,jpj)                &
477         & )
478
479
480      !=======================================================================
481      ! 1) dx = ( tb_tl, ta_tl, sb_tl, sa_tl, gtu_tl, gtv_tl, gsu_tl, gsv_tl )
482      !    dy = ( ta_tl, sa_tl )
483      !=======================================================================
484
485      !--------------------------------------------------------------------
486      ! Reset the tangent and adjoint variables
487      !--------------------------------------------------------------------
488
489      ztb_tlin(:,:,:)  = 0.0_wp
490      zsb_tlin(:,:,:)  = 0.0_wp
491      zta_tlin(:,:,:)  = 0.0_wp
492      zsa_tlin(:,:,:)  = 0.0_wp
493      zgtu_tlin(:,:)   = 0.0_wp
494      zgsu_tlin(:,:)   = 0.0_wp
495      zgtv_tlin(:,:)   = 0.0_wp
496      zgsv_tlin(:,:)   = 0.0_wp
497      zta_tlout(:,:,:) = 0.0_wp
498      zsa_tlout(:,:,:) = 0.0_wp
499      zta_adin(:,:,:)  = 0.0_wp
500      zsa_adin(:,:,:)  = 0.0_wp
501      ztb_adout(:,:,:) = 0.0_wp
502      zsb_adout(:,:,:) = 0.0_wp
503      zta_adout(:,:,:) = 0.0_wp
504      zsa_adout(:,:,:) = 0.0_wp
505      zgtu_adout(:,:)  = 0.0_wp
506      zgsu_adout(:,:)  = 0.0_wp
507      zgtv_adout(:,:)  = 0.0_wp
508      zgsv_adout(:,:)  = 0.0_wp
509
510      tsb_tl(:,:,:,:) = 0.0_wp
511      tsa_tl(:,:,:,:) = 0.0_wp
512      gtsu_tl(:,:,:)  = 0.0_wp
513      gtsv_tl(:,:,:)  = 0.0_wp
514      tsb_ad(:,:,:,:) = 0.0_wp
515      tsa_ad(:,:,:,:) = 0.0_wp
516      gtsu_ad(:,:,:)  = 0.0_wp
517      gtsv_ad(:,:,:)  = 0.0_wp
518
519      !--------------------------------------------------------------------
520      ! Initialize the tangent input with random noise: dx
521      !--------------------------------------------------------------------
522
523      CALL grid_random(  z3r, 'T', 0.0_wp, stdt )
524      DO jk = 1, jpk
525        DO jj = nldj, nlej
526           DO ji = nldi, nlei
527              ztb_tlin(ji,jj,jk) = z3r(ji,jj,jk)
528            END DO
529         END DO
530      END DO
531
532      CALL grid_random(  z3r, 'T', 0.0_wp, stds )
533      DO jk = 1, jpk
534        DO jj = nldj, nlej
535           DO ji = nldi, nlei
536              zsb_tlin(ji,jj,jk) = z3r(ji,jj,jk)
537            END DO
538         END DO
539      END DO
540
541      CALL grid_random(  z3r, 'T', 0.0_wp, stdt )
542      DO jk = 1, jpk
543        DO jj = nldj, nlej
544           DO ji = nldi, nlei
545              zta_tlin(ji,jj,jk) = z3r(ji,jj,jk)
546            END DO
547         END DO
548      END DO
549
550      CALL grid_random(  z3r, 'T', 0.0_wp, stds )
551      DO jk = 1, jpk
552        DO jj = nldj, nlej
553           DO ji = nldi, nlei
554              zsa_tlin(ji,jj,jk) = z3r(ji,jj,jk)
555            END DO
556         END DO
557      END DO
558
559      CALL grid_random(  z2r, 'U', 0.0_wp, stds )
560      DO jj = nldj, nlej
561         DO ji = nldi, nlei
562            zgtu_tlin(ji,jj) = z2r(ji,jj)
563         END DO
564      END DO
565
566      CALL grid_random(  z2r, 'U', 0.0_wp, stds )
567      DO jj = nldj, nlej
568         DO ji = nldi, nlei
569            zgsu_tlin(ji,jj) = z2r(ji,jj)
570         END DO
571      END DO
572
573      CALL grid_random(  z2r, 'V', 0.0_wp, stds )
574      DO jj = nldj, nlej
575         DO ji = nldi, nlei
576            zgtv_tlin(ji,jj) = z2r(ji,jj)
577         END DO
578      END DO
579
580      CALL grid_random(  z2r, 'V', 0.0_wp, stds )
581      DO jj = nldj, nlej
582         DO ji = nldi, nlei
583            zgsv_tlin(ji,jj) = z2r(ji,jj)
584         END DO
585      END DO
586
587      tsb_tl(:,:,:,jp_tem) = ztb_tlin(:,:,:)
588      tsb_tl(:,:,:,jp_sal) = zsb_tlin(:,:,:)
589      tsa_tl(:,:,:,jp_tem) = zta_tlin(:,:,:)
590      tsa_tl(:,:,:,jp_sal) = zsa_tlin(:,:,:)
591      gtsu_tl(:,:,jp_tem)  = zgtu_tlin(:,:)
592      gtsu_tl(:,:,jp_sal)  = zgsu_tlin(:,:)
593      gtsv_tl(:,:,jp_tem)  = zgtv_tlin(:,:)
594      gtsv_tl(:,:,jp_sal)  = zgsv_tlin(:,:)
595
596      CALL tra_ldf_bilap_tan( nit000, nit000, 'TRA', gtsu_tl, gtsv_tl,tsb_tl, tsa_tl, jpts  )
597
598      zta_tlout(:,:,:) = tsa_tl(:,:,:,jp_tem)
599      zsa_tlout(:,:,:) = tsa_tl(:,:,:,jp_sal)
600
601      !--------------------------------------------------------------------
602      ! Initialize the adjoint variables: dy^* = W dy
603      !--------------------------------------------------------------------
604
605      DO jk = 1, jpk
606        DO jj = nldj, nlej
607           DO ji = nldi, nlei
608              zsa_adin(ji,jj,jk) = zsa_tlout(ji,jj,jk) &
609                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
610                 &               * tmask(ji,jj,jk) * wesp_s(jk)
611              zta_adin(ji,jj,jk) = zta_tlout(ji,jj,jk) &
612                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
613                 &               * tmask(ji,jj,jk) * wesp_t(jk)
614            END DO
615         END DO
616      END DO
617
618      !--------------------------------------------------------------------
619      ! Compute the scalar product: ( L dx )^T W dy
620      !-------------------------------------------------------------------
621
622      zsp1_T = DOT_PRODUCT( zta_tlout, zta_adin )
623      zsp1_S = DOT_PRODUCT( zsa_tlout, zsa_adin )
624      zsp1 = zsp1_T + zsp1_S
625
626      !--------------------------------------------------------------------
627      ! Call the adjoint routine: dx^* = L^T dy^*
628      !--------------------------------------------------------------------
629
630      tsa_ad(:,:,:,jp_tem) = zta_adin(:,:,:)
631      tsa_ad(:,:,:,jp_sal) = zsa_adin(:,:,:)
632
633      CALL tra_ldf_bilap_adj( nit000 , nit000, 'TRA', gtsu_ad, gtsv_ad,tsb_ad, tsa_ad, jpts)
634
635      ztb_adout(:,:,:) = tsb_ad(:,:,:,jp_tem)
636      zsb_adout(:,:,:) = tsb_ad(:,:,:,jp_sal)
637      zta_adout(:,:,:) = tsa_ad(:,:,:,jp_tem)
638      zsa_adout(:,:,:) = tsa_ad(:,:,:,jp_sal)
639      zgtu_adout(:,:)  = gtsu_ad(:,:,jp_tem)
640      zgsu_adout(:,:)  = gtsu_ad(:,:,jp_sal)
641      zgtv_adout(:,:)  = gtsv_ad(:,:,jp_tem)
642      zgsv_adout(:,:)  = gtsv_ad(:,:,jp_sal)
643
644      zsp2_1 = DOT_PRODUCT( ztb_tlin , ztb_adout  )
645      zsp2_2 = DOT_PRODUCT( zta_tlin , zta_adout  )
646      zsp2_3 = DOT_PRODUCT( zgtu_tlin, zgtu_adout )
647      zsp2_4 = DOT_PRODUCT( zgtv_tlin, zgtv_adout )
648      zsp2_5 = DOT_PRODUCT( zsb_tlin , zsb_adout  )
649      zsp2_6 = DOT_PRODUCT( zsa_tlin , zsa_adout  )
650      zsp2_7 = DOT_PRODUCT( zgsu_tlin, zgsu_adout )
651      zsp2_8 = DOT_PRODUCT( zgsv_tlin, zgsv_adout )
652
653      zsp2_T = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4
654      zsp2_S = zsp2_5 + zsp2_6 + zsp2_7 + zsp2_8
655      zsp2   = zsp2_T + zsp2_S
656
657      cl_name = 'tra_ldf_bilap'
658      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
659
660      DEALLOCATE(         &
661         & ztb_tlin,      & ! Tangent input
662         & zsb_tlin,      & ! Tangent input
663         & zta_tlin,      & ! Tangent input
664         & zsa_tlin,      & ! Tangent input
665         & zgtu_tlin,     & ! Tangent input
666         & zgsu_tlin,     & ! Tangent input
667         & zgtv_tlin,     & ! Tangent input
668         & zgsv_tlin,     & ! Tangent input
669         & zta_tlout,     & ! Tangent output
670         & zsa_tlout,     & ! Tangent output
671         & zta_adin,      & ! Adjoint input
672         & zsa_adin,      & ! Adjoint input
673         & ztb_adout,     & ! Adjoint output
674         & zsb_adout,     & ! Adjoint output
675         & zta_adout,     & ! Adjoint output
676         & zsa_adout,     & ! Adjoint output
677         & zgtu_adout,    & ! Adjoint output
678         & zgsu_adout,    & ! Adjoint output
679         & zgtv_adout,    & ! Adjoint output
680         & zgsv_adout,    & ! Adjoint output
681         & z3r,           & ! 3D random field
682         & z2r            &
683         & )
684
685
686   END SUBROUTINE tra_ldf_bilap_adj_tst
687
688#endif
689END MODULE traldf_bilap_tam
Note: See TracBrowser for help on using the repository browser.