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_lap_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_lap_tam.F90 @ 3640

Last change on this file since 3640 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: 25.3 KB
Line 
1MODULE traldf_lap_tam
2#if defined key_tam
3   !!==============================================================================
4   !!                       ***  MODULE  traldf_lap  ***
5   !! Ocean active tracers:  horizontal component of the lateral tracer mixing trend
6   !!                        Tangent and adjoint module
7   !!==============================================================================
8   !! History :   9.0  !  ?????
9   !! History of the T&A module:
10   !!             9.0  !  09-03  (F. Vigilant) original version
11   !!       NEMO  3.4  ! 12-07   (P.-A. Bouttier)
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   tra_ldf_lap  : update the tracer trend with the horizontal diffusion
16   !!                 using a iso-level harmonic (laplacien) operator.
17   !!----------------------------------------------------------------------
18   !! * Modules used
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 timing
31   USE wrk_nemo
32
33   IMPLICIT NONE
34   PRIVATE
35
36   !! * Routine accessibility
37   PUBLIC tra_ldf_lap_tan      ! routine called by tradldf_tam.F90
38   PUBLIC tra_ldf_lap_adj      ! routine called by tradldf_tam.F90
39   PUBLIC tra_ldf_lap_adj_tst  ! routine called by tradldf_tam.F90
40
41   REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) ::   e1ur, e2vr   ! scale factor coefficients
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "ldftra_substitute.h90"
46#  include "vectopt_loop_substitute.h90"
47   !!----------------------------------------------------------------------
48   !!   OPA 9.0 , LOCEAN-IPSL (2005)
49   !! $Id: traldf_lap.F90 1152 2008-06-26 14:11:13Z rblod $
50   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
51   !!----------------------------------------------------------------------
52
53CONTAINS
54
55   SUBROUTINE tra_ldf_lap_tan( kt, kit000, cdtype, pgu_tl, pgv_tl,      &
56      &                                ptb_tl, pta_tl, kjpt )
57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE tra_ldf_lap_tan  ***
59      !!
60      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
61      !!      trend and add it to the general trend of tracer equation.
62      !!
63      !! ** Method  :   Second order diffusive operator evaluated using before
64      !!      fields (forward time scheme). The horizontal diffusive trends of
65      !!      temperature (idem for salinity) is given by:
66      !!          difft = 1/(e1t*e2t*e3t) {  di-1[ aht e2u*e3u/e1u di(tb) ]
67      !!                                   + dj-1[ aht e1v*e3v/e2v dj(tb) ] }
68      !!     Note: key_zco defined, the e3t=e3u=e3v, the trend becomes:
69      !!          difft = 1/(e1t*e2t) {  di-1[ aht e2u/e1u di(tb) ]
70      !!                               + dj-1[ aht e1v/e2v dj(tb) ] }
71      !!      Add this trend to the general tracer trend (ta,sa):
72      !!          (ta,sa) = (ta,sa) + ( difft , diffs )
73      !!
74      !! ** Action  : - Update (ta,sa) arrays with the before iso-level
75      !!                harmonic mixing trend.
76      !!
77      !! History :
78      !!   1.0  !  87-06  (P. Andrich, D. L Hostis)  Original code
79      !!        !  91-11  (G. Madec)
80      !!        !  95-11  (G. Madec)  suppress volumetric scale factors
81      !!        !  96-01  (G. Madec)  statement function for e3
82      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
83      !!   9.0  !  04-08  (C. Talandier) New trends organization
84      !!        !  05-11  (G. Madec)  add zps case
85      !! History of the tangent routine
86      !!   9.0  !  03-09 (F. Vigilant) tangent of 9.0
87      !!----------------------------------------------------------------------
88
89      !! * Arguments
90      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
91      INTEGER                              , INTENT(in   ) ::   kit000           ! first time step index
92      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype           ! =TRA or TRC (tracer indicator)
93      INTEGER                              , INTENT(in   ) ::   kjpt             ! number of tracers
94      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(in   ) ::   pgu_tl, pgv_tl   ! tracer gradient at pstep levels
95      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb_tl           ! before and now tracer fields
96      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta_tl           ! tracer trend
97      !! * Local declarations
98      INTEGER  ::   ji, jj, jk, jn, ierr         ! dummy loop indices
99      INTEGER  ::   iku, ikv               ! temporary integers
100      REAL(wp) ::  zabe1, zabe2, zbtr      ! temporary scalars
101      REAL(wp) ::  ztatl, zsatl            ! temporary scalars
102      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztutl, ztvtl     ! 3D workspace
103      !!----------------------------------------------------------------------
104      !
105      CALL timing_start('tra_ldf_lap_tan')
106      !
107      CALL wrk_alloc(jpi,jpj,jpk, ztutl, ztvtl)
108      !
109      IF( kt == nit000 ) THEN
110         IF(lwp) WRITE(numout,*)
111         IF(lwp) WRITE(numout,*) 'tra_ldf_lap_tan : iso-level laplacian diffusion on ', cdtype
112         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
113         IF( .NOT. ALLOCATED( e1ur ) ) THEN
114            ! This routine may be called for both active and passive tracers.
115            ! Allocate and set saved arrays on first call only.
116            ALLOCATE( e1ur(jpi,jpj), e2vr(jpi,jpj), STAT=ierr )
117            IF( lk_mpp    )   CALL mpp_sum( ierr )
118            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'tra_ldf_lap : unable to allocate arrays' )
119            !
120            e1ur(:,:) = e2u(:,:) / e1u(:,:)
121            e2vr(:,:) = e1v(:,:) / e2v(:,:)
122         ENDIF
123      ENDIF
124
125      !
126      DO jn = 1, kjpt
127         !                                                  ! =============
128         DO jk = 1, jpkm1                                   ! Vertical slab
129            !                                               ! =============
130            ! 1. First derivative (gradient)
131            ! -------------------
132            DO jj = 1, jpjm1
133               DO ji = 1, fs_jpim1   ! vector opt.
134                  zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * e1ur(ji,jj) * fse3u(ji,jj,jk)
135                  zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * e2vr(ji,jj) * fse3v(ji,jj,jk)
136                  ztutl(ji,jj,jk) = zabe1 * ( ptb_tl(ji+1,jj  ,jk,jn) - ptb_tl(ji,jj,jk,jn) )
137                  ztvtl(ji,jj,jk) = zabe2 * ( ptb_tl(ji  ,jj+1,jk,jn) - ptb_tl(ji,jj,jk,jn) )
138               END DO
139            END DO
140            IF( ln_zps ) THEN      ! set gradient at partial step level
141               DO jj = 1, jpjm1
142                  DO ji = 1, fs_jpim1   ! vector opt.
143                     ! last level
144                     iku = mbku(ji,jj)
145                     ikv = mbku(ji,jj)
146                     IF( iku == jk ) THEN
147                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e1ur(ji,jj) * fse3u(ji,jj,iku)
148                        ztutl(ji,jj,jk) = zabe1 * pgu_tl(ji,jj,jn)
149                     ENDIF
150                     IF( ikv == jk ) THEN
151                        zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e2vr(ji,jj) * fse3v(ji,jj,ikv)
152                        ztvtl(ji,jj,jk) = zabe2 * pgv_tl(ji,jj,jn)
153                     ENDIF
154                  END DO
155               END DO
156            ENDIF
157
158
159            ! 2. Second derivative (divergence)
160            ! --------------------
161            DO jj = 2, jpjm1
162               DO ji = fs_2, fs_jpim1   ! vector opt.
163                  zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk))
164                  ! horizontal diffusive trends
165                  ztatl = zbtr * (  ztutl(ji,jj,jk) - ztutl(ji-1,jj  ,jk)   &
166                     &            + ztvtl(ji,jj,jk) - ztvtl(ji  ,jj-1,jk)  )
167                  ! add it to the general tracer trends
168                  pta_tl(ji,jj,jk,jn) = pta_tl(ji,jj,jk,jn) + ztatl
169               END DO
170            END DO
171            !                                               ! =============
172         END DO                                             !  End of slab
173         !                                                  ! =============
174      END DO
175      !
176      CALL timing_stop('tra_ldf_lap_tan')
177      !
178      CALL wrk_dealloc(jpi,jpj,jpk, ztutl, ztvtl)
179      !
180
181   END SUBROUTINE tra_ldf_lap_tan
182
183
184   SUBROUTINE tra_ldf_lap_adj( kt, kit000, cdtype, pgu_ad, pgv_ad,      &
185      &                                ptb_ad, pta_ad, kjpt )
186      !!----------------------------------------------------------------------
187      !!                  ***  ROUTINE tra_ldf_lap_adj  ***
188      !!
189      !! ** Purpose :   Compute the before horizontal tracer (t & s) diffusive
190      !!      trend and add it to the general trend of tracer equation.
191      !!
192      !! ** Method  :   Second order diffusive operator evaluated using before
193      !!      fields (forward time scheme). The horizontal diffusive trends of
194      !!      temperature (idem for salinity) is given by:
195      !!          difft = 1/(e1t*e2t*e3t) {  di-1[ aht e2u*e3u/e1u di(tb) ]
196      !!                                   + dj-1[ aht e1v*e3v/e2v dj(tb) ] }
197      !!     Note: key_zco defined, the e3t=e3u=e3v, the trend becomes:
198      !!          difft = 1/(e1t*e2t) {  di-1[ aht e2u/e1u di(tb) ]
199      !!                               + dj-1[ aht e1v/e2v dj(tb) ] }
200      !!      Add this trend to the general tracer trend (ta,sa):
201      !!          (ta,sa) = (ta,sa) + ( difft , diffs )
202      !!
203      !! ** Action  : - Update (ta,sa) arrays with the before iso-level
204      !!                harmonic mixing trend.
205      !!
206      !! History :
207      !!   1.0  !  87-06  (P. Andrich, D. L Hostis)  Original code
208      !!        !  91-11  (G. Madec)
209      !!        !  95-11  (G. Madec)  suppress volumetric scale factors
210      !!        !  96-01  (G. Madec)  statement function for e3
211      !!   8.5  !  02-06  (G. Madec)  F90: Free form and module
212      !!   9.0  !  04-08  (C. Talandier) New trends organization
213      !!        !  05-11  (G. Madec)  add zps case
214      !! History of the tangent routine
215      !!   9.0  !  03-09 (F. Vigilant) adjoint of 9.0
216      !!----------------------------------------------------------------------
217
218      !! * Arguments
219      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
220            INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
221      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
222      INTEGER                              , INTENT(in   ) ::   kjpt       ! number of tracers
223      REAL(wp), DIMENSION(jpi,jpj    ,kjpt), INTENT(inout) ::   pgu_ad, pgv_ad   ! tracer gradient at pstep levels
224      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   ptb_ad        ! before and now tracer fields
225      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta_ad        ! tracer trend
226      !! * Local declarations
227      INTEGER ::   ji, jj, jk, jn             ! dummy loop indices
228      INTEGER ::   iku, ikv, ierr               ! temporary integers
229      REAL(wp) ::  zabe1, zabe2, zbtr     ! temporary scalars
230      REAL(wp) ::  ztaad, zsaad           ! temporary scalars
231      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztuad, ztvad    ! 3D workspace
232      !!----------------------------------------------------------------------
233      !
234      CALL timing_start('tra_ldf_lap_adj')
235      !
236      CALL wrk_alloc(jpi,jpj,jpk, ztuad, ztvad)
237      !
238      ztvad(:,:,:) = 0.0_wp
239      ztuad(:,:,:) = 0.0_wp
240      ztaad        = 0.0_wp
241      zsaad        = 0.0_wp
242
243      IF( kt == nitend ) THEN
244         IF(lwp) WRITE(numout,*)
245         IF(lwp) WRITE(numout,*) 'tra_ldf_lap_adj : iso-level laplacian diffusion on ', cdtype
246         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
247         IF( .NOT. ALLOCATED( e1ur ) ) THEN
248            ! This routine may be called for both active and passive tracers.
249            ! Allocate and set saved arrays on first call only.
250            ALLOCATE( e1ur(jpi,jpj), e2vr(jpi,jpj), STAT=ierr )
251            IF( lk_mpp    )   CALL mpp_sum( ierr )
252            IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'tra_ldf_lap : unable to allocate arrays' )
253            !
254            e1ur(:,:) = e2u(:,:) / e1u(:,:)
255            e2vr(:,:) = e1v(:,:) / e2v(:,:)
256         ENDIF
257      ENDIF
258
259
260      DO jn = 1, kjpt
261         !                                                  ! =============
262         DO jk = 1, jpkm1                                   ! Vertical slab
263            !                                               ! =============
264
265            ! 2. Second derivative (divergence)
266            ! --------------------
267            DO jj = jpjm1, 2, -1
268               DO ji = fs_jpim1, fs_2, -1   ! vector opt.
269                  zbtr = 1._wp / ( e1t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk) )
270                  ! horizontal diffusive trends
271                  ztaad             = zbtr * pta_ad(ji,jj,jk,jn)
272                  ztuad(ji  ,jj  ,jk) = ztuad(ji  ,jj  ,jk) + ztaad
273                  ztuad(ji-1,jj  ,jk) = ztuad(ji-1,jj  ,jk) - ztaad
274                  ztvad(ji  ,jj  ,jk) = ztvad(ji  ,jj  ,jk) + ztaad
275                  ztvad(ji  ,jj-1,jk) = ztvad(ji  ,jj-1,jk) - ztaad
276                  !ztaad   = 0.0_wp
277                  !zsaad   = 0.0_wp
278               END DO
279            END DO
280
281            ! 1. First derivative (gradient)
282            ! -------------------
283
284            IF( ln_zps ) THEN      ! set gradient at partial step level
285               DO jj = jpjm1, 1, -1
286                  DO ji = fs_jpim1, 1, -1   ! vector opt.?
287                     ! last level
288                     iku = mbku(ji,jj)
289                     ikv = mbkv(ji,jj)
290                     IF( iku == jk ) THEN
291                        zabe1 = fsahtu(ji,jj,iku) * umask(ji,jj,iku) * e1ur(ji,jj) * fse3u(ji,jj,iku)
292                        pgu_ad(ji,jj,jn) = zabe1 * ztuad(ji,jj,jk) + pgu_ad(ji,jj,jn)
293                        ztuad(ji,jj,jk)  = 0.0_wp
294                     ENDIF
295                     IF( ikv == jk ) THEN
296                        zabe2 = fsahtv(ji,jj,ikv) * vmask(ji,jj,ikv) * e2vr(ji,jj) * fse3v(ji,jj,ikv)
297                        pgv_ad(ji,jj,jn) = zabe2 * ztvad(ji,jj,jk) + pgv_ad(ji,jj,jn)
298                        ztvad(ji,jj,jk)  = 0.0_wp
299                     ENDIF
300                  END DO
301               END DO
302            ENDIF
303            DO jj = jpjm1, 1, -1
304               DO ji = fs_jpim1, 1, -1   ! vector opt. ?
305                  zabe1 = fsahtu(ji,jj,jk) * umask(ji,jj,jk) * e1ur(ji,jj) * fse3u(ji,jj,jk)
306                  zabe2 = fsahtv(ji,jj,jk) * vmask(ji,jj,jk) * e2vr(ji,jj) * fse3v(ji,jj,jk)
307                  ptb_ad(ji  ,jj  ,jk,jn) = ptb_ad(ji  ,jj  ,jk,jn) - (zabe1 * ztuad(ji,jj,jk) + zabe2 * ztvad(ji,jj,jk))
308                  ptb_ad(ji+1,jj  ,jk,jn) = ptb_ad(ji+1,jj  ,jk,jn) +  zabe1 * ztuad(ji,jj,jk)
309                  ptb_ad(ji  ,jj+1,jk,jn) = ptb_ad(ji  ,jj+1,jk,jn) +  zabe2 * ztvad(ji,jj,jk)
310                  ztuad(ji,jj,jk)   = 0.0_wp
311                  ztvad(ji,jj,jk)   = 0.0_wp
312               END DO
313            END DO
314             !                                               ! =============
315         END DO                                              !  End of slab
316      !                                                      ! =============
317      END DO
318      !
319      CALL timing_stop('tra_ldf_lap_adj')
320      !
321      CALL wrk_dealloc(jpi,jpj,jpk, ztuad, ztvad)
322      !
323   END SUBROUTINE tra_ldf_lap_adj
324
325
326   SUBROUTINE tra_ldf_lap_adj_tst ( kumadt )
327      !!-----------------------------------------------------------------------
328      !!
329      !!                  ***  ROUTINE example_adj_tst ***
330      !!
331      !! ** Purpose : Test the adjoint routine.
332      !!
333      !! ** Method  : Verify the scalar product
334      !!
335      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
336      !!
337      !!              where  L   = tangent routine
338      !!                     L^T = adjoint routine
339      !!                     W   = diagonal matrix of scale factors
340      !!                     dx  = input perturbation (random field)
341      !!                     dy  = L dx
342      !!
343      !! History :
344      !!        ! 08-08 (A. Vidard)
345      !!-----------------------------------------------------------------------
346      !! * Modules used
347
348      !! * Arguments
349      INTEGER, INTENT(IN) :: &
350         & kumadt             ! Output unit
351
352      !! * Local declarations
353      INTEGER ::  &
354         & ji,    &        ! dummy loop indices
355         & jj,    &
356         & jk
357      INTEGER, DIMENSION(jpi,jpj) :: &
358         & iseed_2d        ! 2D seed for the random number generator
359      REAL(KIND=wp) :: &
360         & zsp1,         & ! scalar product involving the tangent routine
361         & zsp1_T,       &
362         & zsp1_S,       &
363         & zsp2,         & ! scalar product involving the adjoint routine
364         & zsp2_1,       &
365         & zsp2_2,       &
366         & zsp2_3,       &
367         & zsp2_4,       &
368         & zsp2_5,       &
369         & zsp2_6,       &
370         & zsp2_7,       &
371         & zsp2_8,       &
372         & zsp2_T,       &
373         & zsp2_S
374      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
375         & ztb_tlin ,      & ! Tangent input
376         & zsb_tlin ,      & ! Tangent input
377         & zta_tlin ,      & ! Tangent input
378         & zsa_tlin ,      & ! Tangent input
379         & zta_tlout,      & ! Tangent output
380         & zsa_tlout,      & ! Tangent output
381         & zta_adin,       & ! Adjoint input
382         & zsa_adin,       & ! Adjoint input
383         & ztb_adout ,     & ! Adjoint output
384         & zsb_adout ,     & ! Adjoint output
385         & zta_adout ,     & ! Adjoint output
386         & zsa_adout ,     & ! Adjoint output
387         & z3r               ! 3D random field
388      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
389         & zgtu_tlin ,     & ! Tangent input
390         & zgsu_tlin ,     & ! Tangent input
391         & zgtv_tlin ,     & ! Tangent input
392         & zgsv_tlin ,     & ! Tangent input
393         & zgtu_adout ,    & ! Adjoint output
394         & zgsu_adout ,    & ! Adjoint output
395         & zgtv_adout ,    & ! Adjoint output
396         & zgsv_adout ,    & ! Adjoint output
397         & z2r               ! 2D random field
398      CHARACTER(LEN=14) :: cl_name
399      ! Allocate memory
400
401      ALLOCATE( &
402         & ztb_tlin(jpi,jpj,jpk),      &
403         & zsb_tlin(jpi,jpj,jpk),      &
404         & zta_tlin(jpi,jpj,jpk),      &
405         & zsa_tlin(jpi,jpj,jpk),      &
406         & zgtu_tlin(jpi,jpj),         &
407         & zgsu_tlin(jpi,jpj),         &
408         & zgtv_tlin(jpi,jpj),         &
409         & zgsv_tlin(jpi,jpj),         &
410         & zta_tlout(jpi,jpj,jpk),     &
411         & zsa_tlout(jpi,jpj,jpk),     &
412         & zta_adin(jpi,jpj,jpk),      &
413         & zsa_adin(jpi,jpj,jpk),      &
414         & ztb_adout(jpi,jpj,jpk),     &
415         & zsb_adout(jpi,jpj,jpk),     &
416         & zta_adout(jpi,jpj,jpk),     &
417         & zsa_adout(jpi,jpj,jpk),     &
418         & zgtu_adout(jpi,jpj),        &
419         & zgsu_adout(jpi,jpj),        &
420         & zgtv_adout(jpi,jpj),        &
421         & zgsv_adout(jpi,jpj),        &
422         & z3r(jpi,jpj,jpk),           &
423         & z2r(jpi,jpj)                &
424         & )
425
426
427      !=======================================================================
428      ! 1) dx = ( tb_tl, ta_tl, sb_tl, sa_tl, gtu_tl, gtv_tl, gsu_tl, gsv_tl )
429      !    dy = ( ta_tl, sa_tl )
430      !=======================================================================
431
432      !--------------------------------------------------------------------
433      ! Reset the tangent and adjoint variables
434      !--------------------------------------------------------------------
435
436      ztb_tlin(:,:,:)  = 0.0_wp
437      zsb_tlin(:,:,:)  = 0.0_wp
438      zta_tlin(:,:,:)  = 0.0_wp
439      zsa_tlin(:,:,:)  = 0.0_wp
440      zgtu_tlin(:,:)   = 0.0_wp
441      zgsu_tlin(:,:)   = 0.0_wp
442      zgtv_tlin(:,:)   = 0.0_wp
443      zgsv_tlin(:,:)   = 0.0_wp
444      zta_tlout(:,:,:) = 0.0_wp
445      zsa_tlout(:,:,:) = 0.0_wp
446      zta_adin(:,:,:)  = 0.0_wp
447      zsa_adin(:,:,:)  = 0.0_wp
448      ztb_adout(:,:,:) = 0.0_wp
449      zsb_adout(:,:,:) = 0.0_wp
450      zta_adout(:,:,:) = 0.0_wp
451      zsa_adout(:,:,:) = 0.0_wp
452      zgtu_adout(:,:)  = 0.0_wp
453      zgsu_adout(:,:)  = 0.0_wp
454      zgtv_adout(:,:)  = 0.0_wp
455      zgsv_adout(:,:)  = 0.0_wp
456
457      tsb_tl(:,:,:,:) = 0.0_wp
458      tsa_tl(:,:,:,:) = 0.0_wp
459      gtsu_tl(:,:,:)  = 0.0_wp
460      gtsv_tl(:,:,:)  = 0.0_wp
461      tsb_ad(:,:,:,:) = 0.0_wp
462      tsa_ad(:,:,:,:) = 0.0_wp
463      gtsu_ad(:,:,:)  = 0.0_wp
464      gtsv_ad(:,:,:)  = 0.0_wp
465
466      !--------------------------------------------------------------------
467      ! Initialize the tangent input with random noise: dx
468      !--------------------------------------------------------------------
469      CALL grid_random(  z3r, 'T', 0.0_wp, stdt )
470      ztb_tlin(:,:,:) = z3r(:,:,:)
471      CALL grid_random(  z3r, 'T', 0.0_wp, stds )
472      zsb_tlin(:,:,:) = z3r(:,:,:)
473      CALL grid_random(  z3r, 'T', 0.0_wp, stdt )
474      zta_tlin(:,:,:) = z3r(:,:,:)
475      CALL grid_random(  z3r, 'T', 0.0_wp, stds )
476      zsa_tlin(:,:,:) = z3r(:,:,:)
477      CALL grid_random(  z2r, 'U', 0.0_wp, stds )
478      zgtu_tlin(:,:) = z2r(:,:)
479      CALL grid_random(  z2r, 'U', 0.0_wp, stds )
480      zgsu_tlin(:,:) = z2r(:,:)
481      CALL grid_random(  z2r, 'V', 0.0_wp, stds )
482      zgtv_tlin(:,:) = z2r(:,:)
483      CALL grid_random(  z2r, 'V', 0.0_wp, stds )
484      zgsv_tlin(:,:) = z2r(:,:)
485
486      tsb_tl(:,:,:,jp_tem) = ztb_tlin(:,:,:)
487      tsb_tl(:,:,:,jp_sal) = zsb_tlin(:,:,:)
488      tsa_tl(:,:,:,jp_tem) = zta_tlin(:,:,:)
489      tsa_tl(:,:,:,jp_sal) = zsa_tlin(:,:,:)
490      gtsu_tl(:,:,jp_tem)  = zgtu_tlin(:,:)
491      gtsu_tl(:,:,jp_sal)  = zgsu_tlin(:,:)
492      gtsv_tl(:,:,jp_tem)  = zgtv_tlin(:,:)
493      gtsv_tl(:,:,jp_sal)  = zgsv_tlin(:,:)
494
495      CALL tra_ldf_lap_tan( nit000, nit000, 'TRA', gtsu_tl, gtsv_tl, tsb_tl, tsa_tl, jpts )
496
497      zta_tlout(:,:,:) = tsa_tl(:,:,:,jp_tem)
498      zsa_tlout(:,:,:) = tsa_tl(:,:,:,jp_sal)
499
500      !--------------------------------------------------------------------
501      ! Initialize the adjoint variables: dy^* = W dy
502      !--------------------------------------------------------------------
503
504      DO jk = 1, jpk
505        DO jj = nldj, nlej
506           DO ji = nldi, nlei
507              zsa_adin(ji,jj,jk) = zsa_tlout(ji,jj,jk) &
508                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
509                 &               * tmask(ji,jj,jk) * wesp_s(jk)
510              zta_adin(ji,jj,jk) = zta_tlout(ji,jj,jk) &
511                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
512                 &               * tmask(ji,jj,jk) * wesp_t(jk)
513            END DO
514         END DO
515      END DO
516
517      !--------------------------------------------------------------------
518      ! Compute the scalar product: ( L dx )^T W dy
519      !-------------------------------------------------------------------
520
521      zsp1_T = DOT_PRODUCT( zta_tlout, zta_adin )
522      zsp1_S = DOT_PRODUCT( zsa_tlout, zsa_adin )
523      zsp1 = zsp1_T + zsp1_S
524
525      !--------------------------------------------------------------------
526      ! Call the adjoint routine: dx^* = L^T dy^*
527      !--------------------------------------------------------------------
528
529      tsa_ad(:,:,:,jp_tem) = zta_adin(:,:,:)
530      tsa_ad(:,:,:,jp_sal) = zsa_adin(:,:,:)
531
532      CALL tra_ldf_lap_adj( nit000, nit000, 'TRA', gtsu_ad, gtsv_ad, tsb_ad, tsa_ad, jpts )
533
534      ztb_adout(:,:,:) = tsb_ad(:,:,:,jp_tem)
535      zsb_adout(:,:,:) = tsb_ad(:,:,:,jp_sal)
536      zta_adout(:,:,:) = tsa_ad(:,:,:,jp_tem)
537      zsa_adout(:,:,:) = tsa_ad(:,:,:,jp_sal)
538      zgtu_adout(:,:)  = gtsu_ad(:,:,jp_tem)
539      zgsu_adout(:,:)  = gtsu_ad(:,:,jp_sal)
540      zgtv_adout(:,:)  = gtsv_ad(:,:,jp_tem)
541      zgsv_adout(:,:)  = gtsv_ad(:,:,jp_sal)
542
543      zsp2_1 = DOT_PRODUCT( ztb_tlin , ztb_adout  )
544      zsp2_2 = DOT_PRODUCT( zta_tlin , zta_adout  )
545      zsp2_3 = DOT_PRODUCT( zgtu_tlin, zgtu_adout )
546      zsp2_4 = DOT_PRODUCT( zgtv_tlin, zgtv_adout )
547      zsp2_5 = DOT_PRODUCT( zsb_tlin , zsb_adout  )
548      zsp2_6 = DOT_PRODUCT( zsa_tlin , zsa_adout  )
549      zsp2_7 = DOT_PRODUCT( zgsu_tlin, zgsu_adout )
550      zsp2_8 = DOT_PRODUCT( zgsv_tlin, zgsv_adout )
551
552      zsp2_T = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4
553      zsp2_S = zsp2_5 + zsp2_6 + zsp2_7 + zsp2_8
554      zsp2   = zsp2_T + zsp2_S
555
556      cl_name = 'tra_ldf_lap_ad'
557      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
558
559      DEALLOCATE(         &
560         & ztb_tlin,      & ! Tangent input
561         & zsb_tlin,      & ! Tangent input
562         & zta_tlin,      & ! Tangent input
563         & zsa_tlin,      & ! Tangent input
564         & zgtu_tlin,     & ! Tangent input
565         & zgsu_tlin,     & ! Tangent input
566         & zgtv_tlin,     & ! Tangent input
567         & zgsv_tlin,     & ! Tangent input
568         & zta_tlout,     & ! Tangent output
569         & zsa_tlout,     & ! Tangent output
570         & zta_adin,      & ! Adjoint input
571         & zsa_adin,      & ! Adjoint input
572         & ztb_adout,     & ! Adjoint output
573         & zsb_adout,     & ! Adjoint output
574         & zta_adout,     & ! Adjoint output
575         & zsa_adout,     & ! Adjoint output
576         & zgtu_adout,    & ! Adjoint output
577         & zgsu_adout,    & ! Adjoint output
578         & zgtv_adout,    & ! Adjoint output
579         & zgsv_adout,    & ! Adjoint output
580         & z3r,           & ! 3D random field
581         & z2r            &
582         & )
583
584
585   END SUBROUTINE tra_ldf_lap_adj_tst
586#endif
587   !!==============================================================================
588END MODULE traldf_lap_tam
Note: See TracBrowser for help on using the repository browser.