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.
dynbfr_tam.F90 in branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/DYN/dynbfr_tam.F90 @ 3611

Last change on this file since 3611 was 3611, checked in by pabouttier, 11 years ago

Add TAM code and ORCA2_TAM configuration - see Ticket #1007

  • Property svn:executable set to *
File size: 15.5 KB
Line 
1MODULE dynbfr_tam
2#ifdef key_tam
3   !!==============================================================================
4   !!                 ***  MODULE  dynbfr  ***
5   !! Ocean dynamics :  bottom friction component of the momentum mixing trend
6   !!==============================================================================
7   !! History of the drect module:
8   !!            9.0  !  2008-11  (A. C. Coward)  Original code
9   !! History of the TAM module:
10   !!  NEMO      3.2  ! 2010-04 (F. Vigilant) Original code
11   !!  NEMO      3.4  ! 2012-07 (P.-A. bouttier) Phasing with 3.4
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   dyn_bfr      : Update the momentum trend with the bottom friction contribution
16   !!----------------------------------------------------------------------
17   USE oce
18   USE oce_tam
19   USE par_oce
20   USE dom_oce
21   USE zdf_oce         ! ocean vertical physics variables
22   USE zdfbfr
23   USE zdf_oce_tam
24   USE in_out_manager
25   USE gridrandom
26   USE dotprodfld
27   USE tstool_tam
28   USE timing
29   USE wrk_nemo
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   dyn_bfr_tan     ! routine called by step_tam.F90
35   PUBLIC   dyn_bfr_adj     ! routine called by step_tam.F90
36   PUBLIC   dyn_bfr_adj_tst ! routine called by the tst.F90
37
38   !! * Substitutions
39#  include "domzgr_substitute.h90"
40#  include "zdfddm_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42
43CONTAINS
44
45   SUBROUTINE dyn_bfr_tan( kt )
46      !!----------------------------------------------------------------------
47      !!                  ***  ROUTINE dyn_bfr_tan  ***
48      !!
49      !! ** Purpose of direct routine:
50      !!                compute the bottom friction ocean dynamics physics.
51      !!
52      !! ** Action  :   (ua,va)   momentum trend increased by bottom friction trend
53      !!---------------------------------------------------------------------
54      !!
55      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
56      !!
57      INTEGER  ::   ji, jj          ! dummy loop indexes
58      INTEGER  ::   ikbu  , ikbv    ! temporary integers
59      REAL(wp) ::   zm1_2dt, zbfru, zbfrv   ! temporary scalar
60      REAL(wp) ::   zbfrutl, zbfrvtl     ! temporary scalar
61      !!---------------------------------------------------------------------
62      !
63      !
64      IF( nn_timing == 1 )  CALL timing_start('dyn_bfr_tan')
65      !
66      IF ( .NOT. ln_bfrimp ) THEN
67         zm1_2dt = -1._wp / ( 2._wp * rdt )
68# if defined key_vectopt_loop
69         DO jj = 1, 1
70            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
71# else
72         DO jj = 2, jpjm1
73            DO ji = 2, jpim1
74# endif
75               ikbu = mbku(ji,jj)
76               ikbv = mbkv(ji,jj)
77               !
78               ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr)
79               IF ( bfrua(ji,jj) >= fse3u(ji,jj,ikbu)*zm1_2dt ) THEN
80                  zbfru   = bfrua(   ji,jj)
81                  zbfrutl  = bfrua_tl(ji,jj)
82               ELSE
83                  zbfru   = fse3u(ji,jj,ikbu)*zm1_2dt
84                  zbfrutl  = 0.0_wp
85               END IF
86
87               IF ( bfrva(ji,jj) >= fse3v(ji,jj,ikbv)*zm1_2dt ) THEN
88                  zbfrv   = bfrva(   ji,jj)
89                  zbfrvtl = bfrva_tl(ji,jj)
90               ELSE
91                  zbfrv   = fse3v(ji,jj,ikbv)*zm1_2dt
92                  zbfrvtl = 0.0_wp
93               END IF
94               !
95               ua_tl(ji,jj,ikbu) = ua_tl(ji,jj,ikbu) &
96                  &                + ( zbfru * ub_tl(ji,jj,ikbu) + zbfrutl * ub(ji,jj,ikbu) ) &
97                  &                / fse3u(ji,jj,ikbu)
98               va_tl(ji,jj,ikbv) = va_tl(ji,jj,ikbv) &
99                  &                + ( zbfrv * vb_tl(ji,jj,ikbv) + zbfrvtl * vb(ji,jj,ikbv) ) &
100                  &                / fse3v(ji,jj,ikbv)
101               !
102            END DO
103         END DO
104      ENDIF
105      !
106      IF( nn_timing == 1 )  CALL timing_stop('dyn_bfr_tan')
107      !
108   END SUBROUTINE dyn_bfr_tan
109
110   SUBROUTINE dyn_bfr_adj( kt )
111      !!----------------------------------------------------------------------
112      !!                  ***  ROUTINE dyn_bfr_adj  ***
113      !!
114      !! ** Purpose of direct routine:
115      !!                compute the bottom friction ocean dynamics physics.
116      !!
117      !! ** Action  :   (ua,va)   momentum trend increased by bottom friction trend
118      !!---------------------------------------------------------------------
119      !!
120      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
121      !!
122      INTEGER  ::   ji, jj          ! dummy loop indexes
123      INTEGER  ::   ikbu  , ikbv    ! temporary integers
124      REAL(wp) ::   zm1_2dt, zbfru, zbfrv   ! temporary scalar
125      REAL(wp) ::   zbfruad, zbfrvad     ! temporary scalar
126      !!---------------------------------------------------------------------
127      !
128      IF( nn_timing == 1 )  CALL timing_start('dyn_bfr_adj')
129      !
130      IF( .NOT.ln_bfrimp) THEN
131         zm1_2dt = -1._wp / ( 2._wp * rdt )
132         zbfruad = 0.0_wp  ;  zbfrvad = 0.0_wp
133# if defined key_vectopt_loop
134         DO jj = 1, 1
135            DO ji = jpij-jpj-1, jpi+2, -1   ! vector opt. (forced unrolling)
136# else
137         DO jj = jpjm1, 2, -1
138            DO ji = jpim1, 2, -1
139# endif
140               ikbu = mbku(ji,jj)
141               ikbv = mbkv(ji,jj)
142               !
143               ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr)
144               zbfru = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbu)*zm1_2dt )
145               zbfrv = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbv)*zm1_2dt )
146               !
147               ub_ad(ji,jj,ikbu) = ub_ad(ji,jj,ikbu) + zbfru * ua_ad(ji,jj,ikbu) / fse3u(ji,jj,ikbu)
148               zbfruad = zbfruad + ub(ji,jj,ikbu) * ua_ad(ji,jj,ikbu) / fse3u(ji,jj,ikbu)
149               vb_ad(ji,jj,ikbv) = vb_ad(ji,jj,ikbv) + zbfrv * va_ad(ji,jj,ikbv) / fse3v(ji,jj,ikbv)
150               zbfrvad = zbfrvad + vb(ji,jj,ikbv) * va_ad(ji,jj,ikbv) / fse3v(ji,jj,ikbv)
151               !
152               ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr)
153               IF ( bfrua(ji,jj) >= fse3u(ji,jj,ikbu)*zm1_2dt ) THEN
154                  bfrua_ad(ji,jj) = bfrua_ad(ji,jj) + zbfruad
155                  zbfruad  = 0.0_wp
156               ELSE
157                  zbfruad  = 0.0_wp
158               END IF
159               !
160               IF ( bfrva(ji,jj) >= fse3v(ji,jj,ikbv)*zm1_2dt ) THEN
161                  bfrva_ad(ji,jj) = bfrva_ad(ji,jj) + zbfrvad
162                  zbfrvad = 0.0_wp
163               ELSE
164                  zbfrvad = 0.0_wp
165               END IF
166               !
167            END DO
168         END DO
169      ENDIF
170      !
171      IF( nn_timing == 1 )  CALL timing_stop('dyn_bfr_adj')
172      !
173   END SUBROUTINE dyn_bfr_adj
174
175   SUBROUTINE dyn_bfr_adj_tst( kumadt )
176      !!-----------------------------------------------------------------------
177      !!
178      !!                  ***  ROUTINE dyn_bfr_adj_tst ***
179      !!
180      !! ** Purpose : Test the adjoint routine.
181      !!
182      !! ** Method  : Verify the scalar product
183      !!
184      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
185      !!
186      !!              where  L   = tangent routine
187      !!                     L^T = adjoint routine
188      !!                     W   = diagonal matrix of scale factors
189      !!                     dx  = input perturbation (random field)
190      !!                     dy  = L dx
191      !!
192      !! ** Action  :
193      !!
194      !! History :
195      !!        ! 2010-04 (F. Vigilant)
196      !!-----------------------------------------------------------------------
197      !! * Modules used
198
199      !! * Arguments
200      INTEGER, INTENT(IN) :: &
201         & kumadt        ! Output unit
202
203      !! * Local declarations
204      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
205         & zua_tlin,    & ! Tangent input: ua_tl
206         & zva_tlin,    & ! Tangent input: va_tl
207         & zub_tlin,    & ! Tangent input: ub_tl
208         & zvb_tlin,    & ! Tangent input: vb_tl
209         & zua_tlout,   & ! Tangent output: ua_tl
210         & zva_tlout,   & ! Tangent output: va_tl
211         & zua_adin,    & ! Adjoint input: ua_ad
212         & zva_adin,    & ! Adjoint input: va_ad
213         & zua_adout,   & ! Adjoint output: ua_ad
214         & zva_adout,   & ! Adjoint output: va_ad
215         & zub_adout,   & ! Adjoint oputput: ub_ad
216         & zvb_adout,   & ! Adjoint output: vb_ad
217         & znu            ! 3D random field for u
218      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
219         & zbu_tlin,     & ! Tangent input: bfrua_tl
220         & zbv_tlin,     & ! Tangent input: bfrva_tl
221         & zbu_adout,    & ! Adjoint output: bfrua_ad
222         & zbv_adout       ! Adjoint output: bfrva_ad
223      REAL(wp) :: &
224         & zsp1,    &   ! scalar product involving the tangent routine
225         & zsp2         ! scalar product involving the adjoint routine
226      INTEGER, DIMENSION(jpi,jpj) :: &
227         & iseed_2d    ! 2D seed for the random number generator
228      INTEGER :: &
229         & ji, &
230         & jj, &
231         & jk
232      CHARACTER (LEN=14) :: &
233         & cl_name
234
235      ! Allocate memory
236
237      ALLOCATE( &
238         & zua_tlin(jpi,jpj,jpk),  &
239         & zva_tlin(jpi,jpj,jpk),  &
240         & zub_tlin(jpi,jpj,jpk),  &
241         & zvb_tlin(jpi,jpj,jpk),  &
242         & zua_tlout(jpi,jpj,jpk), &
243         & zva_tlout(jpi,jpj,jpk), &
244         & zua_adin(jpi,jpj,jpk),  &
245         & zva_adin(jpi,jpj,jpk),  &
246         & zua_adout(jpi,jpj,jpk), &
247         & zva_adout(jpi,jpj,jpk), &
248         & zub_adout(jpi,jpj,jpk), &
249         & zvb_adout(jpi,jpj,jpk), &
250         & znu(jpi,jpj,jpk)        &
251         & )
252      ALLOCATE( &
253         & zbu_tlin(jpi,jpj),      &
254         & zbv_tlin(jpi,jpj),      &
255         & zbu_adout(jpi,jpj),     &
256         & zbv_adout(jpi,jpj)      &
257         & )
258
259      !=========================================================================
260      !     dx = ( ub_tl, ua_tl, vb_tl, va_tl )
261      ! and dy = ( ua_tl, va_tl )
262      !=========================================================================
263
264      !--------------------------------------------------------------------
265      ! Reset the tangent and adjoint variables
266      !--------------------------------------------------------------------
267
268      zub_tlin (:,:,:) = 0.0_wp
269      zvb_tlin (:,:,:) = 0.0_wp
270      zua_tlin (:,:,:) = 0.0_wp
271      zva_tlin (:,:,:) = 0.0_wp
272      zua_tlout(:,:,:) = 0.0_wp
273      zva_tlout(:,:,:) = 0.0_wp
274      zua_adin (:,:,:) = 0.0_wp
275      zva_adin (:,:,:) = 0.0_wp
276      zub_adout(:,:,:) = 0.0_wp
277      zvb_adout(:,:,:) = 0.0_wp
278      zua_adout(:,:,:) = 0.0_wp
279      zva_adout(:,:,:) = 0.0_wp
280
281      ub_tl(:,:,:) = 0.0_wp
282      vb_tl(:,:,:) = 0.0_wp
283      ua_tl(:,:,:) = 0.0_wp
284      va_tl(:,:,:) = 0.0_wp
285      ub_ad(:,:,:) = 0.0_wp
286      vb_ad(:,:,:) = 0.0_wp
287      ua_ad(:,:,:) = 0.0_wp
288      va_ad(:,:,:) = 0.0_wp
289      !--------------------------------------------------------------------
290      ! Initialize the tangent input with random noise: dx
291      !--------------------------------------------------------------------
292
293      CALL grid_random(  znu, 'U', 0.0_wp, stdu )
294
295      DO jk = 1, jpk
296         DO jj = nldj, nlej
297            DO ji = nldi, nlei
298               zua_tlin(ji,jj,jk) = znu(ji,jj,jk)
299            END DO
300         END DO
301      END DO
302
303      CALL grid_random(  znu, 'V', 0.0_wp, stdv )
304
305      DO jk = 1, jpk
306         DO jj = nldj, nlej
307            DO ji = nldi, nlei
308              zva_tlin(ji,jj,jk) = znu(ji,jj,jk)
309            END DO
310         END DO
311      END DO
312      CALL grid_random(  znu, 'U', 0.0_wp, stdu )
313
314      DO jk = 1, jpk
315         DO jj = nldj, nlej
316            DO ji = nldi, nlei
317             zub_tlin(ji,jj,jk) = znu(ji,jj,jk)
318            END DO
319         END DO
320      END DO
321      CALL grid_random(  znu, 'V', 0.0_wp, stdv )
322
323      DO jk = 1, jpk
324         DO jj = nldj, nlej
325            DO ji = nldi, nlei
326             zvb_tlin(ji,jj,jk) = znu(ji,jj,jk)
327            END DO
328         END DO
329      END DO
330      CALL grid_random(  znu, 'U', 0.0_wp, stdv )
331
332      DO jj = nldj, nlej
333         DO ji = nldi, nlei
334            zbu_tlin(ji,jj) = znu(ji,jj,1)
335         END DO
336      END DO
337      CALL grid_random(  znu, 'V', 0.0_wp, stdv )
338
339      DO jj = nldj, nlej
340         DO ji = nldi, nlei
341            zbv_tlin(ji,jj) = znu(ji,jj,1)
342         END DO
343      END DO
344
345      !--------------------------------------------------------------------
346      ! Call the tangent routine: dy = L dx
347      !--------------------------------------------------------------------
348
349      ua_tl(:,:,:) = zua_tlin(:,:,:)
350      va_tl(:,:,:) = zva_tlin(:,:,:)
351      ub_tl(:,:,:) = zub_tlin(:,:,:)
352      vb_tl(:,:,:) = zvb_tlin(:,:,:)
353      bfrua_tl(:,:) = zbu_tlin(:,:)
354      bfrva_tl(:,:) = zbv_tlin(:,:)
355
356      CALL dyn_bfr_tan( nit000 )
357
358      zua_tlout(:,:,:) = ua_tl(:,:,:)
359      zva_tlout(:,:,:) = va_tl(:,:,:)
360      !--------------------------------------------------------------------
361      ! Initialize the adjoint variables: dy^* = W dy
362      !--------------------------------------------------------------------
363
364      DO jk = 1, jpk
365        DO jj = nldj, nlej
366           DO ji = nldi, nlei
367              zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) &
368                 &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
369                 &               * umask(ji,jj,jk)
370              zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) &
371                 &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
372                 &               * vmask(ji,jj,jk)
373            END DO
374         END DO
375      END DO
376
377      !--------------------------------------------------------------------
378      ! Compute the scalar product: ( L dx )^T W dy
379      !--------------------------------------------------------------------
380
381      zsp1 =   DOT_PRODUCT( zua_tlout  , zua_adin   ) + DOT_PRODUCT( zva_tlout  , zva_adin   )
382
383
384      !--------------------------------------------------------------------
385      ! Call the adjoint routine: dx^* = L^T dy^*
386      !--------------------------------------------------------------------
387
388      ua_ad(:,:,:) = zua_adin(:,:,:)
389      va_ad(:,:,:) = zva_adin(:,:,:)
390
391      CALL dyn_bfr_adj( nit000 )
392
393      zua_adout(:,:,:) = ua_ad(:,:,:)
394      zva_adout(:,:,:) = va_ad(:,:,:)
395      zub_adout(:,:,:) = ub_ad(:,:,:)
396      zvb_adout(:,:,:) = vb_ad(:,:,:)
397      zbu_adout(:,:) = bfrua_ad(:,:)
398      zbv_adout(:,:) = bfrva_ad(:,:)
399
400      !--------------------------------------------------------------------
401      ! Compute the scalar product: dx^T L^T W dy
402      !--------------------------------------------------------------------
403
404      zsp2 =   DOT_PRODUCT( zua_tlin  , zua_adout   ) &
405            &   + DOT_PRODUCT( zva_tlin  , zva_adout   ) &
406            &   + DOT_PRODUCT( zub_tlin  , zub_adout   ) &
407            &   + DOT_PRODUCT( zvb_tlin  , zvb_adout   ) &
408            &   + DOT_PRODUCT( zbu_tlin  , zbu_adout   ) &
409            &   + DOT_PRODUCT( zbv_tlin  , zbv_adout   )
410
411      ! Compare the scalar products
412
413      ! 14 char:'12345678901234'
414      cl_name = 'dyn_bfr_adj   '
415      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
416
417      ! Deallocate memory
418
419      DEALLOCATE( &
420         & zua_tlin,  &
421         & zva_tlin,  &
422         & zub_tlin,  &
423         & zvb_tlin,  &
424         & zua_tlout, &
425         & zva_tlout, &
426         & zua_adin,  &
427         & zva_adin,  &
428         & zua_adout, &
429         & zva_adout, &
430         & zub_adout, &
431         & zvb_adout, &
432         & znu        &
433         & )
434      DEALLOCATE( &
435         & zbu_tlin,      &
436         & zbv_tlin,      &
437         & zbu_adout,     &
438         & zbv_adout      &
439         & )
440
441   END SUBROUTINE dyn_bfr_adj_tst
442   !!==============================================================================
443#endif
444END MODULE dynbfr_tam
Note: See TracBrowser for help on using the repository browser.