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_2_2/NEMOTAM/OPATAM_SRC/TRA – NEMO

source: branches/TAM_V3_2_2/NEMOTAM/OPATAM_SRC/TRA/traldf_iso_tam.F90 @ 3032

Last change on this file since 3032 was 2578, checked in by rblod, 13 years ago

first import of NEMOTAM 3.2.2

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