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

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/DYN/dynzdf_tam.F90 @ 1885

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

add TAM sources

File size: 17.2 KB
Line 
1MODULE dynzdf_tam
2#ifdef key_tam
3   !!==============================================================================
4   !!                 ***  MODULE  dynzdf_tam  ***
5   !! Ocean dynamics :  vertical component of the momentum mixing trend
6   !!                   Tangent and adjoint module
7   !!==============================================================================
8   !! History of the direct module:
9   !!         9.0  !  05-11  (G. Madec)  Original code
10   !! History of the T&A module:
11   !!         9.0  !  08-06  (A. Vidard) Skeleton
12   !!         9.0  !  08-08  (A. Vidard) tam of the 05-11 version
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   dyn_zdf      : Update the momentum trend with the vertical diffusion
17   !!       zdf_ctl  : initializations of the vertical diffusion scheme
18   !!----------------------------------------------------------------------
19   USE par_kind      , ONLY: & ! Precision variables
20      & wp
21   USE par_oce       , ONLY: & ! Ocean space and time domain variables
22      & jpi,                 &
23      & jpj,                 & 
24      & jpk,                 & 
25      & jpiglo
26   USE oce_tam       , ONLY: & ! ocean dynamics and tracers tam variables
27      & ub_tl,               &
28      & vb_tl,               &
29      & ub_ad,               &
30      & vb_ad,               &
31      & ua_tl,               &
32      & va_tl,               &
33      & ua_ad,               &
34      & va_ad
35   USE dom_oce       , ONLY: & ! ocean space and time domain variables
36      & rdt,                 &
37      & neuler,              &
38      & e1u,                 &
39      & e2u,                 &
40      & e1v,                 &
41      & e2v,                 &
42#if defined key_zco
43      & e3t_0,               &
44#else
45      & e3u,                 &
46      & e3v,                 &
47#endif
48      & mig,                 &
49      & mjg,                 &
50      & nldi,                &
51      & nldj,                &
52      & nlei,                &
53      & nlej,                &
54      & umask,               &
55      & vmask,               &
56      & ln_sco,              &
57      & lk_esopa
58   USE zdf_oce       , ONLY: & ! ocean vertical physics
59      & avmu,                &
60      & avmv,                &
61      & avm0,                &
62      & ln_zdfexp
63   USE dynzdf_exp_tam, ONLY: & ! vertical diffusion: explicit (dyn_zdf_exp     routine)
64      & dyn_zdf_exp_tan,     &
65      & dyn_zdf_exp_adj
66   USE dynzdf_imp_tam, ONLY: & ! vertical diffusion: implicit (dyn_zdf_imp     routine)
67      & dyn_zdf_imp_tan,     &
68      & dyn_zdf_imp_adj
69   USE ldfdyn_oce    , ONLY: & ! ocean dynamics: lateral physics
70      & ln_dynldf_iso,       &
71      & ln_dynldf_hor
72   USE in_out_manager, ONLY: & ! I/O manager
73      & numout,              &
74      & nit000,              &
75      & nitend,              &
76      & lwp
77   USE gridrandom    , ONLY: & ! Random Gaussian noise on grids
78      & grid_random
79   USE dotprodfld,     ONLY: & ! Computes dot product for 3D and 2D fields
80      & dot_product
81   USE tstool_tam    , ONLY: &
82      & prntst_adj,          & !
83                               ! random field standard deviation for:
84      & stdu,                & !   u-velocity
85      & stdv
86   IMPLICIT NONE
87   PRIVATE
88
89   PUBLIC   dyn_zdf_tan    !  routine called by step_tam.F90
90   PUBLIC   dyn_zdf_adj    !  routine called by step_tam.F90
91   PUBLIC   dyn_zdf_adj_tst!  routine called by tst.F90
92
93   INTEGER  ::   nzdf = 0              ! type vertical diffusion algorithm used
94      !                                ! defined from ln_zdf...  namlist logicals)
95
96   REAL(wp) ::   r2dt                  ! time-step, = 2 rdttra
97      !                                ! except at nit000 (=rdttra) if neuler=0
98   LOGICAL :: lfirst = .TRUE.
99
100   !! * Substitutions
101#  include "domzgr_substitute.h90"
102#  include "zdfddm_substitute.h90"
103#  include "vectopt_loop_substitute.h90"
104
105CONTAINS
106   
107   SUBROUTINE dyn_zdf_tan( kt )
108      !!----------------------------------------------------------------------
109      !!                  ***  ROUTINE dyn_zdf_tan  ***
110      !!
111      !! ** Purpose of the direct routine:
112      !!            compute the vertical ocean dynamics physics.
113      !!---------------------------------------------------------------------
114      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
115      !!
116      IF( kt == nit000 .AND. lfirst )   CALL zdf_ctl_tam          ! initialisation & control of options
117
118      !                                          ! set time step
119      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restarting with Euler time stepping)
120      ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog)
121      ENDIF
122      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
123      !
124      CASE ( 0 )   ;   CALL dyn_zdf_exp_tan    ( kt, r2dt )      ! explicit scheme
125      CASE ( 1 )   ;   CALL dyn_zdf_imp_tan    ( kt, r2dt )      ! implicit scheme (k-j-i loop)
126      !
127      END SELECT
128   END SUBROUTINE dyn_zdf_tan
129
130   SUBROUTINE dyn_zdf_adj( kt )
131      !!----------------------------------------------------------------------
132      !!                  ***  ROUTINE dyn_zdf_adj  ***
133      !!
134      !! ** Purpose of the direct routine:
135      !!            compute the vertical ocean dynamics physics.
136      !!---------------------------------------------------------------------
137      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
138      !!
139      IF( kt == nitend .AND. lfirst )   CALL zdf_ctl_tam          ! initialisation & control of options
140
141      !                                          ! set time step
142      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restarting with Euler time stepping)
143      ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog)
144      ELSEIF( kt == nitend ) THEN ; r2dt = 2. * rdt
145      ENDIF
146      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
147      !
148      CASE ( 0 )   ;   CALL dyn_zdf_exp_adj    ( kt, r2dt )      ! explicit scheme
149      CASE ( 1 )   ;   CALL dyn_zdf_imp_adj    ( kt, r2dt )      ! implicit scheme (k-j-i loop)
150      !
151      END SELECT
152
153   END SUBROUTINE dyn_zdf_adj
154   SUBROUTINE zdf_ctl_tam
155      !!----------------------------------------------------------------------
156      !!                 ***  ROUTINE zdf_ctl_tam  ***
157      !!
158      !! ** Purpose :   initializations of the vertical diffusion scheme
159      !!
160      !! ** Method  :   implicit (euler backward) scheme (default)
161      !!                explicit (time-splitting) scheme if ln_zdfexp=T
162      !!----------------------------------------------------------------------
163      USE zdftke, ONLY : lk_zdftke
164      USE zdfkpp, ONLY : lk_zdfkpp
165      !!----------------------------------------------------------------------
166
167      ! Choice from ln_zdfexp read in namelist in zdfini
168      IF( ln_zdfexp ) THEN   ;   nzdf = 0           ! use explicit scheme
169      ELSE                   ;   nzdf = 1           ! use implicit scheme
170      ENDIF
171
172      ! Force implicit schemes
173      IF( lk_zdftke .OR. lk_zdfkpp   )   nzdf = 1   ! TKE or KPP physics
174      IF( ln_dynldf_iso              )   nzdf = 1   ! iso-neutral lateral physics
175      IF( ln_dynldf_hor .AND. ln_sco )   nzdf = 1   ! horizontal lateral physics in s-coordinate
176
177
178      IF( lk_esopa )    nzdf = -1                   ! Esopa key: All schemes used
179
180      IF(lwp) THEN                                  ! Print the choice
181         WRITE(numout,*)
182         WRITE(numout,*) 'dyn:zdf_ctl_tam : vertical dynamics physics scheme'
183         WRITE(numout,*) '~~~~~~~~~~~~~~~'
184         IF( nzdf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used'
185         IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme'
186         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme'
187      ENDIF
188      !
189      lfirst = .FALSE. 
190   END SUBROUTINE zdf_ctl_tam
191   SUBROUTINE dyn_zdf_adj_tst( kumadt )
192      !!-----------------------------------------------------------------------
193      !!
194      !!                  ***  ROUTINE dyn_zdf_adj_tst ***
195      !!
196      !! ** Purpose : Test the adjoint routine.
197      !!
198      !! ** Method  : Verify the scalar product
199      !!           
200      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
201      !!
202      !!              where  L   = tangent routine
203      !!                     L^T = adjoint routine
204      !!                     W   = diagonal matrix of scale factors
205      !!                     dx  = input perturbation (random field)
206      !!                     dy  = L dx
207      !!
208      !! ** Action  : Separate tests are applied for the following dx and dy:
209      !!               
210      !!              1) dx = ( SSH ) and dy = ( SSH )
211      !!                   
212      !! History :
213      !!        ! 08-08 (A. Vidard)
214      !!-----------------------------------------------------------------------
215      !! * Modules used
216
217      !! * Arguments
218      INTEGER, INTENT(IN) :: &
219         & kumadt          ! Output unit
220      INTEGER :: &
221         & ji,    &        ! dummy loop indices
222         & jj,    &       
223         & jk     
224      INTEGER, DIMENSION(jpi,jpj) :: &
225         & iseed_2d        ! 2D seed for the random number generator
226
227      !! * Local declarations
228       REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
229         & zub_tlin,     & ! Tangent input: before u-velocity
230         & zvb_tlin,     & ! Tangent input: before v-velocity
231         & zua_tlin,     & ! Tangent input: after u-velocity
232         & zva_tlin,     & ! Tangent input: after v-velocity
233         & zua_tlout,    & ! Tangent output: after u-velocity
234         & zva_tlout,    & ! Tangent output: after v-velocity
235         & zua_adin,     & ! Adjoint input: after u-velocity
236         & zva_adin,     & ! Adjoint input: after v-velocity
237         & zub_adout,    & ! Adjoint output: before u-velocity
238         & zvb_adout,    & ! Adjoint output: before v-velocity
239         & zua_adout,    & ! Adjoint output: after u-velocity
240         & zva_adout,    & ! Adjoint output: after v-velocity
241         & zau,          & ! 3D random field for ua
242         & zav,          & ! 3D random field for va
243         & zbu,          & ! 3D random field for ub
244         & zbv             ! 3D random field for vb
245      REAL(KIND=wp) :: &
246         & zsp1,         & ! scalar product involving the tangent routine
247         & zsp1_1,       & !   scalar product components
248         & zsp1_2,       & 
249         & zsp2,         & ! scalar product involving the adjoint routine
250         & zsp2_1,       & !   scalar product components
251         & zsp2_2,       & 
252         & zsp2_3,       & 
253         & zsp2_4
254      CHARACTER(LEN=14) :: cl_name
255      ! Allocate memory
256
257      ALLOCATE( &
258         & zua_tlin(jpi,jpj,jpk),     & 
259         & zva_tlin(jpi,jpj,jpk),     & 
260         & zub_tlin(jpi,jpj,jpk),     & 
261         & zvb_tlin(jpi,jpj,jpk),     & 
262         & zua_tlout(jpi,jpj,jpk),    & 
263         & zva_tlout(jpi,jpj,jpk),    & 
264         & zua_adin(jpi,jpj,jpk),     & 
265         & zva_adin(jpi,jpj,jpk),     & 
266         & zua_adout(jpi,jpj,jpk),    & 
267         & zva_adout(jpi,jpj,jpk),    & 
268         & zub_adout(jpi,jpj,jpk),    & 
269         & zvb_adout(jpi,jpj,jpk),    & 
270         & zau(jpi,jpj,jpk),          & 
271         & zav(jpi,jpj,jpk),          & 
272         & zbu(jpi,jpj,jpk),          & 
273         & zbv(jpi,jpj,jpk)           &
274         & )
275
276      ! Initialize the direct trajectory
277      avmu(:,:,:) = avm0 * umask(:,:,:)
278      avmv(:,:,:) = avm0 * vmask(:,:,:)
279
280      !==================================================================
281      ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
282      !    dy = ( hdivb_tl, hdivn_tl )
283      !==================================================================
284
285      !--------------------------------------------------------------------
286      ! Reset the tangent and adjoint variables
287      !--------------------------------------------------------------------
288     
289      zua_tlin(:,:,:)  = 0.0_wp     
290      zva_tlin(:,:,:)  = 0.0_wp     
291      zub_tlin(:,:,:)  = 0.0_wp     
292      zvb_tlin(:,:,:)  = 0.0_wp     
293      zua_tlout(:,:,:) = 0.0_wp     
294      zva_tlout(:,:,:) = 0.0_wp     
295      zua_adin(:,:,:)  = 0.0_wp     
296      zva_adin(:,:,:)  = 0.0_wp     
297      zua_adout(:,:,:) = 0.0_wp     
298      zva_adout(:,:,:) = 0.0_wp     
299      zub_adout(:,:,:) = 0.0_wp     
300      zvb_adout(:,:,:) = 0.0_wp     
301      zau(:,:,:)       = 0.0_wp           
302      zav(:,:,:)       = 0.0_wp           
303      zbu(:,:,:)       = 0.0_wp           
304      zbv(:,:,:)       = 0.0_wp           
305
306      ub_tl(:,:,:)     = 0.0_wp
307      vb_tl(:,:,:)     = 0.0_wp
308      ua_tl(:,:,:)     = 0.0_wp
309      va_tl(:,:,:)     = 0.0_wp
310      ub_ad(:,:,:)     = 0.0_wp
311      vb_ad(:,:,:)     = 0.0_wp
312      ua_ad(:,:,:)     = 0.0_wp
313      va_ad(:,:,:)     = 0.0_wp
314         
315      !--------------------------------------------------------------------
316      ! Initialize the tangent input with random noise: dx
317      !--------------------------------------------------------------------
318
319      DO jj = 1, jpj
320         DO ji = 1, jpi
321            iseed_2d(ji,jj) = - ( 596035 + &
322               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
323         END DO
324      END DO
325      CALL grid_random( iseed_2d, zbu, 'U', 0.0_wp, stdu )
326
327      DO jj = 1, jpj
328         DO ji = 1, jpi
329            iseed_2d(ji,jj) = - ( 523432 + &
330               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
331         END DO
332      END DO
333      CALL grid_random( iseed_2d, zbv, 'V', 0.0_wp, stdv )
334
335      DO jj = 1, jpj
336         DO ji = 1, jpi
337            iseed_2d(ji,jj) = - ( 432545 + &
338               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
339         END DO
340      END DO
341      CALL grid_random( iseed_2d, zau, 'U', 0.0_wp, stdu )
342
343      DO jj = 1, jpj
344         DO ji = 1, jpi
345            iseed_2d(ji,jj) = - ( 287503 + &
346               &                  mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
347         END DO
348      END DO
349      CALL grid_random( iseed_2d, zav, 'V', 0.0_wp, stdv )
350      DO jk = 1, jpk
351         DO jj = nldj, nlej
352            DO ji = nldi, nlei
353               zub_tlin(ji,jj,jk) = zbu(ji,jj,jk) 
354               zvb_tlin(ji,jj,jk) = zbv(ji,jj,jk) 
355               zua_tlin(ji,jj,jk) = zau(ji,jj,jk) 
356               zva_tlin(ji,jj,jk) = zav(ji,jj,jk) 
357            END DO
358         END DO
359      END DO
360
361      ub_tl(:,:,:) = zub_tlin(:,:,:)
362      vb_tl(:,:,:) = zvb_tlin(:,:,:)
363      ua_tl(:,:,:) = zua_tlin(:,:,:)
364      va_tl(:,:,:) = zva_tlin(:,:,:)
365
366      CALL dyn_zdf_tan( nit000 )
367
368      zua_tlout(:,:,:) = ua_tl(:,:,:)
369      zva_tlout(:,:,:) = va_tl(:,:,:)
370
371      !--------------------------------------------------------------------
372      ! Initialize the adjoint variables: dy^* = W dy
373      !--------------------------------------------------------------------
374
375      DO jk = 1, jpk
376        DO jj = nldj, nlej
377           DO ji = nldi, nlei
378              zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) &
379                 &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
380                 &               * umask(ji,jj,jk)
381              zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) &
382                 &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
383                 &               * vmask(ji,jj,jk)
384            END DO
385         END DO
386      END DO
387      !--------------------------------------------------------------------
388      ! Compute the scalar product: ( L dx )^T W dy
389      !--------------------------------------------------------------------
390
391      zsp1_1 = DOT_PRODUCT( zua_tlout, zua_adin )
392      zsp1_2 = DOT_PRODUCT( zva_tlout, zva_adin )
393      zsp1   = zsp1_1 + zsp1_2
394
395      !--------------------------------------------------------------------
396      ! Call the adjoint routine: dx^* = L^T dy^*
397      !--------------------------------------------------------------------
398
399      ua_ad(:,:,:) = zua_adin(:,:,:)
400      va_ad(:,:,:) = zva_adin(:,:,:)
401
402      CALL dyn_zdf_adj ( nit000 )
403      zub_adout(:,:,:)   = ub_ad(:,:,:)
404      zvb_adout(:,:,:)   = vb_ad(:,:,:)
405      zua_adout(:,:,:)   = ua_ad(:,:,:)
406      zva_adout(:,:,:)   = va_ad(:,:,:)
407     
408      zsp2_1 = DOT_PRODUCT( zub_tlin, zub_adout )
409      zsp2_2 = DOT_PRODUCT( zvb_tlin, zvb_adout )
410      zsp2_3 = DOT_PRODUCT( zua_tlin, zua_adout )
411      zsp2_4 = DOT_PRODUCT( zva_tlin, zva_adout )
412      zsp2   = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 
413
414      ! Compare the scalar products
415
416      ! 14 char:'12345678901234'
417      cl_name = 'dyn_zdf_adj   '
418      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
419
420      DEALLOCATE( &
421         & zua_tlin,     & 
422         & zva_tlin,     & 
423         & zub_tlin,     & 
424         & zvb_tlin,     & 
425         & zua_tlout,    & 
426         & zva_tlout,    & 
427         & zua_adin,     & 
428         & zva_adin,     & 
429         & zua_adout,    & 
430         & zva_adout,    & 
431         & zub_adout,    & 
432         & zvb_adout,    & 
433         & zau,          & 
434         & zav,          & 
435         & zbu,          & 
436         & zbv           &
437         & )
438
439  END SUBROUTINE dyn_zdf_adj_tst
440   !!==============================================================================
441#endif
442END MODULE dynzdf_tam
Note: See TracBrowser for help on using the repository browser.