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_exp_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/dynzdf_exp_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: 11.0 KB
Line 
1MODULE dynzdf_exp_tam
2#ifdef key_tam
3   !!==============================================================================
4   !!                     ***  MODULE  dynzdf_exp_tam  ***
5   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
6   !!                 Tangent and Adjoint Module
7   !!==============================================================================
8   !! History of the direct module:
9   !!                !  90-10  (B. Blanke)  Original code
10   !!                !  97-05  (G. Madec)  vertical component of isopycnal
11   !!           8.5  !  02-08  (G. Madec)  F90: Free form and module
12   !! History of the TAM module:
13   !!           9.0  !  08-0!  (A. Vidard) Skeleton
14   !!           3.4  !  12-07  (P.-A. Bouttier) Phasing with 3.4
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   dyn_zdf_exp  : update the momentum trend with the vertical diffu-
19   !!                  sion using an explicit time-stepping scheme.
20   !!----------------------------------------------------------------------
21   !! * Modules used
22   USE par_oce
23   USE oce_tam
24   USE zdf_oce
25   USE dom_oce
26   USE phycst
27   USE in_out_manager
28   USE lib_mpp         ! MPP library
29   USE wrk_nemo        ! Memory Allocation
30   USE timing          ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   !! * Routine accessibility
36   PUBLIC dyn_zdf_exp_tan     ! called by dynzdf_tam.F90
37   PUBLIC dyn_zdf_exp_adj     ! called by dynzdf_tam.F90
38
39   !! * Substitutions
40#  include "domzgr_substitute.h90"
41#  include "vectopt_loop_substitute.h90"
42   !!----------------------------------------------------------------------
43
44CONTAINS
45
46   SUBROUTINE dyn_zdf_exp_tan( kt, p2dt )
47      !!----------------------------------------------------------------------
48      !!                  ***  ROUTINE dyn_zdf_exp_tan  ***
49      !!
50      !! ** Purpose of the direct routine:
51      !!      Compute the trend due to the vert. momentum diffusion
52      !!
53      !! ** Method of the direct routine:
54      !!      Explicit forward time stepping with a time splitting
55      !!      technique. The vertical diffusion of momentum is given by:
56      !!         diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ub) )
57      !!      Surface boundary conditions: wind stress input
58      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F90)
59      !!      Add this trend to the general trend ua :
60      !!         ua = ua + dz( avmu dz(u) )
61      !!
62      !! ** Action : - Update (ua,va) with the vertical diffusive trend
63      !!---------------------------------------------------------------------
64      !! * Arguments
65      INTEGER , INTENT( in ) ::   kt                           ! ocean time-step index
66      REAL(wp), INTENT( in ) ::   p2dt                         ! time-step
67
68      !! * Local declarations
69      INTEGER ::   ji, jj, jk, jl                              ! dummy loop indices
70      REAL(wp) ::   zrau0r, zlavmr, zuatl, zvatl                   ! temporary scalars
71      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwxtl, zwytl, zwztl, zwwtl ! temporary workspace arrays
72      !!----------------------------------------------------------------------
73      !
74      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_exp_adj')
75      !
76      CALL wrk_alloc( jpi,jpj,jpk, zwxtl, zwytl, zwztl, zwwtl )
77      !
78      IF( kt == nit000 .AND. lwp) THEN
79         IF(lwp) WRITE(numout,*)
80         IF(lwp) WRITE(numout,*) 'dyn_zdf_exp_tan : vertical momentum diffusion explicit operator'
81         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ '
82      ENDIF
83
84      ! Local constant initialization
85      ! -----------------------------
86      zrau0r = 1. / rau0                                   ! inverse of the reference density
87      zlavmr = 1. / REAL( nn_zdfexp )                      ! inverse of the number of sub time step
88      !                                                ! ===============
89                                                       !  Vertical slab
90      !                                                ! ===============
91      ! Surface boundary condition
92      DO jj = 2, jpjm1
93         DO ji = 2, jpim1
94            zwytl(ji,jj,1) = 0.0_wp
95            zwwtl(ji,jj,1) = 0.0_wp
96         END DO
97      END DO
98      ! Initialization of x, z and contingently trends array
99      DO jk = 1, jpk
100         DO jj = 2, jpjm1
101            DO ji = 2, jpim1
102               zwxtl(ji,jj,jk) = ub_tl(ji,jj,jk)
103               zwztl(ji,jj,jk) = vb_tl(ji,jj,jk)
104            END DO
105         END DO
106      END DO
107      ! Time splitting loop
108      DO jl = 1, nn_zdfexp
109         !
110         ! First vertical derivative
111         DO jk = 2, jpk
112            DO jj = 2, jpjm1
113               DO ji = 2, jpim1
114                  zwytl(ji,jj,jk) = avmu(ji,jj,jk) * ( zwxtl(ji,jj,jk-1) - zwxtl(ji,jj,jk) ) / fse3uw(ji,jj,jk)
115                  zwwtl(ji,jj,jk) = avmv(ji,jj,jk) * ( zwztl(ji,jj,jk-1) - zwztl(ji,jj,jk) ) / fse3vw(ji,jj,jk)
116               END DO
117            END DO
118         END DO
119         ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp
120         DO jk = 1, jpkm1
121            DO jj = 2, jpjm1
122               DO ji = 2, jpim1
123                  zuatl = zlavmr*( zwytl(ji,jj,jk) - zwytl(ji,jj,jk+1) ) / fse3u(ji,jj,jk)
124                  zvatl = zlavmr*( zwwtl(ji,jj,jk) - zwwtl(ji,jj,jk+1) ) / fse3v(ji,jj,jk)
125                  ua_tl(ji,jj,jk) = ua_tl(ji,jj,jk) + zuatl
126                  va_tl(ji,jj,jk) = va_tl(ji,jj,jk) + zvatl
127                  zwxtl(ji,jj,jk) = zwxtl(ji,jj,jk) + p2dt*zuatl*umask(ji,jj,jk)
128                  zwztl(ji,jj,jk) = zwztl(ji,jj,jk) + p2dt*zvatl*vmask(ji,jj,jk)
129               END DO
130            END DO
131         END DO
132         !
133      END DO
134      !                                                ! ===============
135      !                                                !   End of slab
136      !                                                ! ===============
137      !
138      CALL wrk_dealloc( jpi,jpj,jpk, zwxtl, zwytl, zwztl, zwwtl )
139      !
140      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_exp_tan')
141      !
142   END SUBROUTINE dyn_zdf_exp_tan
143   SUBROUTINE dyn_zdf_exp_adj( kt, p2dt )
144      !!----------------------------------------------------------------------
145      !!                  ***  ROUTINE dyn_zdf_exp_adj  ***
146      !!
147      !! ** Purpose of the direct routine:
148      !!      Compute the trend due to the vert. momentum diffusion
149      !!
150      !! ** Method of the direct routine:
151      !!      Explicit forward time stepping with a time splitting
152      !!      technique. The vertical diffusion of momentum is given by:
153      !!         diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ub) )
154      !!      Surface boundary conditions: wind stress input
155      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F90)
156      !!      Add this trend to the general trend ua :
157      !!         ua = ua + dz( avmu dz(u) )
158      !!
159      !! ** Action : - Update (ua,va) with the vertical diffusive trend
160      !!---------------------------------------------------------------------
161      !! * Arguments
162      INTEGER , INTENT( in ) ::   kt                           ! ocean time-step index
163      REAL(wp), INTENT( in ) ::   p2dt                         ! time-step
164
165      !! * Local declarations
166      INTEGER ::   ji, jj, jk, jl                                  ! dummy loop indices
167      REAL(wp) ::   zrau0r, zlavmr, zuaad, zvaad                   ! temporary scalars
168      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zwxad, zwyad, zwzad, zwwad ! temporary workspace arrays
169      !!----------------------------------------------------------------------
170      !
171      IF( nn_timing == 1 )  CALL timing_start('dyn_zdf_exp_adj')
172      !
173      CALL wrk_alloc( jpi,jpj,jpk, zwxad, zwyad, zwzad, zwwad )
174      !
175      IF( kt == nitend ) THEN
176         IF(lwp) WRITE(numout,*)
177         IF(lwp) WRITE(numout,*) 'dyn_zdf_exp_adj : vertical momentum diffusion explicit operator'
178         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ '
179      ENDIF
180      ! Local constant initialization
181      ! -----------------------------
182      zrau0r = 1. / rau0                                   ! inverse of the reference density
183      zlavmr = 1. / float( nn_zdfexp )                      ! inverse of the number of sub time step
184      zwxad(:,:,:) = 0.0_wp ; zwyad(:,:,:) = 0.0_wp ; zwzad(:,:,:) = 0.0_wp ; zwwad(:,:,:) = 0.0_wp
185      !                                                ! ===============
186      !                                                !  Vertical slab
187      !                                                ! ===============
188      ! Time splitting loop
189      DO jl = 1, nn_zdfexp
190         !
191         ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp
192         DO jk = 1, jpkm1
193            DO jj = 2, jpjm1
194               DO ji = 2, jpim1
195                  zuaad = p2dt * zwxad(ji,jj,jk) * umask(ji,jj,jk)
196                  zvaad = p2dt * zwzad(ji,jj,jk) * vmask(ji,jj,jk)
197                  zuaad = zuaad + ua_ad(ji,jj,jk)
198                  zvaad = zvaad + va_ad(ji,jj,jk)
199                  zwyad(ji,jj,jk  ) = zwyad(ji,jj,jk  ) + zlavmr * zuaad  / fse3u(ji,jj,jk)
200                  zwyad(ji,jj,jk+1) = zwyad(ji,jj,jk+1) - zlavmr * zuaad  / fse3u(ji,jj,jk)
201                  zwwad(ji,jj,jk  ) = zwwad(ji,jj,jk  ) + zlavmr * zvaad  / fse3v(ji,jj,jk)
202                  zwwad(ji,jj,jk+1) = zwwad(ji,jj,jk+1) - zlavmr * zvaad  / fse3v(ji,jj,jk)
203               END DO
204            END DO
205         END DO
206         ! First vertical derivative
207         DO jk = 2, jpk
208            DO jj = 2, jpjm1
209               DO ji = 2, jpim1
210                  zwxad(ji,jj,jk-1) = zwxad(ji,jj,jk-1) + avmu(ji,jj,jk) * zwyad(ji,jj,jk) / fse3uw(ji,jj,jk)
211                  zwxad(ji,jj,jk  ) = zwxad(ji,jj,jk  ) - avmu(ji,jj,jk) * zwyad(ji,jj,jk) / fse3uw(ji,jj,jk)
212                  zwyad(ji,jj,jk  ) = 0.0_wp
213
214                  zwzad(ji,jj,jk-1) = zwzad(ji,jj,jk-1) + avmv(ji,jj,jk) * zwwad(ji,jj,jk) / fse3vw(ji,jj,jk)
215                  zwzad(ji,jj,jk  ) = zwzad(ji,jj,jk  ) - avmv(ji,jj,jk) * zwwad(ji,jj,jk) / fse3vw(ji,jj,jk)
216                  zwwad(ji,jj,jk  ) = 0.0_wp
217               END DO
218            END DO
219         END DO
220         !
221      END DO
222      !
223      ! Initialization of x, z and contingently trends array
224      DO jk = 1, jpk
225         DO jj = 2, jpjm1
226            DO ji = 2, jpim1
227               ub_ad(ji,jj,jk) = ub_ad(ji,jj,jk) + zwxad(ji,jj,jk)
228               vb_ad(ji,jj,jk) = vb_ad(ji,jj,jk) + zwzad(ji,jj,jk)
229               zwxad(ji,jj,jk) = 0.0_wp
230               zwzad(ji,jj,jk) = 0.0_wp
231            END DO
232         END DO
233      END DO
234      !                                                ! ===============
235      !                                                !   End of slab
236      !                                                ! ===============
237      !
238      CALL wrk_dealloc( jpi,jpj,jpk, zwxad, zwyad, zwzad, zwwad )
239      !
240      IF( nn_timing == 1 )  CALL timing_stop('dyn_zdf_exp_adj')
241      !
242   END SUBROUTINE dyn_zdf_exp_adj
243#endif
244   !!==============================================================================
245END MODULE dynzdf_exp_tam
Note: See TracBrowser for help on using the repository browser.