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

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/ZDF/zdfbfr_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: 18.3 KB
Line 
1MODULE zdfbfr_tam
2   !!======================================================================
3   !!                       ***  MODULE  zdfbfr_tam  ***
4   !! Ocean physics: Bottom friction
5   !!======================================================================
6   !! History of the direct module:
7   !!            OPA  ! 1997-06  (G. Madec, A.-M. Treguier)  Original code
8   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
9   !!            3.2  ! 2009-09  (A.C.Coward)  Correction to include barotropic contribution
10   !! History of the T&A module:
11   !!   NEMO     3.2.2! 2011-02  (A. Vidard) Original version
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   zdf_bfr_tan  : update momentum Kz at the ocean bottom due to the type of bottom friction chosen
16   !!   zdf_bfr_adj  : update momentum Kz at the ocean bottom due to the type of bottom friction chosen
17   !!                  parameters.
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and tracers variables
20   USE dom_oce         ! ocean space and time domain variables
21   USE zdf_oce         ! ocean vertical physics variables
22   USE in_out_manager  ! I/O manager
23   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
24   USE lib_mpp         ! distributed memory computing
25   USE zdfbfr
26   USE oce_tam
27   USE lbclnk_tam
28   USE timing
29   USE zdf_oce_tam
30   USE gridrandom
31   USE dotprodfld
32   USE paresp
33   USE tstool_tam
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   zdf_bfr_tan    ! called by step_tam.F90
39   PUBLIC   zdf_bfr_adj    ! called by step_tam.F90
40   PUBLIC   zdf_bfr_adj_tst
41   PUBLIC   zdf_bfr_init_tam
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
48   !! $Id$
49   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51
52CONTAINS
53   SUBROUTINE zdf_bfr_init_tam
54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE zdf_bfr_init  ***
56      !!
57      !! ** Purpose :   Initialization of the bottom friction
58      !!
59      !! ** Method  :   Read the nammbf namelist and check their consistency
60      !!              called at the first timestep (nit000)
61      !!----------------------------------------------------------------------
62      !!----------------------------------------------------------------------
63      !
64      IF( nn_timing == 1 )  CALL timing_start('zdf_bfr_init_tam')
65      !
66      bfrua_tl = 0._wp
67      bfrva_tl = 0._wp
68      bfrua_ad = 0._wp
69      bfrva_ad = 0._wp
70
71      IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr_init_tam')
72      !
73   END SUBROUTINE zdf_bfr_init_tam
74   SUBROUTINE zdf_bfr_tan( kt )
75      !!----------------------------------------------------------------------
76      !!                   ***  ROUTINE zdf_bfr_tan  ***
77      !!
78      !! ** Purpose :   tangent of the computation of the bottom friction coefficient.
79      !!
80      !! ** Method  :   Calculate and store part of the momentum trend due
81      !!              to bottom friction following the chosen friction type
82      !!              (free-slip, linear, or quadratic). The component
83      !!              calculated here is multiplied by the bottom velocity in
84      !!              dyn_bfr to provide the trend term.
85      !!                The coefficients are updated at each time step only
86      !!              in the quadratic case.
87      !!
88      !! ** Action  :   bfrua , bfrva   bottom friction coefficients
89      !!----------------------------------------------------------------------
90      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
91      !!
92      INTEGER  ::   ji, jj         ! dummy loop indices
93      INTEGER  ::   ikbu           ! temporary integers
94      INTEGER  ::   ikbv           !    -          -
95      REAL(wp) ::   zvu, zuv, zecu, zecv   ! temporary scalars
96      REAL(wp) ::   zvutl, zuvtl, zecutl, zecvtl   ! temporary scalars
97      !!----------------------------------------------------------------------
98      !
99      IF( nn_timing == 1 )  CALL timing_start('zdf_bfr_tan')
100      !
101      IF( kt == nit000 )   THEN
102         bfrua_tl = 0.0_wp
103         bfrva_tl = 0.0_wp
104      END IF
105      IF( nn_bfr == 2 ) THEN                 ! quadratic botton friction
106         ! Calculate and store the quadratic bottom friction coefficient bfrua and bfrva
107         ! where bfrUa = C_d*SQRT(u_bot^2 + v_bot^2 + e_b) {U=[u,v]}
108         ! from these the trend due to bottom friction:  -F_h/e3U  can be calculated
109         ! where -F_h/e3U_bot = bfrUa*Ub/e3U_bot {U=[u,v]}
110         !
111# if defined key_vectopt_loop
112         DO jj = 1, 1
113!CDIR NOVERRCHK
114            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
115# else
116!CDIR NOVERRCHK
117         DO jj = 2, jpjm1
118!CDIR NOVERRCHK
119            DO ji = 2, jpim1
120# endif
121               ikbu   = mbku(ji,jj)
122               ikbv   = mbkv(ji,jj)
123               !
124               zvu   = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     &
125                  &            + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  )
126               zuv   = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     &
127                  &            + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  )
128               zvutl = 0.25 * (  vn_tl(ji,jj  ,ikbu) + vn_tl(ji+1,jj  ,ikbu)     &
129                  &           + vn_tl(ji,jj-1,ikbu) + vn_tl(ji+1,jj-1,ikbu)  )
130               zuvtl = 0.25 * (  un_tl(ji,jj  ,ikbv) + un_tl(ji-1,jj  ,ikbv)     &
131                  &           + un_tl(ji,jj+1,ikbv) + un_tl(ji-1,jj+1,ikbv)  )
132               !
133               zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2  )
134               zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2  )
135               zecutl = ( un(ji,jj,ikbu) * un_tl(ji,jj,ikbu) + zvu*zvutl )  &
136                  &     / SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2  )
137               zecvtl = ( vn(ji,jj,ikbv) * vn_tl(ji,jj,ikbv) + zuv*zuvtl )  &
138                  &     / SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2  )
139               !
140               bfrua_tl(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji+1,jj  ) ) * zecutl
141               bfrva_tl(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji  ,jj+1) ) * zecvtl
142            END DO
143         END DO
144         !
145         CALL lbc_lnk( bfrua_tl, 'U', 1. )   ;   CALL lbc_lnk( bfrva_tl, 'V', 1. ) ! Lateral boundary condition
146         !
147      ENDIF
148      !
149      IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr_tan')
150      !
151   END SUBROUTINE zdf_bfr_tan
152   SUBROUTINE zdf_bfr_adj( kt )
153      !!----------------------------------------------------------------------
154      !!                   ***  ROUTINE zdf_bfr_adj  ***
155      !!
156      !! ** Purpose :   adjoint of the computation of the bottom friction coefficient.
157      !!
158      !! ** Method  :   Calculate and store part of the momentum trend due
159      !!              to bottom friction following the chosen friction type
160      !!              (free-slip, linear, or quadratic). The component
161      !!              calculated here is multiplied by the bottom velocity in
162      !!              dyn_bfr to provide the trend term.
163      !!                The coefficients are updated at each time step only
164      !!              in the quadratic case.
165      !!
166      !! ** Action  :   bfrua , bfrva   bottom friction coefficients
167      !!----------------------------------------------------------------------
168      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
169      !!
170      INTEGER  ::   ji, jj         ! dummy loop indices
171      INTEGER  ::   ikbu, ikbum1   ! temporary integers
172      INTEGER  ::   ikbv, ikbvm1   !    -          -
173      REAL(wp) ::   zvu, zuv, zecu, zecv   ! temporary scalars
174      REAL(wp) ::   zvuad, zuvad, zecuad, zecvad   ! temporary scalars
175      !!----------------------------------------------------------------------
176      !
177      IF( nn_timing == 1 )  CALL timing_start('zdf_bfr_adj')
178      !
179      zvuad = 0.0_wp   ;   zuvad = 0.0_wp   ;   zecuad = 0.0_wp   ;   zecvad = 0.0_wp
180
181      IF( kt == nitend )   THEN
182         bfrua_ad = 0.0_wp
183         bfrva_ad = 0.0_wp
184      END IF
185      IF( nn_bfr == 2 ) THEN                 ! quadratic botton friction
186         ! Calculate and store the quadratic bottom friction coefficient bfrua and bfrva
187         ! where bfrUa = C_d*SQRT(u_bot^2 + v_bot^2 + e_b) {U=[u,v]}
188         ! from these the trend due to bottom friction:  -F_h/e3U  can be calculated
189         ! where -F_h/e3U_bot = bfrUa*Ub/e3U_bot {U=[u,v]}
190         !
191         CALL lbc_lnk_adj( bfrua_ad, 'U', 1. )   ;   CALL lbc_lnk_adj( bfrva_ad, 'V', 1. ) ! Lateral boundary condition
192         !
193# if defined key_vectopt_loop
194         DO jj = 1, 1
195!CDIR NOVERRCHK
196            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
197# else
198!CDIR NOVERRCHK
199         DO jj = 2, jpjm1
200!CDIR NOVERRCHK
201            DO ji = 2, jpim1
202# endif
203               ikbu   = mbku(ji,jj)
204               ikbv   = mbkv(ji,jj)
205               ! direct computation
206               zvu   = 0.25 * (  vn(ji,jj  ,ikbu) + vn(ji+1,jj  ,ikbu)     &
207                  &            + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu)  )
208               zuv   = 0.25 * (  un(ji,jj  ,ikbv) + un(ji-1,jj  ,ikbv)     &
209                  &            + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv)  )
210               zecu = SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2  )
211               zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2  )
212               ! Adjoint counterpart
213
214               zecuad = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji+1,jj  ) ) * bfrua_ad(ji,jj)
215               bfrua_ad(ji,jj) = 0.0_wp
216               zecvad = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji  ,jj+1) ) * bfrva_ad(ji,jj)
217               bfrva_ad(ji,jj) = 0.0_wp
218               !
219               un_ad(ji,jj,ikbu) = un_ad(ji,jj,ikbu) + zecuad * un(ji,jj,ikbu) &
220                  &                  / SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2  )
221               zvuad  = zecuad * zvu / SQRT(  un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2  )
222
223               vn_ad(ji,jj,ikbv) = vn_ad(ji,jj,ikbv) + zecvad * vn(ji,jj,ikbv) &
224                  &                  / SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2  )
225               zuvad  = zecvad * zuv / SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2  )
226               !
227               vn_ad(ji  ,jj  ,ikbu) = vn_ad(ji,jj    ,ikbu) + zvuad * 0.25
228               vn_ad(ji+1,jj  ,ikbu) = vn_ad(ji+1,jj  ,ikbu) + zvuad * 0.25
229               vn_ad(ji  ,jj-1,ikbu) = vn_ad(ji,jj-1  ,ikbu) + zvuad * 0.25
230               vn_ad(ji+1,jj-1,ikbu) = vn_ad(ji+1,jj-1,ikbu) + zvuad * 0.25
231
232               un_ad(ji  ,jj  ,ikbv) = un_ad(ji  ,jj  ,ikbv) + zuvad * 0.25
233               un_ad(ji-1,jj  ,ikbv) = un_ad(ji-1,jj  ,ikbv) + zuvad * 0.25
234               un_ad(ji  ,jj+1,ikbv) = un_ad(ji  ,jj+1,ikbv) + zuvad * 0.25
235               un_ad(ji-1,jj+1,ikbv) = un_ad(ji-1,jj+1,ikbv) + zuvad * 0.25
236               !
237            END DO
238         END DO
239         !
240      ENDIF
241      !
242      IF ( kt == nit000 ) THEN
243         bfrua_ad = 0.0_wp
244         bfrva_ad = 0.0_wp
245      END IF
246      !
247      IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr_adj')
248      !
249   END SUBROUTINE zdf_bfr_adj
250   SUBROUTINE zdf_bfr_adj_tst( kumadt )
251      !!-----------------------------------------------------------------------
252      !!
253      !!                  ***  ROUTINE zdf_bfr_adj_tst ***
254      !!
255      !! ** Purpose : Test the adjoint routine.
256      !!
257      !! ** Method  : Verify the scalar product
258      !!
259      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
260      !!
261      !!              where  L   = tangent routine
262      !!                     L^T = adjoint routine
263      !!                     W   = diagonal matrix of scale factors
264      !!                     dx  = input perturbation (random field)
265      !!                     dy  = L dx
266      !!
267      !! ** Action  :
268      !!
269      !! History :
270      !!        ! 09-01 (A. Weaver)
271      !!-----------------------------------------------------------------------
272      !! * Modules used
273
274      !! * Arguments
275      INTEGER, INTENT(IN) :: &
276         & kumadt        ! Output unit
277
278      !! * Local declarations
279      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
280         & zua_tlin,    & ! Tangent input: ua_tl
281         & zva_tlin,    & ! Tangent input: va_tl
282         & zua_tlout,   & ! Tangent output: ua_tl
283         & zva_tlout,   & ! Tangent output: va_tl
284         & zua_adin,    & ! Adjoint input: ua_ad
285         & zva_adin,    & ! Adjoint input: va_ad
286         & zua_adout,   & ! Adjoint output: ua_ad
287         & zva_adout,   & ! Adjoint output: va_ad
288         & znu            ! 3D random field for u
289
290      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
291         & zspgu_tlout, zspgv_tlout, zspgu_adin, zspgv_adin
292
293      REAL(wp) :: &
294         & zsp1,    &   ! scalar product involving the tangent routine
295         & zsp2         ! scalar product involving the adjoint routine
296      INTEGER :: &
297         & ji, &
298         & jj, &
299         & jk
300      CHARACTER (LEN=14) :: &
301         & cl_name
302
303      ALLOCATE( &
304         & zua_tlin(jpi,jpj,jpk),  &
305         & zva_tlin(jpi,jpj,jpk),  &
306         & zua_tlout(jpi,jpj,jpk), &
307         & zva_tlout(jpi,jpj,jpk), &
308         & zua_adin(jpi,jpj,jpk),  &
309         & zva_adin(jpi,jpj,jpk),  &
310         & zua_adout(jpi,jpj,jpk), &
311         & zva_adout(jpi,jpj,jpk), &
312         & znu(jpi,jpj,jpk)        &
313         & )
314
315      ALLOCATE( zspgu_tlout (jpi,jpj), zspgv_tlout (jpi,jpj), zspgu_adin (jpi,jpj), zspgv_adin (jpi,jpj))
316
317      !--------------------------------------------------------------------
318      ! Reset the tangent and adjoint variables
319      !--------------------------------------------------------------------
320
321      zua_tlin (:,:,:) = 0.0_wp
322      zva_tlin (:,:,:) = 0.0_wp
323      zua_tlout(:,:,:) = 0.0_wp
324      zva_tlout(:,:,:) = 0.0_wp
325      zua_adin (:,:,:) = 0.0_wp
326      zva_adin (:,:,:) = 0.0_wp
327      zua_adout(:,:,:) = 0.0_wp
328      zva_adout(:,:,:) = 0.0_wp
329
330      zspgu_adin (:,:) = 0.0_wp
331      zspgv_adin (:,:) = 0.0_wp
332      zspgu_tlout(:,:) = 0.0_wp
333      zspgv_tlout(:,:) = 0.0_wp
334
335      ua_tl(:,:,:) = 0.0_wp
336      va_tl(:,:,:) = 0.0_wp
337      spgu_tl(:,:) = 0.0_wp
338      spgv_tl(:,:) = 0.0_wp
339      ua_ad(:,:,:) = 0.0_wp
340      va_ad(:,:,:) = 0.0_wp
341      spgu_ad(:,:) = 0.0_wp
342      spgv_ad(:,:) = 0.0_wp
343      !--------------------------------------------------------------------
344      ! Initialize the tangent input with random noise: dx
345      !--------------------------------------------------------------------
346      CALL grid_random(  znu, 'U', 0.0_wp, stdu )
347
348      DO jk = 1, jpk
349         DO jj = nldj, nlej
350            DO ji = nldi, nlei
351               zua_tlin(ji,jj,jk) = znu(ji,jj,jk)
352            END DO
353         END DO
354      END DO
355
356      CALL grid_random(  znu, 'V', 0.0_wp, stdv )
357
358      DO jk = 1, jpk
359         DO jj = nldj, nlej
360            DO ji = nldi, nlei
361               zva_tlin(ji,jj,jk) = znu(ji,jj,jk)
362            END DO
363         END DO
364      END DO
365      !--------------------------------------------------------------------
366      ! Call the tangent routine: dy = L dx
367      !--------------------------------------------------------------------
368
369      ua_tl(:,:,:) = zua_tlin(:,:,:)
370      va_tl(:,:,:) = zva_tlin(:,:,:)
371
372      CALL zdf_bfr_tan( nit000 )
373
374      zua_tlout(:,:,:) = ua_tl(:,:,:)   ;   zva_tlout(:,:,:) = va_tl(:,:,:)
375      zspgu_tlout(:,:) = spgu_tl(:,:)   ;   zspgv_tlout(:,:) = spgv_tl(:,:)
376
377      !--------------------------------------------------------------------
378      ! Initialize the adjoint variables: dy^* = W dy
379      !--------------------------------------------------------------------
380
381      DO jk = 1, jpk
382         DO jj = nldj, nlej
383            DO ji = nldi, nlei
384               zua_adin(ji,jj,jk) = zua_tlout(ji,jj,jk) &
385                  &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
386                  &               * umask(ji,jj,jk)
387               zva_adin(ji,jj,jk) = zva_tlout(ji,jj,jk) &
388                  &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
389                  &               * vmask(ji,jj,jk)
390            END DO
391         END DO
392      END DO
393      DO jj = nldj, nlej
394         DO ji = nldi, nlei
395            zspgu_adin (ji,jj) = zspgu_tlout (ji,jj) &
396               &              * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1) * umask(ji,jj,1)
397            zspgv_adin(ji,jj) = zspgv_tlout(ji,jj) &
398               &              * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1) * vmask(ji,jj,1)
399         END DO
400      END DO
401
402      !--------------------------------------------------------------------
403      ! Compute the scalar product: ( L dx )^T W dy
404      !--------------------------------------------------------------------
405
406      zsp1 =   DOT_PRODUCT( zua_tlout  , zua_adin   ) &
407         &   + DOT_PRODUCT( zspgu_tlout , zspgu_adin  ) &
408         &   + DOT_PRODUCT( zspgv_tlout , zspgv_adin  ) &
409         &   + DOT_PRODUCT( zva_tlout  , zva_adin   )
410
411
412      !--------------------------------------------------------------------
413      ! Call the adjoint routine: dx^* = L^T dy^*
414      !--------------------------------------------------------------------
415
416      ua_ad(:,:,:) = zua_adin(:,:,:)
417      va_ad(:,:,:) = zva_adin(:,:,:)
418
419      spgu_ad(:,:)   = zspgu_adin(:,:)
420      spgv_ad(:,:)   = zspgv_adin(:,:)
421
422      CALL zdf_bfr_adj( nitend )
423
424      zua_adout(:,:,:) = ua_ad(:,:,:)
425      zva_adout(:,:,:) = va_ad(:,:,:)
426
427      !--------------------------------------------------------------------
428      ! Compute the scalar product: dx^T L^T W dy
429      !--------------------------------------------------------------------
430
431      zsp2 =   DOT_PRODUCT( zua_tlin  , zua_adout   ) &
432         &   + DOT_PRODUCT( zva_tlin  , zva_adout   )
433
434      cl_name = 'zdf_bfr_adj   '
435      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
436
437      DEALLOCATE( &
438         & zua_tlin,  &
439         & zva_tlin,  &
440         & zua_tlout, &
441         & zva_tlout, &
442         & zua_adin,  &
443         & zva_adin,  &
444         & zua_adout, &
445         & zva_adout, &
446         & znu        &
447         & )
448   END SUBROUTINE zdf_bfr_adj_tst
449
450   !!======================================================================
451END MODULE zdfbfr_tam
Note: See TracBrowser for help on using the repository browser.