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

source: branches/TAM_V3_2_2/NEMOTAM/OPATAM_SRC/TRA/tradmp_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: 31.4 KB
Line 
1MODULE tradmp_tam
2#ifdef key_tam
3   !!======================================================================
4   !!                       ***  MODULE  tradmp_tam  ***
5   !! Ocean physics: internal restoring trend on active tracers (T and S)
6   !!                Tangent and Adjoint Module
7   !!======================================================================
8   !! History of the direct module:
9   !!            OPA  ! 1991-03  (O. Marti, G. Madec)  Original code
10   !!                 ! 1992-06  (M. Imbard)  doctor norme
11   !!                 ! 1996-01  (G. Madec)  statement function for e3
12   !!                 ! 1997-05  (G. Madec)  macro-tasked on jk-slab
13   !!                 ! 1998-07  (M. Imbard, G. Madec) ORCA version
14   !!            7.0  ! 2001-02  (M. Imbard)  cofdis, Original code
15   !!            8.1  ! 2001-02  (G. Madec, E. Durand)  cleaning
16   !!  NEMO      1.0  ! 2002-08  (G. Madec, E. Durand)  free form + modules
17   !!                 ! 2003-08  (M. Balmaseda)  mods to allow eq. damping
18   !!            3.2  ! 2009-08  (G. Madec, C. Talandier)  DOCTOR norm for namelist parameter
19   !! History of the TAM:
20   !!                 ! 2008-09  (A. Vidard) tangent and adjoint module
21   !!                                       of the 03-08 version
22   !!        NEMO 3.2 ! 2010-04 (F. Vigilant) 3.2 version   
23   !!----------------------------------------------------------------------
24#if   defined key_tradmp   ||   defined key_esopa
25   !!----------------------------------------------------------------------
26   !!   key_tradmp                                         internal damping
27   !!----------------------------------------------------------------------
28   !!----------------------------------------------------------------------
29   !!   tra_dmp      : update the tracer trend with the internal damping
30   !!   tra_dmp_init : initialization, namlist read, parameters control
31   !!   dtacof_zoom  : restoring coefficient for zoom domain
32   !!   dtacof       : restoring coefficient for global domain
33   !!   cofdis       : compute the distance to the coastline
34   !!----------------------------------------------------------------------
35   USE par_kind      , ONLY: & ! Precision variables
36      & wp
37   USE par_oce       , ONLY: & ! Ocean space and time domain variables
38      & jpi,                 &
39      & jpj,                 & 
40      & jpk,                 &
41      & jpim1,               &
42      & jpjm1,               &
43      & jpkm1,               &
44      & jpiglo,              &
45      & jpjglo,              &
46      & jpizoom,             &
47      & jpjzoom
48   USE oce_tam       , ONLY: & ! ocean dynamics and tracers variables
49      & tb_tl,               &
50      & sb_tl,               &
51      & ta_tl,               &
52      & sa_tl,               &
53      & tb_ad,               &
54      & sb_ad,               &
55      & ta_ad,               &
56      & sa_ad     
57   USE dom_oce       , ONLY: & ! Ocean space and time domain
58      & lzoom,               & 
59      & lzoom_e,             & 
60      & lzoom_w,             & 
61      & lzoom_s,             & 
62      & lzoom_n,             & 
63      & lzoom_arct,          & 
64      & lzoom_anta,          &
65      & nperio,              &
66      & rdt,                 &
67      & mi0,                 &
68      & mj0,                 &
69      & mi1,                 &
70      & mj1,                 &
71      & mig,                 &
72      & mjg,                 &
73      & nldi,                &
74      & nldj,                &
75      & nlei,                &
76      & nlej,                &
77      & e1t,                 &
78      & e2t,                 &
79#if defined key_zco
80      & e3t_0,               &
81      & gdept_0,             &
82#else
83      & e3t,                 &
84      & gdept,               &
85#endif
86      & gdept_0,             &
87      & umask,               &
88      & vmask,               &
89      & tmask,               &
90      & fmask,               &
91      & gphiu,               &
92      & glamu,               &
93      & gphiv,               &
94      & glamv,               &
95      & gphit,               &
96      & glamt,               &
97      & gphif,               &
98      & glamf
99   USE zdf_oce, ONLY: &! ocean vertical physics
100      & avt
101   USE in_out_manager, ONLY: & ! I/O manager
102      & ctl_stop,            &
103      & lwp,                 &
104      & numout,              & 
105      & numnam,              & 
106      & ctmp1,               & 
107      & nit000,              &
108      & nitend
109   USE phycst        , ONLY: & ! Define parameters for the routines
110      & rday,                & 
111      & rpi,                 &
112      & ra,                  &
113      & rad
114   USE lib_mpp       , ONLY: & ! distributed memory computing
115      & lk_mpp
116   USE tradmp        , ONLY : &
117      & dtacof_zoom, dtacof,  &
118      & cofdis
119   USE prtctl        , ONLY: & ! Print control
120      & prt_ctl
121   USE gridrandom    , ONLY: & ! Random Gaussian noise on grids
122      & grid_random
123   USE dotprodfld    , ONLY: & ! Computes dot product for 3D and 2D fields
124      & dot_product
125   USE dtatem        , ONLY: & ! temperature data
126      & lk_dtatem
127   USE dtasal        , ONLY: & ! salinity data
128      & lk_dtasal
129   USE zdfmxl        , ONLY: & ! mixed layer depth
130      & hmlp
131   USE tstool_tam    , ONLY: &
132      & prntst_adj,          & !
133      & stdt,                & ! stdev for temperature
134      & stds                   !           salinity
135   USE paresp        , ONLY: &
136      & wesp_t,              &
137      & wesp_s
138   IMPLICIT NONE
139   PRIVATE
140
141   PUBLIC tra_dmp_tan      ! routine called by step_tam.F90
142   PUBLIC tra_dmp_adj      ! routine called by step_tam.F90
143   PUBLIC tra_dmp_adj_tst  ! routine called by tst.F90
144#if ! defined key_agrif
145   LOGICAL, PUBLIC, PARAMETER ::   lk_tradmp = .TRUE.     !: internal damping flag
146#else
147   LOGICAL, PUBLIC            ::   lk_tradmp = .TRUE.     !: internal damping flag
148#endif
149   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: &
150     & strdmp_tl,   & !: damping salinity trend (psu/s)
151     & ttrdmp_tl,   & !: damping temperature trend (psu/s)
152     & ttrdmp_ad,   & !: damping temperature trend (psu/s)
153     & strdmp_ad,   & !: damping salinity trend (psu/s)
154     & resto          !: restoring coeff. on T and S (s-1)
155   
156   
157   LOGICAL :: lfirst = .TRUE.     !: flag for initialisation
158   !                             !!* Namelist namtra_dmp : T & S newtonian damping *
159   INTEGER  ::   nn_hdmp =   -1   ! = 0/-1/'latitude' for damping over T and S
160   INTEGER  ::   nn_zdmp =    0   ! = 0/1/2 flag for damping in the mixed layer
161   REAL(wp) ::   rn_surf =   50.  ! surface time scale for internal damping        [days]
162   REAL(wp) ::   rn_bot  =  360.  ! bottom time scale for internal damping         [days]
163   REAL(wp) ::   rn_dep  =  800.  ! depth of transition between rn_surf and rn_bot [meters]
164   INTEGER  ::   nn_file =    2   ! = 1 create a damping.coeff NetCDF file
165   !! * Substitutions
166#  include "domzgr_substitute.h90"
167#  include "vectopt_loop_substitute.h90"
168
169
170CONTAINS
171
172   SUBROUTINE tra_dmp_tan( kt )
173      !!----------------------------------------------------------------------
174      !!                   ***  ROUTINE tra_dmp_tan  ***
175      !!                 
176      !! ** Purpose of direct routine:
177      !!      Compute the tracer trend due to a newtonian damping
178      !!      of the tracer field towards given data field and add it to the
179      !!      general tracer trends.
180      !!
181      !! ** Method  :   Newtonian damping towards t_dta and s_dta computed
182      !!      and add to the general tracer trends:
183      !!                     ta = ta + resto * (t_dta - tb)
184      !!                     sa = sa + resto * (s_dta - sb)
185      !!         The trend is computed either throughout the water column
186      !!      (nlmdmp=0) or in area of weak vertical mixing (nlmdmp=1) or
187      !!      below the well mixed layer (nlmdmp=2)
188      !!
189      !! ** Action  : - (ta,sa)   tracer trends updated with the damping trend
190      !!
191      !! ASSUME key_zdfcst_tam
192      !!----------------------------------------------------------------------
193      !!
194      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
195      !!
196      INTEGER  ::   ji, jj, jk            ! dummy loop indices
197      REAL(wp) ::   ztest, ztatl, zsatl       ! temporary scalars
198      !!----------------------------------------------------------------------
199      IF( kt == nit000 )   CALL tra_dmp_init_tam  ! Initialization
200
201      ! 1. Newtonian damping trends on tracer fields
202      ! --------------------------------------------
203      !    compute the newtonian damping trends depending on nmldmp
204
205      SELECT CASE ( nn_zdmp )
206      !
207      CASE( 0 )                ! newtonian damping throughout the water column
208         DO jk = 1, jpkm1
209            DO jj = 2, jpjm1
210               DO ji = fs_2, fs_jpim1   ! vector opt.
211                  ztatl = - resto(ji,jj,jk) *  tb_tl(ji,jj,jk) 
212                  zsatl = - resto(ji,jj,jk) *  sb_tl(ji,jj,jk) 
213                  ! add the trends to the general tracer trends
214                  ta_tl(ji,jj,jk) = ta_tl(ji,jj,jk) + ztatl
215                  sa_tl(ji,jj,jk) = sa_tl(ji,jj,jk) + zsatl
216                  ! save the salinity trend (used in flx to close the salt budget)
217                  strdmp_tl(ji,jj,jk) = zsatl
218                  ttrdmp_tl(ji,jj,jk) = ztatl
219               END DO
220            END DO
221         END DO
222         !
223      CASE ( 1 )                ! no damping in the turbocline (avt > 5 cm2/s)
224         DO jk = 1, jpkm1
225            DO jj = 2, jpjm1
226               DO ji = fs_2, fs_jpim1   ! vector opt.
227                  ztest = avt(ji,jj,jk) - 5.e-4
228      !! ASSUME key_zdfcst_tam
229                  IF( ztest < 0. ) THEN
230                     ztatl = - resto(ji,jj,jk) * tb_tl(ji,jj,jk) 
231                     zsatl = - resto(ji,jj,jk) * sb_tl(ji,jj,jk)
232                  ELSE
233                     ztatl = 0.0_wp
234                     zsatl = 0.0_wp
235                  ENDIF
236                  ! add the trends to the general tracer trends
237                  ta_tl(ji,jj,jk) = ta_tl(ji,jj,jk) + ztatl
238                  sa_tl(ji,jj,jk) = sa_tl(ji,jj,jk) + zsatl
239                  ! save the salinity trend (used in flx to close the salt budget)
240                  strdmp_tl(ji,jj,jk) = zsatl
241                  ttrdmp_tl(ji,jj,jk) = ztatl
242               END DO
243            END DO
244         END DO
245         !
246      CASE ( 2 )                ! no damping in the mixed layer 
247         DO jk = 1, jpkm1
248            DO jj = 2, jpjm1
249               DO ji = fs_2, fs_jpim1   ! vector opt.
250                  IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
251                     ztatl = - resto(ji,jj,jk) * tb_tl(ji,jj,jk)
252                     zsatl = - resto(ji,jj,jk) * sb_tl(ji,jj,jk)
253                  ELSE
254                     ztatl = 0.e0
255                     zsatl = 0.e0
256                  ENDIF
257                  ! add the trends to the general tracer trends
258                  ta_tl(ji,jj,jk) = ta_tl(ji,jj,jk) + ztatl
259                  sa_tl(ji,jj,jk) = sa_tl(ji,jj,jk) + zsatl
260                  ! save the salinity trend (used in flx to close the salt budget)
261                  strdmp_tl(ji,jj,jk) = zsatl
262                  ttrdmp_tl(ji,jj,jk) = ztatl
263               END DO
264            END DO
265         END DO
266         !
267      END SELECT
268
269   END SUBROUTINE tra_dmp_tan
270   SUBROUTINE tra_dmp_adj( kt )
271      !!----------------------------------------------------------------------
272      !!                   ***  ROUTINE tra_dmp_adj  ***
273      !!                 
274      !! ** Purpose of direct routine:
275      !!      Compute the tracer trend due to a newtonian damping
276      !!      of the tracer field towards given data field and add it to the
277      !!      general tracer trends.
278      !!
279      !! ** Method  :   Newtonian damping towards t_dta and s_dta computed
280      !!      and add to the general tracer trends:
281      !!                     ta = ta + resto * (t_dta - tb)
282      !!                     sa = sa + resto * (s_dta - sb)
283      !!         The trend is computed either throughout the water column
284      !!      (nlmdmp=0) or in area of weak vertical mixing (nlmdmp=1) or
285      !!      below the well mixed layer (nlmdmp=2)
286      !!
287      !! ** Action  : - (ta,sa)   tracer trends updated with the damping trend
288      !!              - save the trends in (ttrd,strd) ('key_trdtra')
289      !! ASSUME key_zdfcst_tam
290      !!----------------------------------------------------------------------
291      !!
292      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
293      !!
294      INTEGER  ::   ji, jj, jk            ! dummy loop indices
295      REAL(wp) ::   ztest, ztaad, zsaad       ! temporary scalars
296      !!----------------------------------------------------------------------
297      IF( kt == nitend )   CALL tra_dmp_init_tam  ! Initialization
298
299      ! 1. Newtonian damping trends on tracer fields
300      ! --------------------------------------------
301      !    compute the newtonian damping trends depending on nmldmp
302      ztaad = 0.0_wp
303      zsaad = 0.0_wp
304
305      SELECT CASE ( nn_zdmp )
306      !
307      CASE( 0 )                ! newtonian damping throughout the water column
308         DO jk = 1, jpkm1
309            DO jj = 2, jpjm1
310               DO ji = fs_2, fs_jpim1   ! vector opt.
311                  ! save the salinity trend (used in flx to close the salt budget)
312                  ztaad = ztaad + ttrdmp_ad(ji,jj,jk) 
313                  ttrdmp_ad(ji,jj,jk) = 0.0_wp
314                  zsaad = zsaad + strdmp_ad(ji,jj,jk) 
315                  strdmp_ad(ji,jj,jk) = 0.0_wp
316                  ! add the trends to the general tracer trends
317                  ztaad = ta_ad(ji,jj,jk) + ztaad
318                  zsaad = sa_ad(ji,jj,jk) + zsaad
319
320                  tb_ad(ji,jj,jk)  = tb_ad(ji,jj,jk) - ztaad * resto(ji,jj,jk)
321                  sb_ad(ji,jj,jk)  = sb_ad(ji,jj,jk) - zsaad * resto(ji,jj,jk)
322                  ztaad = 0.0_wp
323                  zsaad = 0.0_wp
324                END DO
325            END DO
326         END DO
327         !
328      CASE ( 1 )                ! no damping in the turbocline (avt > 5 cm2/s)
329         DO jk = 1, jpkm1
330            DO jj = 2, jpjm1
331               DO ji = fs_2, fs_jpim1   ! vector opt.
332                  ! save the salinity trend (used in flx to close the salt budget)
333                  ztaad = ztaad + ttrdmp_ad(ji,jj,jk) 
334                  ttrdmp_ad(ji,jj,jk) = 0.0_wp
335                  zsaad = zsaad + strdmp_ad(ji,jj,jk) 
336                  strdmp_ad(ji,jj,jk) = 0.0_wp
337                  ! add the trends to the general tracer trends
338                  ztaad = ta_ad(ji,jj,jk) + ztaad
339                  zsaad = sa_ad(ji,jj,jk) + zsaad
340                  ztest = avt(ji,jj,jk) - 5.e-4
341      !! ASSUME key_zdfcst_tam
342                  IF( ztest < 0. ) THEN
343                     tb_ad(ji,jj,jk)  = tb_ad(ji,jj,jk) - ztaad * resto(ji,jj,jk)
344                     sb_ad(ji,jj,jk)  = sb_ad(ji,jj,jk) - zsaad * resto(ji,jj,jk)
345                     ztaad = 0.0_wp
346                     zsaad = 0.0_wp
347                  ELSE
348                     ztaad = 0.0_wp
349                     zsaad = 0.0_wp
350                  ENDIF
351               END DO
352            END DO
353         END DO
354         !
355      CASE ( 2 )                ! no damping in the mixed layer
356         DO jk = 1, jpkm1
357            DO jj = 2, jpjm1
358               DO ji = fs_2, fs_jpim1   ! vector opt.
359                  ! save the salinity trend (used in flx to close the salt budget)
360                  ztaad = ztaad + ttrdmp_ad(ji,jj,jk) 
361                  ttrdmp_ad(ji,jj,jk) = 0.0_wp
362                  zsaad = zsaad + strdmp_ad(ji,jj,jk) 
363                  strdmp_ad(ji,jj,jk) = 0.0_wp
364                  ! add the trends to the general tracer trends
365                  ztaad = ta_ad(ji,jj,jk) + ztaad
366                  zsaad = sa_ad(ji,jj,jk) + zsaad
367                  IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
368                     tb_ad(ji,jj,jk)  = tb_ad(ji,jj,jk) - ztaad * resto(ji,jj,jk)
369                     sb_ad(ji,jj,jk)  = sb_ad(ji,jj,jk) - zsaad * resto(ji,jj,jk)
370                     ztaad = 0.0_wp
371                     zsaad = 0.0_wp
372                  ELSE
373                     ztaad = 0.0_wp
374                     zsaad = 0.0_wp
375                  ENDIF
376               END DO
377            END DO
378         END DO
379         !
380      END SELECT
381
382
383   END SUBROUTINE tra_dmp_adj
384   SUBROUTINE tra_dmp_init_tam
385      !!----------------------------------------------------------------------
386      !!                  ***  ROUTINE tra_dmp_init_tam  ***
387      !!
388      !! ** Purpose :   Initialization for the newtonian damping
389      !!
390      !! ** Method  :   read the nammbf namelist and check the parameters
391      !!      called by tra_dmp at the first timestep (nit000)
392      !!----------------------------------------------------------------------
393
394      NAMELIST/namtra_dmp/ nn_hdmp, nn_file, nn_zdmp, rn_surf, rn_bot, rn_dep
395      !!----------------------------------------------------------------------
396
397      IF (lfirst) THEN
398
399         REWIND ( numnam )                  ! Read Namelist namtdp : temperature and salinity damping term
400         READ   ( numnam, namtra_dmp )
401         IF( lzoom )   nn_zdmp = 0           ! restoring to climatology at closed north or south boundaries
402         
403         IF(lwp) THEN                       ! Namelist print
404            WRITE(numout,*)
405            WRITE(numout,*) 'tra_dmp_tam : T and S newtonian damping'
406            WRITE(numout,*) '~~~~~~~~~~~'
407            WRITE(numout,*) '      Namelist namtra_dmp : set damping parameter'
408            WRITE(numout,*) '      T and S damping option         nn_hdmp = ', nn_hdmp
409            WRITE(numout,*) '      mixed layer damping option     nn_zdmp = ', nn_zdmp, '(zoom: forced to 0)'
410            WRITE(numout,*) '      surface time scale (days)      rn_surf = ', rn_surf
411            WRITE(numout,*) '      bottom time scale (days)       rn_bot  = ', rn_bot
412            WRITE(numout,*) '      depth of transition (meters)   rn_dep  = ', rn_dep
413            WRITE(numout,*) '      create a damping.coeff file    nn_file = ', nn_file
414         ENDIF
415         
416         !! ** modified to allow damping at the equator
417         
418         SELECT CASE ( nn_hdmp )
419         CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   tracer damping in the Med & Red seas only'
420         CASE ( 0:90 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping poleward of', nn_hdmp, ' degrees'
421         CASE DEFAULT
422            IF(lwp) WRITE(numout,*) '          tracer damping disabled', nn_hdmp
423            !         WRITE(ctmp1,*) '          bad flag value for ndmp = ', ndmp
424            !         CALL ctl_stop(ctmp1)
425         END SELECT
426         
427         SELECT CASE ( nn_zdmp )
428         CASE ( 0 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping throughout the water column'
429         CASE ( 1 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline (avt > 5 cm2/s)'
430         CASE ( 2 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer'
431         CASE DEFAULT
432            WRITE(ctmp1,*) '   bad flag value for nn_zdmp = ', nn_zdmp
433            CALL ctl_stop(ctmp1)
434         END SELECT
435         
436         IF( .NOT.lk_dtasal .OR. .NOT.lk_dtatem )   &
437            &   CALL ctl_stop( 'no temperature and/or salinity data define key_dtatem and key_dtasal' )
438         
439         !                          ! Damping coefficients initialization
440         IF( lzoom ) THEN   ;   CALL dtacof_zoom
441         ELSE               ;   CALL dtacof
442         ENDIF
443         !
444         lfirst = .FALSE.
445      END IF
446      ttrdmp_tl(:,:,:) = 0.0_wp       ! internal damping salinity trend (used in ocesbc)
447      ttrdmp_ad(:,:,:) = 0.0_wp       ! internal damping salinity trend (used in ocesbc)
448      strdmp_tl(:,:,:) = 0.0_wp       ! internal damping salinity trend (used in ocesbc)
449      strdmp_ad(:,:,:) = 0.0_wp       ! internal damping salinity trend (used in ocesbc)
450   END SUBROUTINE tra_dmp_init_tam
451
452   SUBROUTINE tra_dmp_adj_tst ( kumadt ) 
453      !!-----------------------------------------------------------------------
454      !!
455      !!          ***  ROUTINE tra_sbc_adj_tst : TEST OF tra_sbc_adj  ***
456      !!
457      !! ** Purpose : Test the adjoint routine.
458      !!
459      !! ** Method  : Verify the scalar product
460      !!           
461      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
462      !!
463      !!              where  L   = tangent routine
464      !!                     L^T = adjoint routine
465      !!                     W   = diagonal matrix of scale factors
466      !!                     dx  = input perturbation (random field)
467      !!                     dy  = L dx
468      !!
469       !! History :
470      !!        ! 08-08 (A. Vidard)
471      !!-----------------------------------------------------------------------
472      !! * Modules used
473
474      !! * Arguments
475      INTEGER, INTENT(IN) :: &
476         & kumadt             ! Output unit
477 
478      INTEGER ::  &
479         & istp,  &
480         & jstp,  &
481         & ji,    &        ! dummy loop indices
482         & jj,    &       
483         & jk     
484      INTEGER, DIMENSION(jpi,jpj) :: &
485         & iseed_2d        ! 2D seed for the random number generator
486
487      !! * Local declarations
488      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
489         & zsb_tlin,     &! Tangent input : before salinity
490         & ztb_tlin,     &! Tangent input : before temperature
491         & zsa_tlin,     &! Tangent input : after salinity
492         & zta_tlin,     &! Tangent input : after temperature
493         & zsa_tlout,    &! Tangent output: after salinity
494         & zta_tlout,    &! Tangent output: after temperature
495         & zstrdmp_tlout,&! Tangent output: salinity trend
496         & zttrdmp_tlout,&! Tangent output: temperature trend
497         & zsb_adout,    &! Adjoint output : before salinity
498         & ztb_adout,    &! Adjoint output : before temperature
499         & zsa_adout,    &! Adjoint output : after salinity
500         & zta_adout,    &! Adjoint output : after temperature
501         & zsa_adin,     &! Adjoint input : after salinity
502         & zta_adin,     &! Adjoint input : after temperature
503         & zstrdmp_adin, &! Adjoint input : salinity trend
504         & zttrdmp_adin, &! Tangent output: temperature trend
505         & z3r            ! 3D field
506         
507      REAL(KIND=wp) ::       &
508         & zsp1,             & ! scalar product involving the tangent routine
509         & zsp1_1,           & ! scalar product involving the tangent routine
510         & zsp1_2,           & ! scalar product involving the tangent routine
511         & zsp1_3,           & ! scalar product involving the tangent routine
512         & zsp1_4,           & ! scalar product involving the tangent routine
513         & zsp2,             & ! scalar product involving the adjoint routine
514         & zsp2_1,           & ! scalar product involving the adjoint routine
515         & zsp2_2,           & ! scalar product involving the adjoint routine
516         & zsp2_3,           & ! scalar product involving the adjoint routine
517         & zsp2_4              ! scalar product involving the adjoint routine
518      CHARACTER(LEN=14) ::   &
519         & cl_name
520
521      ALLOCATE( & 
522         & ztb_tlin(jpi,jpj,jpk),     &
523         & zsb_tlin(jpi,jpj,jpk),     &
524         & zsa_tlin(jpi,jpj,jpk),     &
525         & zta_tlin(jpi,jpj,jpk),     &
526         & zsa_tlout(jpi,jpj,jpk),    &
527         & zta_tlout(jpi,jpj,jpk),    &
528         & zstrdmp_tlout(jpi,jpj,jpk),&
529         & zttrdmp_tlout(jpi,jpj,jpk),&
530         & ztb_adout(jpi,jpj,jpk),    &
531         & zsb_adout(jpi,jpj,jpk),    &
532         & zsa_adout(jpi,jpj,jpk),    &
533         & zta_adout(jpi,jpj,jpk),    &
534         & zsa_adin(jpi,jpj,jpk),     &
535         & zta_adin(jpi,jpj,jpk),     &
536         & zttrdmp_adin(jpi,jpj,jpk), &
537         & zstrdmp_adin(jpi,jpj,jpk), &
538         & z3r         (jpi,jpj,jpk)  &
539         & )
540      ! Test for time steps nit000 and nit000 + 1 (the matrix changes)
541
542      DO jstp = nit000, nit000 + 2
543         istp = jstp
544         IF ( jstp == nit000 +2 ) istp = nitend
545
546      ! Initialize the reference state
547      avt(:,:,:) = 1.e-1
548      !=============================================================
549      ! 1) dx = ( T ) and dy = ( T )
550      !=============================================================
551
552      !--------------------------------------------------------------------
553      ! Reset the tangent and adjoint variables
554      !--------------------------------------------------------------------
555      zsb_tlin(:,:,:)      = 0.0_wp 
556      ztb_tlin(:,:,:)      = 0.0_wp 
557      zsa_tlin(:,:,:)      = 0.0_wp 
558      zta_tlin(:,:,:)      = 0.0_wp 
559      zsa_tlout(:,:,:)     = 0.0_wp 
560      zta_tlout(:,:,:)     = 0.0_wp 
561      zstrdmp_tlout(:,:,:) = 0.0_wp
562      zttrdmp_tlout(:,:,:) = 0.0_wp 
563      zsb_adout(:,:,:)     = 0.0_wp 
564      ztb_adout(:,:,:)     = 0.0_wp 
565      zsa_adout(:,:,:)     = 0.0_wp 
566      zta_adout(:,:,:)     = 0.0_wp 
567      zsa_adin(:,:,:)      = 0.0_wp 
568      zta_adin(:,:,:)      = 0.0_wp 
569      zstrdmp_adin(:,:,:)  = 0.0_wp 
570      zttrdmp_adin(:,:,:)  = 0.0_wp
571
572      sa_ad(:,:,:)      = 0.0_wp 
573      sb_ad(:,:,:)      = 0.0_wp 
574      ta_ad(:,:,:)      = 0.0_wp 
575      tb_ad(:,:,:)      = 0.0_wp 
576      strdmp_ad(:,:,:)  = 0.0_wp
577      ttrdmp_ad(:,:,:)  = 0.0_wp 
578      sa_tl(:,:,:)      = 0.0_wp 
579      sb_tl(:,:,:)      = 0.0_wp 
580      ta_tl(:,:,:)      = 0.0_wp 
581      tb_tl(:,:,:)      = 0.0_wp 
582      strdmp_tl(:,:,:)  = 0.0_wp
583      ttrdmp_tl(:,:,:)  = 0.0_wp 
584
585      DO jj = 1, jpj
586         DO ji = 1, jpi
587            iseed_2d(ji,jj) = - ( 785483 + &
588               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
589         END DO
590      END DO
591      CALL grid_random( iseed_2d, z3r, 'T', 0.0_wp, stds )
592      DO jk = 1, jpk
593         DO jj = nldj, nlej
594            DO ji = nldi, nlei
595               zsb_tlin(ji,jj,jk) = z3r(ji,jj,jk) 
596            END DO
597         END DO
598      END DO
599      DO jj = 1, jpj
600         DO ji = 1, jpi
601            iseed_2d(ji,jj) = - ( 358606 + &
602               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
603         END DO
604      END DO
605      CALL grid_random( iseed_2d, z3r, 'T', 0.0_wp, stdt )
606      DO jk = 1, jpk
607         DO jj = nldj, nlej
608            DO ji = nldi, nlei
609               ztb_tlin(ji,jj,jk) = z3r(ji,jj,jk) 
610            END DO
611         END DO
612      END DO
613
614      DO jj = 1, jpj
615         DO ji = 1, jpi
616            iseed_2d(ji,jj) = - ( 596035 + &
617               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
618         END DO
619      END DO
620      CALL grid_random( iseed_2d, z3r, 'T', 0.0_wp, stds )
621      DO jk = 1, jpk
622         DO jj = nldj, nlej
623            DO ji = nldi, nlei
624               zsa_tlin(ji,jj,jk) = z3r(ji,jj,jk) 
625            END DO
626         END DO
627      END DO
628
629      DO jj = 1, jpj
630         DO ji = 1, jpi
631            iseed_2d(ji,jj) = - ( 523432 + &
632               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
633         END DO
634      END DO
635      CALL grid_random( iseed_2d, z3r, 'T', 0.0_wp, stdt )
636      DO jk = 1, jpk
637         DO jj = nldj, nlej
638            DO ji = nldi, nlei
639               zta_tlin(ji,jj,jk) = z3r(ji,jj,jk) 
640            END DO
641         END DO
642      END DO
643
644      sb_tl(:,:,:) = zsb_tlin(:,:,:) 
645      tb_tl(:,:,:) = ztb_tlin(:,:,:) 
646      sa_tl(:,:,:) = zsa_tlin(:,:,:) 
647      ta_tl(:,:,:) = zta_tlin(:,:,:)
648     
649      CALL tra_dmp_tan( istp )
650
651      zsa_tlout(:,:,:)     = sa_tl(:,:,:) 
652      zta_tlout(:,:,:)     = ta_tl(:,:,:)
653      zstrdmp_tlout(:,:,:) = strdmp_tl(:,:,:)
654      zttrdmp_tlout(:,:,:) = ttrdmp_tl(:,:,:)
655
656      !--------------------------------------------------------------------
657      ! Initialize the adjoint variables: dy^* = W dy
658      !--------------------------------------------------------------------
659
660      DO jk = 1, jpk
661        DO jj = nldj, nlej
662           DO ji = nldi, nlei
663              zsa_adin(ji,jj,jk)     = zsa_tlout(ji,jj,jk) &
664                 &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
665                 &                   * tmask(ji,jj,jk) * wesp_s(jk)
666              zta_adin(ji,jj,jk)     = zta_tlout(ji,jj,jk) &
667                 &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
668                 &                   * tmask(ji,jj,jk) * wesp_t(jk)
669              zstrdmp_adin(ji,jj,jk) = zstrdmp_tlout(ji,jj,jk) &
670                 &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
671                 &                   * tmask(ji,jj,jk) * wesp_s(jk)
672              zttrdmp_adin(ji,jj,jk) = zttrdmp_tlout(ji,jj,jk) &                 
673                 &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
674                 &                   * tmask(ji,jj,jk) * wesp_t(jk)
675            END DO
676         END DO
677      END DO
678
679      !--------------------------------------------------------------------
680      ! Compute the scalar product: ( L dx )^T W dy
681      !--------------------------------------------------------------------
682
683      zsp1_1 = DOT_PRODUCT( zsa_tlout    , zsa_adin     )
684      zsp1_2 = DOT_PRODUCT( zta_tlout    , zta_adin     )
685      zsp1_3 = DOT_PRODUCT( zstrdmp_tlout, zstrdmp_adin )
686      zsp1_4 = DOT_PRODUCT( zttrdmp_tlout, zttrdmp_adin ) 
687      zsp1   = zsp1_1 + zsp1_2 + zsp1_3 + zsp1_4
688
689      !--------------------------------------------------------------------
690      ! Call the adjoint routine: dx^* = L^T dy^*
691      !--------------------------------------------------------------------
692
693      sa_ad(:,:,:)     = zsa_adin(:,:,:)
694      ta_ad(:,:,:)     = zta_adin(:,:,:)
695      strdmp_ad(:,:,:) = zstrdmp_adin(:,:,:)
696      ttrdmp_ad(:,:,:) = zttrdmp_adin(:,:,:)
697
698      CALL tra_dmp_adj( istp )
699
700      zsb_adout(:,:,:) = sb_ad(:,:,:)
701      ztb_adout(:,:,:) = tb_ad(:,:,:)
702      zsa_adout(:,:,:) = sa_ad(:,:,:)
703      zta_adout(:,:,:) = ta_ad(:,:,:)
704
705      !--------------------------------------------------------------------
706      ! Compute the scalar product: dx^T L^T W dy
707      !--------------------------------------------------------------------
708
709      zsp2_1 = DOT_PRODUCT( zsb_tlin  , zsb_adout   )
710      zsp2_2 = DOT_PRODUCT( ztb_tlin  , ztb_adout   )
711      zsp2_3 = DOT_PRODUCT( zsa_tlin  , zsa_adout   )
712      zsp2_4 = DOT_PRODUCT( zta_tlin  , zta_adout   )
713
714      zsp2 = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4
715
716      ! Compare the scalar products
717
718      ! 14 char:'12345678901234'
719      IF ( istp == nit000 ) THEN
720         cl_name = 'tra_dmp_adj T1'
721      ELSEIF ( istp == nit000 +1 ) THEN
722         cl_name = 'tra_dmp_adj T2'
723      ELSEIF ( istp == nitend ) THEN
724         cl_name = 'tra_dmp_adj T3'
725      END IF
726      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
727
728      END DO
729
730      DEALLOCATE( & 
731         & zsb_tlin,     &
732         & ztb_tlin,     &
733         & zsa_tlin,     &
734         & zta_tlin,     &
735         & zsa_tlout,    &
736         & zta_tlout,    &
737         & zstrdmp_tlout,&
738         & zttrdmp_tlout,&
739         & zsb_adout,    &
740         & ztb_adout,    &
741         & zsa_adout,    &
742         & zta_adout,    &
743         & zsa_adin,     &
744         & zta_adin,     &
745         & zstrdmp_adin, &
746         & zttrdmp_adin, &
747         & z3r           &
748         & )
749
750   END SUBROUTINE tra_dmp_adj_tst
751
752#else
753   !!----------------------------------------------------------------------
754   !!   Default key                                     NO internal damping
755   !!----------------------------------------------------------------------
756   LOGICAL , PUBLIC, PARAMETER ::   lk_tradmp = .FALSE.    !: internal damping flag
757CONTAINS
758   SUBROUTINE tra_dmp_tan( kt )        ! Empty routine
759      WRITE(*,*) 'tra_dmp_tan: You should not have seen this print! error?', kt
760   END SUBROUTINE tra_dmp_tan
761   SUBROUTINE tra_dmp_adj( kt )        ! Empty routine
762      WRITE(*,*) 'tra_dmp_adj: You should not have seen this print! error?', kt
763   END SUBROUTINE tra_dmp_adj
764   SUBROUTINE tra_dmp_adj_tst( kt )        ! Empty routine
765      WRITE(*,*) 'tra_dmp_adj_tst: You should not have seen this print! error?', kt
766   END SUBROUTINE tra_dmp_adj_tst
767#endif
768
769   !!======================================================================
770#endif
771END MODULE tradmp_tam
Note: See TracBrowser for help on using the repository browser.