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.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynzdf_exp.F90 @ 216

Last change on this file since 216 was 216, checked in by opalod, 19 years ago

CT : UPDATE151 : New trends organization
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.9 KB
Line 
1MODULE dynzdf_exp
2   !!==============================================================================
3   !!                     ***  MODULE  dynzdf_exp  ***
4   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   dyn_zdf_exp  : update the momentum trend with the vertical diffu-
9   !!                  sion using an explicit time-stepping scheme.
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE oce             ! ocean dynamics and tracers
13   USE dom_oce         ! ocean space and time domain
14   USE phycst          ! physical constants
15   USE zdf_oce         ! ocean vertical physics
16   USE in_out_manager  ! I/O manager
17   USE taumod          ! surface ocean stress
18   USE trdmod          ! ocean dynamics trends
19   USE trdmod_oce      ! ocean variables trends
20
21   IMPLICIT NONE
22   PRIVATE
23
24   !! * Routine accessibility
25   PUBLIC dyn_zdf_exp    ! called by step.F90
26
27   !! * Substitutions
28#  include "domzgr_substitute.h90"
29#  include "vectopt_loop_substitute.h90"
30   !!----------------------------------------------------------------------
31   !!   OPA 9.0 , LODYC-IPSL  (2003)
32   !!----------------------------------------------------------------------
33
34CONTAINS
35
36   SUBROUTINE dyn_zdf_exp( kt )
37      !!----------------------------------------------------------------------
38      !!                  ***  ROUTINE dyn_zdf_exp  ***
39      !!                   
40      !! ** Purpose :   Compute the trend due to the vert. momentum diffusion
41      !!
42      !! ** Method  :   Explicit forward time stepping with a time splitting
43      !!      technique. The vertical diffusion of momentum is given by:
44      !!         diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ub) )
45      !!      Surface boundary conditions: wind stress input
46      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F90)
47      !!      Add this trend to the general trend ua :
48      !!         ua = ua + dz( avmu dz(u) )
49      !!
50      !! ** Action : - Update (ua,va) with the vertical diffusive trend
51      !!             - Save the trends in (ztdua,ztdva) ('key_trddyn')
52      !!
53      !! History :
54      !!        !  90-10  (B. Blanke)  Original code
55      !!        !  97-05  (G. Madec)  vertical component of isopycnal
56      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
57      !!   9.0  !  04-08  (C. Talandier)  New trends organization
58      !!---------------------------------------------------------------------
59      !! * Modules used     
60      USE oce, ONLY :    ztdua => ta,    & ! use ta as 3D workspace   
61                         ztdva => sa       ! use sa as 3D workspace   
62      !! * Arguments
63      INTEGER, INTENT( in ) ::   kt        ! ocean time-step index
64
65      !! * Local declarations
66      INTEGER ::   &
67         ji, jj, jk, jl,                 & ! dummy loop indices
68         ikbu, ikbum1 , ikbv, ikbvm1       ! temporary integers
69      REAL(wp) ::   &
70         zrau0r, zlavmr, z2dt, zua, zva    ! temporary scalars
71      REAL(wp), DIMENSION(jpi,jpk) ::    &
72         zwx, zwy, zwz, zww                ! temporary workspace arrays
73      REAL(wp), DIMENSION(jpi,jpj) ::    &
74         ztsx, ztsy, ztbx, ztby            ! temporary workspace arrays
75      !!----------------------------------------------------------------------
76
77      IF( kt == nit000 ) THEN
78         IF(lwp) WRITE(numout,*)
79         IF(lwp) WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion explicit operator'
80         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
81      ENDIF
82
83      ! Local constant initialization
84      ! -----------------------------
85      zrau0r = 1. / rau0                                   ! inverse of the reference density
86      zlavmr = 1. / float( n_zdfexp )                      ! inverse of the number of sub time step
87      z2dt = 2. * rdt                                      ! Leap-frog environnement
88      ztsx(:,:) = 0.e0
89      ztsy(:,:) = 0.e0 
90      ztbx(:,:) = 0.e0
91      ztby(:,:) = 0.e0
92
93      ! Save ua and va trends
94      IF( l_trddyn )   THEN
95         ztdua(:,:,:) = ua(:,:,:) 
96         ztdva(:,:,:) = va(:,:,:) 
97      ENDIF
98
99      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt    ! Euler time stepping when starting from rest
100
101      !                                                ! ===============
102      DO jj = 2, jpjm1                                 !  Vertical slab
103         !                                             ! ===============
104
105         ! Surface boundary condition
106         DO ji = 2, jpim1
107            zwy(ji,1) = taux(ji,jj) * zrau0r
108            zww(ji,1) = tauy(ji,jj) * zrau0r
109         END DO 
110
111         ! Initialization of x, z and contingently trends array
112         DO jk = 1, jpk
113            DO ji = 2, jpim1
114               zwx(ji,jk) = ub(ji,jj,jk)
115               zwz(ji,jk) = vb(ji,jj,jk)
116            END DO 
117         END DO 
118
119         ! Time splitting loop
120         DO jl = 1, n_zdfexp
121
122            ! First vertical derivative
123            DO jk = 2, jpk
124               DO ji = 2, jpim1
125                  zwy(ji,jk) = avmu(ji,jj,jk) * ( zwx(ji,jk-1) - zwx(ji,jk) ) / fse3uw(ji,jj,jk) 
126                  zww(ji,jk) = avmv(ji,jj,jk) * ( zwz(ji,jk-1) - zwz(ji,jk) ) / fse3vw(ji,jj,jk)
127               END DO 
128            END DO 
129
130            ! Second vertical derivative and trend estimation at kt+l*rdt/n_zdfexp
131            DO jk = 1, jpkm1
132               DO ji = 2, jpim1
133                  zua = zlavmr*( zwy(ji,jk) - zwy(ji,jk+1) ) / fse3u(ji,jj,jk)
134                  zva = zlavmr*( zww(ji,jk) - zww(ji,jk+1) ) / fse3v(ji,jj,jk)
135                  ua(ji,jj,jk) = ua(ji,jj,jk) + zua
136                  va(ji,jj,jk) = va(ji,jj,jk) + zva
137
138                  zwx(ji,jk) = zwx(ji,jk) + z2dt*zua*umask(ji,jj,jk)
139                  zwz(ji,jk) = zwz(ji,jk) + z2dt*zva*vmask(ji,jj,jk)
140               END DO 
141            END DO 
142
143         END DO 
144
145         !                                             ! ===============
146      END DO                                           !   End of slab
147      !                                                ! ===============
148
149      ! save the vertical diffusion trends for diagnostic
150      ! momentum trends
151      IF( l_trddyn )  THEN 
152         ! save the total vertical momentum diffusive trend
153         ztdua(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
154         ztdva(:,:,:) = va(:,:,:) - ztdva(:,:,:)
155 
156         ! subtract and save surface and momentum fluxes
157         !                                                ! ===============
158         DO jj = 2, jpjm1                                 !  Horizontal slab
159            !                                             ! ===============
160            DO ji = 2, jpim1
161               ! save the surface momentum fluxes
162               ztsx(ji,jj) = zwy(ji,1) / fse3u(ji,jj,1)
163               ztsy(ji,jj) = zww(ji,1) / fse3v(ji,jj,1)
164               ! save bottom friction momentum fluxes
165               ikbu   = MIN( mbathy(ji+1,jj), mbathy(ji,jj) )
166               ikbum1 = MAX( ikbu-1, 1 )
167               ikbv   = MIN( mbathy(ji,jj+1), mbathy(ji,jj) )
168               ikbvm1 = MAX( ikbv-1, 1 )
169               ztbx(ji,jj) = avmu(ji,jj,ikbu) * zwx(ji,ikbum1)   &
170                               / ( fse3u(ji,jj,ikbum1) * fse3uw(ji,jj,ikbu) )
171               ztby(ji,jj) = avmv(ji,jj,ikbv) * zwz(ji,ikbvm1)   &
172                               / ( fse3v(ji,jj,ikbvm1) * fse3vw(ji,jj,ikbv) )
173               ! subtract surface forcing and bottom friction trend from vertical
174               ! diffusive momentum trend
175               ztdua(ji,jj,1     ) = ztdua(ji,jj,1     ) - ztsx(ji,jj)
176               ztdua(ji,jj,ikbum1) = ztdua(ji,jj,ikbum1) - ztbx(ji,jj)
177               ztdva(ji,jj,1     ) = ztdva(ji,jj,1     ) - ztsy(ji,jj)
178               ztdva(ji,jj,ikbvm1) = ztdva(ji,jj,ikbvm1) - ztby(ji,jj)
179            END DO
180            !                                             ! ===============
181         END DO                                           !   End of slab
182         !                                                ! ===============
183
184         CALL trd_mod(ztdua, ztdva, jpdtdzdf, 'DYN', kt)
185         ztdua(:,:,:) = 0.e0
186         ztdva(:,:,:) = 0.e0
187         ztdua(:,:,1) = ztsx(:,:)
188         ztdva(:,:,1) = ztsy(:,:)
189         CALL trd_mod(ztdua , ztdva , jpdtdswf, 'DYN', kt)
190         ztdua(:,:,:) = 0.e0
191         ztdva(:,:,:) = 0.e0
192         ztdua(:,:,1) = ztbx(:,:)
193         ztdva(:,:,1) = ztby(:,:)
194         CALL trd_mod(ztdua , ztdva , jpdtdbfr, 'DYN', kt)
195      ENDIF
196
197      IF(l_ctl) THEN         ! print sum trends (used for debugging)
198         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
199         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
200         WRITE(numout,*) ' zdf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
201         u_ctl = zua   ;   v_ctl = zva
202      ENDIF
203
204   END SUBROUTINE dyn_zdf_exp
205
206   !!==============================================================================
207END MODULE dynzdf_exp
Note: See TracBrowser for help on using the repository browser.