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 @ 3611

Last change on this file since 3611 was 3611, checked in by pabouttier, 11 years ago

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