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 @ 247

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

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