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_iso_tam.F90 in branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/TRA – NEMO

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/TRA/traldf_iso_tam.F90 @ 2587

Last change on this file since 2587 was 1885, checked in by rblod, 14 years ago

add TAM sources

File size: 44.7 KB
Line 
1MODULE traldf_iso_tam
2#if defined key_tam
3   !!======================================================================
4   !!                   ***  MODULE  traldf_iso_tam  ***
5   !! Ocean active tracers:  horizontal component of the lateral tracer mixing trend
6   !!                        Tangent and adjoint module
7   !!======================================================================
8   !! History of the direct module:
9   !!                  !  94-08  (G. Madec, M. Imbard)
10   !!                  !  97-05  (G. Madec)  split into traldf and trazdf
11   !!             8.5  !  02-08  (G. Madec)  Free form, F90
12   !!             9.0  !  05-11  (G. Madec)  merge traldf and trazdf :-)
13   !! History of the T&A module:
14   !!             9.0  !  08-12  (A. Vidard) original version
15   !!              -   !  09-01  (A. Weaver) misc. bug fixes
16   !!----------------------------------------------------------------------
17# if   defined key_ldfslp   ||   defined key_esopa
18   !!----------------------------------------------------------------------
19   !!   'key_ldfslp'               slope of the lateral diffusive direction
20   !!----------------------------------------------------------------------
21   !!----------------------------------------------------------------------
22   !!   tra_ldf_iso  : update the tracer trend with the horizontal
23   !!                  component of a iso-neutral laplacian operator
24   !!                  and with the vertical part of
25   !!                  the isopycnal or geopotential s-coord. operator
26   !!----------------------------------------------------------------------
27   USE par_kind      , ONLY: & ! Precision variables
28      & wp
29   USE par_oce       , ONLY: & ! Ocean space and time domain variables
30      & jpiglo,              &
31      & jpi,                 &
32      & jpj,                 & 
33      & jpk,                 &
34      & jpim1,               &
35      & jpjm1,               &
36      & jpkm1
37   USE oce_tam       , ONLY: & ! ocean dynamics and active tracers
38      & tb_tl,               &
39      & sb_tl,               &
40      & ta_tl,               &
41      & sa_tl,               &
42      & gtu_tl,              &
43      & gsu_tl,              &
44      & gtv_tl,              &
45      & gsv_tl,              &
46      & tb_ad,               &
47      & sb_ad,               &
48      & ta_ad,               &
49      & sa_ad,               &
50      & gtu_ad,              &
51      & gsu_ad,              &
52      & gtv_ad,              &
53      & gsv_ad
54   USE dom_oce       , ONLY: & ! ocean space and time domain
55      & e1u,                 &
56      & e2u,                 &
57      & e1v,                 &
58      & e2v,                 &
59      & e1t,                 &
60      & e2t,                 &
61#if defined key_zco
62      & e3t_0,               &
63#else
64      & e3u,                 &
65      & e3v,                 &
66      & e3t,                 &
67#endif
68      & tmask,               &
69      & umask,               &
70      & vmask,               &
71      & mbathy,              &
72      & ln_zps,              &
73      & mig,                 &
74      & mjg,                 &
75      & nldi,                &
76      & nldj,                &
77      & nlei,                &
78      & nlej
79   USE ldftra_oce    , ONLY: & ! ocean active tracers: lateral physics
80      & ahtv,                &
81      & ahtu,                &
82      & ahtw,                &
83      & ahtb0,               &
84      & aht0
85!   USE zdf_oce         ! ocean vertical physics
86   USE in_out_manager, ONLY: & ! I/O manager
87      & lwp,                 &
88      & numout,              &
89      & nitend,              &
90      & nit000
91   USE ldfslp        , ONLY: & ! iso-neutral slopes
92      & uslp,                &
93      & vslp,                &
94      & wslpi,               &
95      & wslpj
96   USE gridrandom    , ONLY: & ! Random Gaussian noise on grids
97      & grid_random
98   USE dotprodfld    , ONLY: & ! Computes dot product for 3D and 2D fields
99      & dot_product
100   USE tstool_tam    , ONLY: &
101      & prntst_adj,          & !
102      & stdt,                & ! stdev for u-velocity
103      & stds                   !           v-velocity
104   USE paresp        , ONLY: & ! Weights for an energy-type scalar product
105      & wesp_t,              &
106      & wesp_s
107
108   IMPLICIT NONE
109   PRIVATE
110
111   PUBLIC   tra_ldf_iso_tan     ! routine called by tralfd_tam.F90
112   PUBLIC   tra_ldf_iso_adj     ! routine called by traldf_tam.F90
113   PUBLIC   tra_ldf_iso_adj_tst ! routine called by traldf_tam.F90
114
115   !! * Substitutions
116#  include "domzgr_substitute.h90"
117#  include "ldftra_substitute.h90"
118#  include "vectopt_loop_substitute.h90"
119   !!----------------------------------------------------------------------
120   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
121   !!----------------------------------------------------------------------
122
123CONTAINS
124
125   SUBROUTINE tra_ldf_iso_tan( kt )
126      !!----------------------------------------------------------------------
127      !!                  ***  ROUTINE tra_ldf_iso_tan  ***
128      !!
129      !! ** Purpose of the direct routine:
130      !!      Compute the before horizontal tracer (t & s) diffusive
131      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
132      !!      add it to the general trend of tracer equation.
133      !!
134      !! ** Method of the direct routine:
135      !!      The horizontal component of the lateral diffusive trends
136      !!      is provided by a 2nd order operator rotated along neural or geopo-
137      !!      tential surfaces to which an eddy induced advection can be added
138      !!      It is computed using before fields (forward in time) and isopyc-
139      !!      nal or geopotential slopes computed in routine ldfslp.
140      !!
141      !!      1st part :  masked horizontal derivative of T & S ( di[ t ] )
142      !!      ========    with partial cell update if ln_zps=T.
143      !!
144      !!      2nd part :  horizontal fluxes of the lateral mixing operator
145      !!      ========   
146      !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ]
147      !!               - aht       e2u*uslp    dk[ mi(mk(tb)) ]
148      !!         zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ]
149      !!               - aht       e2u*vslp    dk[ mj(mk(tb)) ]
150      !!      take the horizontal divergence of the fluxes:
151      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
152      !!      Add this trend to the general trend (ta,sa):
153      !!         ta = ta + difft
154      !!
155      !!      3rd part: vertical trends of the lateral mixing operator
156      !!      ========  (excluding the vertical flux proportional to dk[t] )
157      !!      vertical fluxes associated with the rotated lateral mixing:
158      !!         zftw =-aht {  e2t*wslpi di[ mi(mk(tb)) ]
159      !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  }
160      !!      take the horizontal divergence of the fluxes:
161      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ]
162      !!      Add this trend to the general trend (ta,sa):
163      !!         ta = ta + difft
164      !!
165      !! ** Action :   Update (ta,sa) arrays with the before rotated diffusion
166      !!            trend (except the dk[ dk[.] ] term)
167      !!----------------------------------------------------------------------
168
169      INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
170      !!
171      INTEGER  ::   ji, jj, jk   ! dummy loop indices
172      INTEGER  ::   iku, ikv     ! temporary integer
173      REAL(wp) ::   zmsku, zabe1, zcof1, zcoef3, ztatl   ! temporary scalars
174      REAL(wp) ::   zmskv, zabe2, zcof2, zcoef4, zsatl   !    "         "
175      REAL(wp) ::   zcoef0, zbtr                       !    "         "
176      REAL(wp), DIMENSION(jpi,jpj)     ::   zdkttl , zdk1ttl, zftutl, zftvtl   ! 2D workspace
177      REAL(wp), DIMENSION(jpi,jpj)     ::   zdkstl , zdk1stl, zfsutl, zfsvtl   !  "     "
178      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdittl, zdjttl, ztfwtl     ! 3D workspace
179      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdistl, zdjstl, zsfwtl     !  "      "
180      !!----------------------------------------------------------------------
181
182      IF( kt == nit000 ) THEN
183         IF(lwp) WRITE(numout,*)
184         IF(lwp) WRITE(numout,*) 'tra_ldf_iso_tan : rotated laplacian diffusion operator'
185         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
186      ENDIF
187
188      !!----------------------------------------------------------------------
189      !!   I - masked horizontal derivative of T & S
190      !!----------------------------------------------------------------------
191!!bug ajout.... why?   ( 1,jpj,:) and (jpi,1,:) should be sufficient....
192      zdittl (1,:,:) = 0.e0     ;     zdittl (jpi,:,:) = 0.e0
193      zdistl (1,:,:) = 0.e0     ;     zdistl (jpi,:,:) = 0.e0
194      zdjttl (1,:,:) = 0.e0     ;     zdjttl (jpi,:,:) = 0.e0
195      zdjstl (1,:,:) = 0.e0     ;     zdjstl (jpi,:,:) = 0.e0
196!!end
197
198      ! Horizontal temperature and salinity gradient
199      DO jk = 1, jpkm1
200         DO jj = 1, jpjm1
201            DO ji = 1, fs_jpim1   ! vector opt.
202               zdittl(ji,jj,jk) = ( tb_tl(ji+1,jj  ,jk) - tb_tl(ji,jj,jk) ) * umask(ji,jj,jk)
203               zdistl(ji,jj,jk) = ( sb_tl(ji+1,jj  ,jk) - sb_tl(ji,jj,jk) ) * umask(ji,jj,jk)
204               zdjttl(ji,jj,jk) = ( tb_tl(ji  ,jj+1,jk) - tb_tl(ji,jj,jk) ) * vmask(ji,jj,jk)
205               zdjstl(ji,jj,jk) = ( sb_tl(ji  ,jj+1,jk) - sb_tl(ji,jj,jk) ) * vmask(ji,jj,jk)
206            END DO
207         END DO
208      END DO
209      IF( ln_zps ) THEN      ! partial steps correction at the last level
210         DO jj = 1, jpjm1
211            DO ji = 1, fs_jpim1   ! vector opt.
212               ! last level
213               iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
214               ikv = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
215               zdittl(ji,jj,iku) = gtu_tl(ji,jj) 
216               zdistl(ji,jj,iku) = gsu_tl(ji,jj)               
217               zdjttl(ji,jj,ikv) = gtv_tl(ji,jj) 
218               zdjstl(ji,jj,ikv) = gsv_tl(ji,jj)               
219            END DO
220         END DO
221      ENDIF
222
223      !!----------------------------------------------------------------------
224      !!   II - horizontal trend of T & S (full)
225      !!----------------------------------------------------------------------
226     
227      !                                                ! ===============
228      DO jk = 1, jpkm1                                 ! Horizontal slab
229         !                                             ! ===============
230         ! 1. Vertical tracer gradient at level jk and jk+1
231         ! ------------------------------------------------
232         ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
233
234         zdk1ttl(:,:) = ( tb_tl(:,:,jk) - tb_tl(:,:,jk+1) ) * tmask(:,:,jk+1)
235         zdk1stl(:,:) = ( sb_tl(:,:,jk) - sb_tl(:,:,jk+1) ) * tmask(:,:,jk+1)
236
237         IF( jk == 1 ) THEN
238            zdkttl(:,:) = zdk1ttl(:,:)
239            zdkstl(:,:) = zdk1stl(:,:)
240         ELSE
241            zdkttl(:,:) = ( tb_tl(:,:,jk-1) - tb_tl(:,:,jk) ) * tmask(:,:,jk)
242            zdkstl(:,:) = ( sb_tl(:,:,jk-1) - sb_tl(:,:,jk) ) * tmask(:,:,jk)
243         ENDIF
244
245         ! 2. Horizontal fluxes
246         ! --------------------
247
248         DO jj = 1 , jpjm1
249            DO ji = 1, fs_jpim1   ! vector opt.
250
251               zabe1 = ( fsahtu(ji,jj,jk) + ahtb0 ) * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
252               zabe2 = ( fsahtv(ji,jj,jk) + ahtb0 ) * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
253
254               zmsku = 1.0_wp / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
255                  &                 + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1.0_wp )
256
257               zmskv = 1.0_wp / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   &
258                  &                 + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1.0_wp )
259
260               ! *** NOTE ***  uslp() and vslp() are not linearized.
261
262               zcof1 = -fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku
263               zcof2 = -fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv
264
265               zftutl(ji,jj) = (  zabe1 * zdittl(ji,jj,jk)   &
266                  &              + zcof1 * (  zdkttl (ji+1,jj) + zdk1ttl(ji,jj)      &
267                  &                         + zdk1ttl(ji+1,jj) + zdkttl (ji,jj)  )  ) * umask(ji,jj,jk)
268               zftvtl(ji,jj) = (  zabe2 * zdjttl(ji,jj,jk)   &
269                  &              + zcof2 * (  zdkttl (ji,jj+1) + zdk1ttl(ji,jj)      &
270                  &                         + zdk1ttl(ji,jj+1) + zdkttl (ji,jj)  )  ) * vmask(ji,jj,jk)
271               zfsutl(ji,jj) = (  zabe1 * zdistl(ji,jj,jk)   &
272                  &              + zcof1 * (  zdkstl (ji+1,jj) + zdk1stl(ji,jj)      &
273                  &                         + zdk1stl(ji+1,jj) + zdkstl (ji,jj)  )  ) * umask(ji,jj,jk)
274               zfsvtl(ji,jj) = (  zabe2 * zdjstl(ji,jj,jk)   &
275                  &              + zcof2 * (  zdkstl (ji,jj+1) + zdk1stl(ji,jj)      &
276                  &                         + zdk1stl(ji,jj+1) + zdkstl (ji,jj)  )  ) * vmask(ji,jj,jk)
277
278            END DO
279         END DO
280
281
282         ! II.4 Second derivative (divergence) and add to the general trend
283         ! ----------------------------------------------------------------
284         DO jj = 2 , jpjm1
285            DO ji = fs_2, fs_jpim1   ! vector opt.
286               zbtr= 1.0_wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
287               ztatl = zbtr * (   zftutl(ji,jj) - zftutl(ji-1,jj  ) &
288                  &             + zftvtl(ji,jj) - zftvtl(ji  ,jj-1)  )
289               zsatl = zbtr * (   zfsutl(ji,jj) - zfsutl(ji-1,jj  ) &
290                  &             + zfsvtl(ji,jj) - zfsvtl(ji  ,jj-1)  )
291
292               ta_tl(ji,jj,jk) = ta_tl(ji,jj,jk) + ztatl
293               sa_tl(ji,jj,jk) = sa_tl(ji,jj,jk) + zsatl
294            END DO
295         END DO
296         !                                          ! ===============
297      END DO                                        !   End of slab 
298      !                                             ! ===============
299
300      !!----------------------------------------------------------------------
301      !!   III - vertical trend of T & S (extra diagonal terms only)
302      !!----------------------------------------------------------------------
303
304      ! Local constant initialization
305      ! -----------------------------
306      ztfwtl(1,:,:) = 0.0_wp     ;     ztfwtl(jpi,:,:) = 0.0_wp
307      zsfwtl(1,:,:) = 0.0_wp     ;     zsfwtl(jpi,:,:) = 0.0_wp
308
309
310      ! Vertical fluxes
311      ! ---------------
312
313      ! Surface and bottom vertical fluxes set to zero
314      ztfwtl(:,:, 1 ) = 0.0_wp   ;     ztfwtl(:,:,jpk) = 0.0_wp
315      zsfwtl(:,:, 1 ) = 0.0_wp   ;     zsfwtl(:,:,jpk) = 0.0_wp
316
317      ! interior (2=<jk=<jpk-1)
318      DO jk = 2, jpkm1
319         DO jj = 2, jpjm1
320            DO ji = fs_2, fs_jpim1   ! vector opt.
321               zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk)
322
323               zmsku = 1.0_wp / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      &
324                  &                  + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1.0_wp  )
325
326               zmskv = 1.0_wp / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)      &
327                  &                  + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.0_wp  )
328
329               ! *** NOTE ***  wslpi() and wslpj() are not linearized.
330
331               zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi(ji,jj,jk)
332               zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj(ji,jj,jk)
333
334               ztfwtl(ji,jj,jk) = zcoef3 * ( zdittl(ji  ,jj  ,jk-1) + zdittl(ji-1,jj  ,jk)      &
335                  &                        + zdittl(ji-1,jj  ,jk-1) + zdittl(ji  ,jj  ,jk)  )   &
336                  &             + zcoef4 * ( zdjttl(ji  ,jj  ,jk-1) + zdjttl(ji  ,jj-1,jk)      &
337                  &                        + zdjttl(ji  ,jj-1,jk-1) + zdjttl(ji  ,jj  ,jk)  )
338
339               zsfwtl(ji,jj,jk) = zcoef3 * ( zdistl(ji  ,jj  ,jk-1) + zdistl(ji-1,jj  ,jk)      &
340                  &                        + zdistl(ji-1,jj  ,jk-1) + zdistl(ji  ,jj  ,jk)  )   &
341                  &             + zcoef4 * ( zdjstl(ji  ,jj  ,jk-1) + zdjstl(ji  ,jj-1,jk)      &
342                  &                        + zdjstl(ji  ,jj-1,jk-1) + zdjstl(ji  ,jj  ,jk)  )
343            END DO
344         END DO
345      END DO
346
347
348      ! I.5 Divergence of vertical fluxes added to the general tracer trend
349      ! -------------------------------------------------------------------
350
351      DO jk = 1, jpkm1
352         DO jj = 2, jpjm1
353            DO ji = fs_2, fs_jpim1   ! vector opt.
354               zbtr =  1.0_wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
355
356               ztatl = ( ztfwtl(ji,jj,jk) - ztfwtl(ji,jj,jk+1) ) * zbtr
357               zsatl = ( zsfwtl(ji,jj,jk) - zsfwtl(ji,jj,jk+1) ) * zbtr
358
359               ta_tl(ji,jj,jk) = ta_tl(ji,jj,jk) + ztatl
360               sa_tl(ji,jj,jk) = sa_tl(ji,jj,jk) + zsatl
361            END DO
362         END DO
363      END DO
364      !
365   END SUBROUTINE tra_ldf_iso_tan
366
367   SUBROUTINE tra_ldf_iso_adj( kt )
368      !!----------------------------------------------------------------------
369      !!                  ***  ROUTINE tra_ldf_iso_adj  ***
370      !!
371      !! ** Purpose of the direct routine:
372      !!      Compute the before horizontal tracer (t & s) diffusive
373      !!      trend for a laplacian tensor (ezxcept the dz[ dz[.] ] term) and
374      !!      add it to the general trend of tracer equation.
375      !!
376      !! ** Method of the direct routine:
377      !!      The horizontal component of the lateral diffusive trends
378      !!      is provided by a 2nd order operator rotated along neural or geopo-
379      !!      tential surfaces to which an eddy induced advection can be added
380      !!      It is computed using before fields (forward in time) and isopyc-
381      !!      nal or geopotential slopes computed in routine ldfslp.
382      !!
383      !!      1st part :  masked horizontal derivative of T & S ( di[ t ] )
384      !!      ========    with partial cell update if ln_zps=T.
385      !!
386      !!      2nd part :  horizontal fluxes of the lateral mixing operator
387      !!      ========   
388      !!         zftu = (aht+ahtb0) e2u*e3u/e1u di[ tb ]
389      !!               - aht       e2u*uslp    dk[ mi(mk(tb)) ]
390      !!         zftv = (aht+ahtb0) e1v*e3v/e2v dj[ tb ]
391      !!               - aht       e2u*vslp    dk[ mj(mk(tb)) ]
392      !!      take the horizontal divergence of the fluxes:
393      !!         difft = 1/(e1t*e2t*e3t) {  di-1[ zftu ] +  dj-1[ zftv ]  }
394      !!      Add this trend to the general trend (ta,sa):
395      !!         ta = ta + difft
396      !!
397      !!      3rd part: vertical trends of the lateral mixing operator
398      !!      ========  (excluding the vertical flux proportional to dk[t] )
399      !!      vertical fluxes associated with the rotated lateral mixing:
400      !!         zftw =-aht {  e2t*wslpi di[ mi(mk(tb)) ]
401      !!                     + e1t*wslpj dj[ mj(mk(tb)) ]  }
402      !!      take the horizontal divergence of the fluxes:
403      !!         difft = 1/(e1t*e2t*e3t) dk[ zftw ]
404      !!      Add this trend to the general trend (ta,sa):
405      !!         ta = ta + difft
406      !!
407      !! ** Action :   Update (ta,sa) arrays with the before rotated diffusion
408      !!            trend (except the dk[ dk[.] ] term)
409      !!----------------------------------------------------------------------
410
411      INTEGER, INTENT( in ) ::   kt    ! ocean time-step index
412      !!
413      INTEGER  ::   ji, jj, jk   ! dummy loop indices
414      INTEGER  ::   iku, ikv     ! temporary integer
415      REAL(wp) ::   zmsku, zabe1, zcof1, zcoef3, ztaad   ! temporary scalars
416      REAL(wp) ::   zmskv, zabe2, zcof2, zcoef4, zsaad   !    "         "
417      REAL(wp) ::   ztf3, ztf4, zsf3, zsf4               !
418      REAL(wp) ::   zcoef0, zbtr                       !    "         "
419      REAL(wp), DIMENSION(jpi,jpj)     ::   zdktad , zdk1tad, zftuad, zftvad   ! 2D workspace
420      REAL(wp), DIMENSION(jpi,jpj)     ::   zdksad , zdk1sad, zfsuad, zfsvad   !  "     "
421      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zditad, zdjtad, ztfwad     ! 3D workspace
422      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdisad, zdjsad, zsfwad     !  "      "
423      !!----------------------------------------------------------------------
424
425      zditad(:,:,:) = 0.0_wp ;  zdjtad(:,:,:) = 0.0_wp ;  ztfwad(:,:,:) = 0.0_wp
426      zdisad(:,:,:) = 0.0_wp ;  zdjsad(:,:,:) = 0.0_wp ;  zsfwad(:,:,:) = 0.0_wp
427
428      zdktad(:,:) = 0.0_wp ;  zdk1tad(:,:) = 0.0_wp   
429      zdksad(:,:) = 0.0_wp ;  zdk1sad(:,:) = 0.0_wp 
430
431      zftuad(:,:) = 0.0_wp ;  zftvad (:,:) = 0.0_wp 
432      zfsuad(:,:) = 0.0_wp ;  zfsvad (:,:) = 0.0_wp
433
434      IF( kt == nitend ) THEN
435         IF(lwp) WRITE(numout,*)
436         IF(lwp) WRITE(numout,*) 'tra_ldf_iso_adj : rotated laplacian diffusion operator'
437         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
438      ENDIF
439
440      !!----------------------------------------------------------------------
441      !!   III - vertical trend of T & S (extra diagonal terms only)
442      !!----------------------------------------------------------------------
443      ! I.5 Divergence of vertical fluxes added to the general tracer trend
444      ! -------------------------------------------------------------------
445
446      DO jk = jpkm1, 1, -1
447         DO jj = jpjm1, 2, -1
448            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
449               zbtr =  1.0_wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
450               ztaad = ta_ad(ji,jj,jk) * zbtr
451               zsaad = sa_ad(ji,jj,jk) * zbtr
452
453               ztfwad(ji,jj,jk  ) = ztfwad(ji,jj,jk  ) + ztaad
454               ztfwad(ji,jj,jk+1) = ztfwad(ji,jj,jk+1) - ztaad
455               zsfwad(ji,jj,jk  ) = zsfwad(ji,jj,jk  ) + zsaad
456               zsfwad(ji,jj,jk+1) = zsfwad(ji,jj,jk+1) - zsaad
457            END DO
458         END DO
459      END DO
460      ! interior (2=<jk=<jpk-1)
461      DO jk = jpkm1, 2, -1
462         DO jj = jpjm1, 2, -1
463            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
464               zcoef0 = - fsahtw(ji,jj,jk) * tmask(ji,jj,jk)
465
466               zmsku = 1.0_wp / MAX(   umask(ji  ,jj,jk-1) + umask(ji-1,jj,jk)      &
467                  &                  + umask(ji-1,jj,jk-1) + umask(ji  ,jj,jk), 1.0_wp  )
468
469               zmskv = 1.0_wp / MAX(   vmask(ji,jj  ,jk-1) + vmask(ji,jj-1,jk)      &
470                  &                  + vmask(ji,jj-1,jk-1) + vmask(ji,jj  ,jk), 1.0_wp  )
471
472               ! *** NOTE ***  wslpi() and wslpj() are not linearized.
473
474               zcoef3 = zcoef0 * e2t(ji,jj) * zmsku * wslpi(ji,jj,jk)
475               zcoef4 = zcoef0 * e1t(ji,jj) * zmskv * wslpj(ji,jj,jk)
476
477               ztf3 = ztfwad(ji,jj,jk) * zcoef3
478               ztf4 = ztfwad(ji,jj,jk) * zcoef4
479               zsf3 = zsfwad(ji,jj,jk) * zcoef3
480               zsf4 = zsfwad(ji,jj,jk) * zcoef4
481
482               zditad(ji  ,jj  ,jk-1) = zditad(ji  ,jj  ,jk-1) + ztf3
483               zditad(ji-1,jj  ,jk  ) = zditad(ji-1,jj  ,jk  ) + ztf3
484               zditad(ji-1,jj  ,jk-1) = zditad(ji-1,jj  ,jk-1) + ztf3
485               zditad(ji  ,jj  ,jk  ) = zditad(ji  ,jj  ,jk  ) + ztf3
486
487               zdjtad(ji  ,jj  ,jk-1) = zdjtad(ji  ,jj  ,jk-1) + ztf4
488               zdjtad(ji  ,jj-1,jk  ) = zdjtad(ji  ,jj-1,jk  ) + ztf4
489               zdjtad(ji  ,jj-1,jk-1) = zdjtad(ji  ,jj-1,jk-1) + ztf4
490               zdjtad(ji  ,jj  ,jk  ) = zdjtad(ji  ,jj  ,jk  ) + ztf4
491
492               zdisad(ji  ,jj  ,jk-1) = zdisad(ji  ,jj  ,jk-1) + zsf3
493               zdisad(ji-1,jj  ,jk  ) = zdisad(ji-1,jj  ,jk  ) + zsf3
494               zdisad(ji-1,jj  ,jk-1) = zdisad(ji-1,jj  ,jk-1) + zsf3
495               zdisad(ji  ,jj  ,jk  ) = zdisad(ji  ,jj  ,jk  ) + zsf3
496
497               zdjsad(ji  ,jj  ,jk-1) = zdjsad(ji  ,jj  ,jk-1) + zsf4
498               zdjsad(ji  ,jj-1,jk  ) = zdjsad(ji  ,jj-1,jk  ) + zsf4
499               zdjsad(ji  ,jj-1,jk-1) = zdjsad(ji  ,jj-1,jk-1) + zsf4
500               zdjsad(ji  ,jj  ,jk  ) = zdjsad(ji  ,jj  ,jk  ) + zsf4
501
502               ztfwad(ji,jj,jk) = 0.0_wp
503               zsfwad(ji,jj,jk) = 0.0_wp
504            END DO
505         END DO
506      END DO
507 
508      ! Local constant initialization
509      ! -----------------------------
510      ztfwad(1,:,:) = 0.0_wp     ;     ztfwad(jpi,:,:) = 0.0_wp
511      zsfwad(1,:,:) = 0.0_wp     ;     zsfwad(jpi,:,:) = 0.0_wp
512
513      ! Vertical fluxes
514      ! ---------------
515
516      ! Surface and bottom vertical fluxes set to zero
517      ztfwad(:,:, 1 ) = 0.0_wp   ;     ztfwad(:,:,jpk) = 0.0_wp
518      zsfwad(:,:, 1 ) = 0.0_wp   ;     zsfwad(:,:,jpk) = 0.0_wp
519
520      !!----------------------------------------------------------------------
521      !!   II - horizontal trend of T & S (full)
522      !!----------------------------------------------------------------------
523     
524      !                                                ! ===============
525      DO jk = jpkm1, 1, -1                             ! Horizontal slab
526         !                                             ! ===============
527         ! II.4 Second derivative (divergence) and add to the general trend
528         ! ----------------------------------------------------------------
529         DO jj = jpjm1, 2, -1
530            DO ji = fs_jpim1, fs_2, -1   ! vector opt.
531
532               zbtr= 1.0_wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
533               ztaad = ta_ad(ji,jj,jk) * zbtr
534               zsaad = sa_ad(ji,jj,jk) * zbtr
535
536               zftuad(ji  ,jj  ) = zftuad(ji  ,jj  ) + ztaad
537               zftuad(ji-1,jj  ) = zftuad(ji-1,jj  ) - ztaad
538               zftvad(ji  ,jj  ) = zftvad(ji  ,jj  ) + ztaad
539               zftvad(ji  ,jj-1) = zftvad(ji  ,jj-1) - ztaad
540
541               zfsuad(ji  ,jj  ) = zfsuad(ji  ,jj  ) + zsaad
542               zfsuad(ji-1,jj  ) = zfsuad(ji-1,jj  ) - zsaad
543               zfsvad(ji  ,jj  ) = zfsvad(ji  ,jj  ) + zsaad
544               zfsvad(ji  ,jj-1) = zfsvad(ji  ,jj-1) - zsaad
545
546            END DO
547         END DO
548
549         ! 2. Horizontal fluxes
550         ! --------------------
551         DO jj = jpjm1, 1, -1
552            DO ji = fs_jpim1, 1, -1   ! vector opt.
553               zabe1 = umask(ji,jj,jk) * ( fsahtu(ji,jj,jk) + ahtb0 ) &
554                  &                    * e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj)
555               zabe2 = vmask(ji,jj,jk) * ( fsahtv(ji,jj,jk) + ahtb0 ) &
556                  &                    * e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj)
557
558               zmsku = 1.0_wp / MAX(  tmask(ji+1,jj,jk  ) + tmask(ji,jj,jk+1)   &
559                  &                 + tmask(ji+1,jj,jk+1) + tmask(ji,jj,jk  ), 1.0_wp )
560
561               zmskv = 1.0_wp / MAX(  tmask(ji,jj+1,jk  ) + tmask(ji,jj,jk+1)   &
562                  &                 + tmask(ji,jj+1,jk+1) + tmask(ji,jj,jk  ), 1.0_wp )
563
564               ! *** NOTE ***  uslp() and vslp() are not linearized.
565
566               zcof1 = -fsahtu(ji,jj,jk) * e2u(ji,jj) * uslp(ji,jj,jk) * zmsku * umask(ji,jj,jk)
567               zcof2 = -fsahtv(ji,jj,jk) * e1v(ji,jj) * vslp(ji,jj,jk) * zmskv * vmask(ji,jj,jk)
568
569               zditad(ji,jj,jk) = zditad(ji,jj,jk) + zftuad(ji,jj) * zabe1
570
571               zdktad (ji+1,jj) = zdktad (ji+1,jj) + zftuad(ji,jj) * zcof1
572               zdktad (ji  ,jj) = zdktad (ji  ,jj) + zftuad(ji,jj) * zcof1
573               zdk1tad(ji  ,jj) = zdk1tad(ji  ,jj) + zftuad(ji,jj) * zcof1
574               zdk1tad(ji+1,jj) = zdk1tad(ji+1,jj) + zftuad(ji,jj) * zcof1
575               zftuad (ji  ,jj) = 0.0_wp
576               !
577               zdjtad(ji,jj,jk) = zdjtad(ji,jj,jk) + zftvad(ji,jj) * zabe2
578
579               zdktad (ji,jj+1) = zdktad (ji,jj+1) + zftvad(ji,jj) * zcof2
580               zdktad (ji,jj  ) = zdktad (ji,jj  ) + zftvad(ji,jj) * zcof2
581               zdk1tad(ji,jj  ) = zdk1tad(ji,jj  ) + zftvad(ji,jj) * zcof2
582               zdk1tad(ji,jj+1) = zdk1tad(ji,jj+1) + zftvad(ji,jj) * zcof2
583               zftvad (ji,jj  ) = 0.0_wp
584               !
585               zdisad(ji,jj,jk) = zdisad(ji,jj,jk) + zfsuad(ji,jj) * zabe1
586
587               zdksad (ji+1,jj) = zdksad (ji+1,jj) + zfsuad(ji,jj) * zcof1 
588               zdksad (ji  ,jj) = zdksad (ji  ,jj) + zfsuad(ji,jj) * zcof1 
589               zdk1sad(ji  ,jj) = zdk1sad(ji  ,jj) + zfsuad(ji,jj) * zcof1
590               zdk1sad(ji+1,jj) = zdk1sad(ji+1,jj) + zfsuad(ji,jj) * zcof1 
591               zfsuad (ji  ,jj) = 0.0_wp
592               !
593               zdjsad(ji,jj,jk) = zdjsad(ji,jj,jk) + zfsvad(ji,jj) * zabe2
594
595               zdksad (ji,jj+1) = zdksad (ji,jj+1) + zfsvad(ji,jj) * zcof2
596               zdksad (ji,jj  ) = zdksad (ji,jj  ) + zfsvad(ji,jj) * zcof2
597               zdk1sad(ji,jj  ) = zdk1sad(ji,jj  ) + zfsvad(ji,jj) * zcof2
598               zdk1sad(ji,jj+1) = zdk1sad(ji,jj+1) + zfsvad(ji,jj) * zcof2
599               zfsvad (ji,jj  ) = 0.0_wp
600
601            END DO
602         END DO
603
604         ! 1. Vertical tracer gradient at level jk and jk+1
605         ! ------------------------------------------------
606         ! surface boundary condition: zdkt(jk=1)=zdkt(jk=2)
607
608         IF( jk == 1 ) THEN
609
610            zdk1tad(:,:) = zdk1tad(:,:) + zdktad(:,:)
611            zdk1sad(:,:) = zdk1sad(:,:) + zdksad(:,:)
612
613            zdktad(:,:)  = 0.0_wp
614            zdksad(:,:)  = 0.0_wp
615
616         ELSE
617
618            tb_ad(:,:,jk-1) = tb_ad(:,:,jk-1) + zdktad(:,:) * tmask(:,:,jk)
619            tb_ad(:,:,jk  ) = tb_ad(:,:,jk  ) - zdktad(:,:) * tmask(:,:,jk)
620
621            sb_ad(:,:,jk-1) = sb_ad(:,:,jk-1) + zdksad(:,:) * tmask(:,:,jk)
622            sb_ad(:,:,jk  ) = sb_ad(:,:,jk  ) - zdksad(:,:) * tmask(:,:,jk)
623
624            zdktad(:,:) = 0.0_wp
625            zdksad(:,:) = 0.0_wp
626
627         ENDIF
628
629         tb_ad(:,:,jk  ) = tb_ad(:,:,jk  ) + zdk1tad(:,:) * tmask(:,:,jk+1)
630         tb_ad(:,:,jk+1) = tb_ad(:,:,jk+1) - zdk1tad(:,:) * tmask(:,:,jk+1)
631         
632         sb_ad(:,:,jk  ) = sb_ad(:,:,jk  ) + zdk1sad(:,:) * tmask(:,:,jk+1)
633         sb_ad(:,:,jk+1) = sb_ad(:,:,jk+1) - zdk1sad(:,:) * tmask(:,:,jk+1)
634         
635         zdk1tad(:,:) = 0.0_wp
636         zdk1sad(:,:) = 0.0_wp
637         !                                          ! ===============
638      END DO                                        !   End of slab 
639      !                                             ! ===============
640      !!----------------------------------------------------------------------
641      !!   I - masked horizontal derivative of T & S
642      !!----------------------------------------------------------------------
643      IF( ln_zps ) THEN      ! partial steps correction at the last level
644         DO jj = jpjm1, 1, -1
645            DO ji = fs_jpim1, 1, -1   ! vector opt.
646
647               ! last level
648               iku = MIN( mbathy(ji,jj), mbathy(ji+1,jj  ) ) - 1
649               ikv = MIN( mbathy(ji,jj), mbathy(ji  ,jj+1) ) - 1
650
651               gtu_ad(ji,jj) = gtu_ad(ji,jj) + zditad(ji,jj,iku)
652               gtv_ad(ji,jj) = gtv_ad(ji,jj) + zdjtad(ji,jj,ikv)
653
654               gsu_ad(ji,jj) = gsu_ad(ji,jj) + zdisad(ji,jj,iku)
655               gsv_ad(ji,jj) = gsv_ad(ji,jj) + zdjsad(ji,jj,ikv)
656
657               zditad(ji,jj,iku) = 0.0_wp
658               zdjtad(ji,jj,ikv) = 0.0_wp
659
660               zdisad(ji,jj,iku) = 0.0_wp
661               zdjsad(ji,jj,ikv) = 0.0_wp
662 
663           END DO
664         END DO
665      ENDIF
666
667      ! Horizontal temperature and salinity gradient
668      DO jk = jpkm1, 1, -1
669         DO jj = jpjm1, 1, -1
670            DO ji = fs_jpim1, 1, -1   ! vector opt.
671
672               zditad(ji,jj,jk) = zditad(ji,jj,jk) * umask(ji,jj,jk)
673               zdjtad(ji,jj,jk) = zdjtad(ji,jj,jk) * vmask(ji,jj,jk)
674               zdisad(ji,jj,jk) = zdisad(ji,jj,jk) * umask(ji,jj,jk)
675               zdjsad(ji,jj,jk) = zdjsad(ji,jj,jk) * vmask(ji,jj,jk)
676
677               tb_ad(ji+1,jj  ,jk) = tb_ad(ji+1,jj  ,jk) + zditad(ji,jj,jk)
678               tb_ad(ji  ,jj  ,jk) = tb_ad(ji  ,jj  ,jk) - zditad(ji,jj,jk)
679               tb_ad(ji  ,jj+1,jk) = tb_ad(ji  ,jj+1,jk) + zdjtad(ji,jj,jk)
680               tb_ad(ji  ,jj  ,jk) = tb_ad(ji  ,jj  ,jk) - zdjtad(ji,jj,jk)
681
682               sb_ad(ji+1,jj  ,jk) = sb_ad(ji+1,jj  ,jk) + zdisad(ji,jj,jk)
683               sb_ad(ji  ,jj  ,jk) = sb_ad(ji  ,jj  ,jk) - zdisad(ji,jj,jk)
684               sb_ad(ji  ,jj+1,jk) = sb_ad(ji  ,jj+1,jk) + zdjsad(ji,jj,jk)
685               sb_ad(ji  ,jj  ,jk) = sb_ad(ji  ,jj  ,jk) - zdjsad(ji,jj,jk)
686
687            END DO
688         END DO
689      END DO
690      !
691   END SUBROUTINE tra_ldf_iso_adj
692
693   SUBROUTINE tra_ldf_iso_adj_tst ( kumadt )
694      !!-----------------------------------------------------------------------
695      !!
696      !!                  ***  ROUTINE example_adj_tst ***
697      !!
698      !! ** Purpose : Test the adjoint routine.
699      !!
700      !! ** Method  : Verify the scalar product
701      !!           
702      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
703      !!
704      !!              where  L   = tangent routine
705      !!                     L^T = adjoint routine
706      !!                     W   = diagonal matrix of scale factors
707      !!                     dx  = input perturbation (random field)
708      !!                     dy  = L dx
709      !!                   
710      !! History :
711      !!        ! 08-08 (A. Vidard)
712      !!-----------------------------------------------------------------------
713      !! * Modules used
714
715      !! * Arguments
716      INTEGER, INTENT(IN) :: &
717         & kumadt             ! Output unit
718 
719      !! * Local declarations
720      INTEGER ::  &
721         & ji,    &        ! dummy loop indices
722         & jj,    &       
723         & jk     
724      INTEGER, DIMENSION(jpi,jpj) :: &
725         & iseed_2d        ! 2D seed for the random number generator
726      REAL(KIND=wp) :: &
727         & zsp1,         & ! scalar product involving the tangent routine
728         & zsp1_T,       &
729         & zsp1_S,       &
730         & zsp2,         & ! scalar product involving the adjoint routine
731         & zsp2_1,       &
732         & zsp2_2,       &
733         & zsp2_3,       &
734         & zsp2_4,       &
735         & zsp2_5,       &
736         & zsp2_6,       &
737         & zsp2_7,       &
738         & zsp2_8,       &
739         & zsp2_T,       &
740         & zsp2_S
741      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
742         & ztb_tlin ,      & ! Tangent input
743         & zsb_tlin ,      & ! Tangent input
744         & zta_tlin ,      & ! Tangent input
745         & zsa_tlin ,      & ! Tangent input
746         & zta_tlout,      & ! Tangent output
747         & zsa_tlout,      & ! Tangent output
748         & zta_adin,       & ! Adjoint input
749         & zsa_adin,       & ! Adjoint input
750         & ztb_adout ,     & ! Adjoint output
751         & zsb_adout ,     & ! Adjoint output
752         & zta_adout ,     & ! Adjoint output
753         & zsa_adout ,     & ! Adjoint output
754         & z3r               ! 3D random field
755      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
756         & zgtu_tlin ,     & ! Tangent input
757         & zgsu_tlin ,     & ! Tangent input
758         & zgtv_tlin ,     & ! Tangent input
759         & zgsv_tlin ,     & ! Tangent input
760         & zgtu_adout ,    & ! Adjoint output
761         & zgsu_adout ,    & ! Adjoint output
762         & zgtv_adout ,    & ! Adjoint output
763         & zgsv_adout ,    & ! Adjoint output
764         & z2r               ! 2D random field
765      CHARACTER(LEN=14) :: cl_name
766      ! Allocate memory
767
768      ALLOCATE( &
769         & ztb_tlin(jpi,jpj,jpk),      & 
770         & zsb_tlin(jpi,jpj,jpk),      & 
771         & zta_tlin(jpi,jpj,jpk),      & 
772         & zsa_tlin(jpi,jpj,jpk),      & 
773         & zgtu_tlin(jpi,jpj),         & 
774         & zgsu_tlin(jpi,jpj),         & 
775         & zgtv_tlin(jpi,jpj),         & 
776         & zgsv_tlin(jpi,jpj),         & 
777         & zta_tlout(jpi,jpj,jpk),     & 
778         & zsa_tlout(jpi,jpj,jpk),     & 
779         & zta_adin(jpi,jpj,jpk),      & 
780         & zsa_adin(jpi,jpj,jpk),      & 
781         & ztb_adout(jpi,jpj,jpk),     & 
782         & zsb_adout(jpi,jpj,jpk),     & 
783         & zta_adout(jpi,jpj,jpk),     & 
784         & zsa_adout(jpi,jpj,jpk),     & 
785         & zgtu_adout(jpi,jpj),        & 
786         & zgsu_adout(jpi,jpj),        & 
787         & zgtv_adout(jpi,jpj),        & 
788         & zgsv_adout(jpi,jpj),        & 
789         & z3r(jpi,jpj,jpk),           &
790         & z2r(jpi,jpj)                &
791         & )
792
793      ! Initialize the reference state
794      uslp (:,:,:) = 2.0_wp
795      vslp (:,:,:) = 3.0_wp
796      wslpi(:,:,:) = 4.0_wp
797      wslpj(:,:,:) = 5.0_wp
798
799      !=======================================================================
800      ! 1) dx = ( tb_tl, ta_tl, sb_tl, sa_tl, gtu_tl, gtv_tl, gsu_tl, gsv_tl )
801      !    dy = ( ta_tl, sa_tl )
802      !=======================================================================
803
804      !--------------------------------------------------------------------
805      ! Reset the tangent and adjoint variables
806      !--------------------------------------------------------------------
807
808      ztb_tlin(:,:,:)  = 0.0_wp
809      zsb_tlin(:,:,:)  = 0.0_wp
810      zta_tlin(:,:,:)  = 0.0_wp
811      zsa_tlin(:,:,:)  = 0.0_wp
812      zgtu_tlin(:,:)   = 0.0_wp
813      zgsu_tlin(:,:)   = 0.0_wp
814      zgtv_tlin(:,:)   = 0.0_wp
815      zgsv_tlin(:,:)   = 0.0_wp
816      zta_tlout(:,:,:) = 0.0_wp
817      zsa_tlout(:,:,:) = 0.0_wp
818      zta_adin(:,:,:)  = 0.0_wp
819      zsa_adin(:,:,:)  = 0.0_wp
820      ztb_adout(:,:,:) = 0.0_wp
821      zsb_adout(:,:,:) = 0.0_wp
822      zta_adout(:,:,:) = 0.0_wp
823      zsa_adout(:,:,:) = 0.0_wp
824      zgtu_adout(:,:)  = 0.0_wp
825      zgsu_adout(:,:)  = 0.0_wp
826      zgtv_adout(:,:)  = 0.0_wp
827      zgsv_adout(:,:)  = 0.0_wp
828
829      tb_tl(:,:,:) = 0.0_wp
830      sb_tl(:,:,:) = 0.0_wp
831      ta_tl(:,:,:) = 0.0_wp
832      sa_tl(:,:,:) = 0.0_wp
833      gtu_tl(:,:)  = 0.0_wp
834      gsu_tl(:,:)  = 0.0_wp
835      gtv_tl(:,:)  = 0.0_wp
836      gsv_tl(:,:)  = 0.0_wp
837      tb_ad(:,:,:) = 0.0_wp
838      sb_ad(:,:,:) = 0.0_wp
839      ta_ad(:,:,:) = 0.0_wp
840      sa_ad(:,:,:) = 0.0_wp
841      gtu_ad(:,:)  = 0.0_wp
842      gsu_ad(:,:)  = 0.0_wp
843      gtv_ad(:,:)  = 0.0_wp
844      gsv_ad(:,:)  = 0.0_wp
845
846      !--------------------------------------------------------------------
847      ! Initialize the tangent input with random noise: dx
848      !--------------------------------------------------------------------
849
850      DO jj = 1, jpj
851         DO ji = 1, jpi
852            iseed_2d(ji,jj) = - ( 596035 + &
853               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
854         END DO
855      END DO
856      CALL grid_random( iseed_2d, z3r, 'T', 0.0_wp, stdt )
857      DO jk = 1, jpk
858        DO jj = nldj, nlej
859           DO ji = nldi, nlei
860              ztb_tlin(ji,jj,jk) = z3r(ji,jj,jk) 
861            END DO
862         END DO
863      END DO
864
865      DO jj = 1, jpj
866         DO ji = 1, jpi
867            iseed_2d(ji,jj) = - ( 727391 + &
868               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
869         END DO
870      END DO
871      CALL grid_random( iseed_2d, z3r, 'T', 0.0_wp, stds ) 
872      DO jk = 1, jpk
873        DO jj = nldj, nlej
874           DO ji = nldi, nlei
875              zsb_tlin(ji,jj,jk) = z3r(ji,jj,jk) 
876            END DO
877         END DO
878      END DO
879
880      DO jj = 1, jpj
881         DO ji = 1, jpi
882            iseed_2d(ji,jj) = - ( 249741 + &
883               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
884         END DO
885      END DO
886      CALL grid_random( iseed_2d, z3r, 'T', 0.0_wp, stdt )
887      DO jk = 1, jpk
888        DO jj = nldj, nlej
889           DO ji = nldi, nlei
890              zta_tlin(ji,jj,jk) = z3r(ji,jj,jk) 
891            END DO
892         END DO
893      END DO
894
895      DO jj = 1, jpj
896         DO ji = 1, jpi
897            iseed_2d(ji,jj) = - ( 182029 + &
898               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
899         END DO
900      END DO
901      CALL grid_random( iseed_2d, z3r, 'T', 0.0_wp, stds )
902      DO jk = 1, jpk
903        DO jj = nldj, nlej
904           DO ji = nldi, nlei
905              zsa_tlin(ji,jj,jk) = z3r(ji,jj,jk) 
906            END DO
907         END DO
908      END DO
909
910      DO jj = 1, jpj
911         DO ji = 1, jpi
912            iseed_2d(ji,jj) = - ( 596035 + &
913               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
914         END DO
915      END DO
916      CALL grid_random( iseed_2d, z2r, 'U', 0.0_wp, stds )
917      DO jj = nldj, nlej
918         DO ji = nldi, nlei
919            zgtu_tlin(ji,jj) = z2r(ji,jj) 
920         END DO
921      END DO
922
923      DO jj = 1, jpj
924         DO ji = 1, jpi
925            iseed_2d(ji,jj) = - ( 727391 + &
926               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
927         END DO
928      END DO
929      CALL grid_random( iseed_2d, z2r, 'U', 0.0_wp, stds ) 
930      DO jj = nldj, nlej
931        DO ji = nldi, nlei
932           zgsu_tlin(ji,jj) = z2r(ji,jj) 
933        END DO
934      END DO
935
936      DO jj = 1, jpj
937         DO ji = 1, jpi
938            iseed_2d(ji,jj) = - ( 249741 + &
939               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
940         END DO
941      END DO
942      CALL grid_random( iseed_2d, z2r, 'V', 0.0_wp, stds )
943      DO jj = nldj, nlej
944        DO ji = nldi, nlei
945           zgtv_tlin(ji,jj) = z2r(ji,jj) 
946        END DO
947      END DO
948
949      DO jj = 1, jpj
950         DO ji = 1, jpi
951            iseed_2d(ji,jj) = - ( 182029 + &
952               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
953         END DO
954      END DO
955      CALL grid_random( iseed_2d, z2r, 'V', 0.0_wp, stds ) 
956      DO jj = nldj, nlej
957        DO ji = nldi, nlei
958           zgsv_tlin(ji,jj) = z2r(ji,jj) 
959        END DO
960      END DO
961
962      tb_tl(:,:,:) = ztb_tlin(:,:,:) 
963      sb_tl(:,:,:) = zsb_tlin(:,:,:) 
964      ta_tl(:,:,:) = zta_tlin(:,:,:) 
965      sa_tl(:,:,:) = zsa_tlin(:,:,:) 
966      gtu_tl(:,:)  = zgtu_tlin(:,:)
967      gsu_tl(:,:)  = zgsu_tlin(:,:)
968      gtv_tl(:,:)  = zgtv_tlin(:,:)
969      gsv_tl(:,:)  = zgsv_tlin(:,:)
970
971      CALL tra_ldf_iso_tan( nit000 ) 
972     
973      zta_tlout(:,:,:) = ta_tl(:,:,:)
974      zsa_tlout(:,:,:) = sa_tl(:,:,:)
975
976      !--------------------------------------------------------------------
977      ! Initialize the adjoint variables: dy^* = W dy
978      !--------------------------------------------------------------------
979
980      DO jk = 1, jpk
981        DO jj = nldj, nlej
982           DO ji = nldi, nlei
983              zsa_adin(ji,jj,jk) = zsa_tlout(ji,jj,jk) &
984                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
985                 &               * tmask(ji,jj,jk) * wesp_s(jk)
986              zta_adin(ji,jj,jk) = zta_tlout(ji,jj,jk) &
987                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
988                 &               * tmask(ji,jj,jk) * wesp_t(jk)
989            END DO
990         END DO
991      END DO
992
993      !--------------------------------------------------------------------
994      ! Compute the scalar product: ( L dx )^T W dy
995      !--------------------------------------------------------------------
996
997      zsp1_T = DOT_PRODUCT( zta_tlout, zta_adin )
998      zsp1_S = DOT_PRODUCT( zsa_tlout, zsa_adin )
999      zsp1 = zsp1_T + zsp1_S
1000
1001      !--------------------------------------------------------------------
1002      ! Call the adjoint routine: dx^* = L^T dy^*
1003      !--------------------------------------------------------------------
1004
1005      ta_ad(:,:,:) = zta_adin(:,:,:)
1006      sa_ad(:,:,:) = zsa_adin(:,:,:)
1007
1008      CALL tra_ldf_iso_adj( nit000 ) 
1009
1010      ztb_adout(:,:,:) = tb_ad(:,:,:)
1011      zsb_adout(:,:,:) = sb_ad(:,:,:)
1012      zta_adout(:,:,:) = ta_ad(:,:,:)
1013      zsa_adout(:,:,:) = sa_ad(:,:,:)
1014      zgtu_adout(:,:)  = gtu_ad(:,:)
1015      zgsu_adout(:,:)  = gsu_ad(:,:)
1016      zgtv_adout(:,:)  = gtv_ad(:,:)
1017      zgsv_adout(:,:)  = gsv_ad(:,:)
1018         
1019      zsp2_1 = DOT_PRODUCT( ztb_tlin , ztb_adout  )
1020      zsp2_2 = DOT_PRODUCT( zta_tlin , zta_adout  )
1021      zsp2_3 = DOT_PRODUCT( zgtu_tlin, zgtu_adout )
1022      zsp2_4 = DOT_PRODUCT( zgtv_tlin, zgtv_adout )
1023      zsp2_5 = DOT_PRODUCT( zsb_tlin , zsb_adout  )
1024      zsp2_6 = DOT_PRODUCT( zsa_tlin , zsa_adout  )
1025      zsp2_7 = DOT_PRODUCT( zgsu_tlin, zgsu_adout )
1026      zsp2_8 = DOT_PRODUCT( zgsv_tlin, zgsv_adout ) 
1027
1028      zsp2_T = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4
1029      zsp2_S = zsp2_5 + zsp2_6 + zsp2_7 + zsp2_8
1030      zsp2   = zsp2_T + zsp2_S
1031
1032      cl_name = 'tra_ldf_iso_ad'
1033      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1034
1035      DEALLOCATE(         &
1036         & ztb_tlin,      & ! Tangent input
1037         & zsb_tlin,      & ! Tangent input
1038         & zta_tlin,      & ! Tangent input
1039         & zsa_tlin,      & ! Tangent input
1040         & zgtu_tlin,     & ! Tangent input
1041         & zgsu_tlin,     & ! Tangent input
1042         & zgtv_tlin,     & ! Tangent input
1043         & zgsv_tlin,     & ! Tangent input
1044         & zta_tlout,     & ! Tangent output
1045         & zsa_tlout,     & ! Tangent output
1046         & zta_adin,      & ! Adjoint input
1047         & zsa_adin,      & ! Adjoint input
1048         & ztb_adout,     & ! Adjoint output
1049         & zsb_adout,     & ! Adjoint output
1050         & zta_adout,     & ! Adjoint output
1051         & zsa_adout,     & ! Adjoint output
1052         & zgtu_adout,    & ! Adjoint output
1053         & zgsu_adout,    & ! Adjoint output
1054         & zgtv_adout,    & ! Adjoint output
1055         & zgsv_adout,    & ! Adjoint output
1056         & z3r,           & ! 3D random field
1057         & z2r            &
1058         & )
1059     
1060   END SUBROUTINE tra_ldf_iso_adj_tst
1061# else
1062   !!----------------------------------------------------------------------
1063   !!   default option :   Dummy code   NO rotation of the diffusive tensor
1064   !!----------------------------------------------------------------------
1065CONTAINS
1066   SUBROUTINE tra_ldf_iso_tan( kt )               ! Empty routine
1067      WRITE(*,*) 'tra_ldf_iso_tan: You should not have seen this print! error?', kt
1068   END SUBROUTINE tra_ldf_iso_tan
1069   SUBROUTINE tra_ldf_iso_adj( kt )               ! Empty routine
1070      WRITE(*,*) 'tra_ldf_iso_adj: You should not have seen this print! error?', kt
1071   END SUBROUTINE tra_ldf_iso_adj
1072   SUBROUTINE tra_ldf_iso_adj_tst ( kumadt )
1073      WRITE(*,*) 'tra_ldf_iso_adj_tst: You should not have seen this print! error?', kt
1074   END SUBROUTINE tra_ldf_iso_adj_tst
1075# endif
1076#endif
1077
1078   !!==============================================================================
1079END MODULE traldf_iso_tam
Note: See TracBrowser for help on using the repository browser.