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.
trazdf_imp.F90 in tags/nemo_v1_02/NEMO/OPA_SRC/TRA – NEMO

source: tags/nemo_v1_02/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 3452

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

nemo_v1_update_002 : CT : Integration of the KPP turbulent closure scheme
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.1 KB
Line 
1MODULE trazdf_imp
2   !!==============================================================================
3   !!                    ***  MODULE  trazdf_imp  ***
4   !! Ocean active tracers:  vertical component of the tracer mixing trend
5   !!==============================================================================
6
7   !!----------------------------------------------------------------------
8   !!   tra_zdf_imp  : update the tracer trend with the vertical diffusion
9   !!                  using an implicit time-stepping.
10   !!----------------------------------------------------------------------
11   !! * Modules used
12   USE oce             ! ocean dynamics and active tracers
13   USE dom_oce         ! ocean space and time domain
14   USE zdf_oce         ! ocean vertical physics
15   USE ldftra_oce      ! ocean active tracers: lateral physics
16   USE zdfddm          ! ocean vertical physics: double diffusion
17   USE zdfkpp          ! KPP parameterisation
18   USE trdmod          ! ocean active tracers trends
19   USE trdmod_oce      ! ocean variables trends
20   USE in_out_manager  ! I/O manager
21
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Routine accessibility
27   PUBLIC tra_zdf_imp          ! routine called by step.F90
28
29   !! * Module variable
30   REAL(wp), DIMENSION(jpk) ::   &
31      r2dt                     ! vertical profile of 2 x tracer time-step
32
33   !! * Substitutions
34#  include "domzgr_substitute.h90"
35#  include "zdfddm_substitute.h90"
36   !!----------------------------------------------------------------------
37   !!   OPA 9.0 , LOCEAN-IPSL (2005)
38   !! $Header$
39   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
40   !!----------------------------------------------------------------------
41
42CONTAINS
43
44   SUBROUTINE tra_zdf_imp( kt )
45      !!----------------------------------------------------------------------
46      !!                  ***  ROUTINE tra_zdf_imp  ***
47      !!
48      !! ** Purpose :   Compute the trend due to the vertical tracer mixing
49      !!      using an implicit time stepping and add it to the general trend
50      !!      of the tracer equations.
51      !!
52      !! ** Method  :   The vertical diffusion of tracers (t & s) is given by:
53      !!          difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(ta) )
54      !!      It is thus evaluated using a backward time scheme
55      !!      Surface and bottom boundary conditions: no diffusive flux on
56      !!      both tracers (bottom, applied through the masked field avt).
57      !!      Add this trend to the general trend ta,sa :
58      !!          ta = ta + dz( avt dz(t) )
59      !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T)
60      !!
61      !! ** Action  : - Update (ta,sa) with the before vertical diffusion trend
62      !!              - save the trends in (ttrd,strd) ('key_trdtra')
63      !!
64      !! History :
65      !!   6.0  !  90-10  (B. Blanke)  Original code
66      !!   7.0  !  91-11  (G. Madec)
67      !!        !  92-06  (M. Imbard)  correction on tracer trend loops
68      !!        !  96-01  (G. Madec)  statement function for e3
69      !!        !  97-05  (G. Madec)  vertical component of isopycnal
70      !!        !  97-07  (G. Madec)  geopotential diffusion in s-coord
71      !!        !  00-08  (G. Madec)  double diffusive mixing
72      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
73      !!   9.0  !  04-08  (C. Talandier) New trends organization
74      !!   9.0  !  05-01  (C. Ethe )  non-local flux in KPP vertical mixing scheme
75      !!---------------------------------------------------------------------
76      !! * Modules used     
77      USE oce, ONLY :    ztdta => ua,       & ! use ua as 3D workspace   
78                         ztdsa => va          ! use va as 3D workspace   
79
80      !! * Arguments
81      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
82
83      !! * Local declarations
84      INTEGER ::   ji, jj, jk                 ! dummy loop indices
85      INTEGER ::   ikst, ikenm2, ikstp1
86      REAL(wp)  ::   zta, zsa                 ! temporary scalars
87      REAL(wp), DIMENSION(jpi,jpk) ::   &
88         zwd, zws, zwi,          &  ! ???
89         zwx, zwy, zwz, zwt         ! ???
90      !!---------------------------------------------------------------------
91
92
93      ! 0. Local constant initialization
94      ! --------------------------------
95
96      ! time step = 2 rdttra ex
97      IF( neuler == 0 .AND. kt == nit000 ) THEN
98         r2dt(:) =  rdttra(:)              ! restarting with Euler time stepping
99      ELSEIF( kt <= nit000 + 1) THEN
100         r2dt(:) = 2. * rdttra(:)          ! leapfrog
101      ENDIF
102
103      ! Save ta and sa trends
104      IF( l_trdtra )   THEN
105         ztdta(:,:,:) = ta(:,:,:) 
106         ztdsa(:,:,:) = sa(:,:,:) 
107      ENDIF
108
109      !                                                ! ===============
110      DO jj = 2, jpjm1                                 !  Vertical slab
111         !                                             ! ===============
112         ! 0. Matrix construction
113         ! ----------------------
114
115         ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked)
116         DO jk = 1, jpkm1
117            DO ji = 2, jpim1
118               zwi(ji,jk) = - r2dt(jk) * avt(ji,jj,jk  )   &
119                                       / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
120               zws(ji,jk) = - r2dt(jk) * avt(ji,jj,jk+1)   &
121                                       / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
122               zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk)
123            END DO
124         END DO
125
126         ! Surface boudary conditions
127         DO ji = 2, jpim1
128            zwi(ji,1) = 0.e0
129            zwd(ji,1) = 1. - zws(ji,1)
130         END DO
131
132         ! 1. Vertical diffusion on temperature
133         ! -------------------------===========
134
135         ! Second member construction
136#if defined key_zdfkpp
137         ! add non-local temperature flux ( in convective case only)
138         DO jk = 1, jpkm1
139            DO ji = 2, jpim1 
140               zwy(ji,jk) = tb(ji,jj,jk) + r2dt(jk) * ta(ji,jj,jk)  &
141                  &  - r2dt(jk) * ( ghats(ji,jj,jk) * avt(ji,jj,jk) - ghats(ji,jj,jk+1) * avt(ji,jj,jk+1) ) &
142                  &               * wt0(ji,jj) / fse3t(ji,jj,jk) 
143            END DO
144         END DO
145#else
146         DO jk = 1, jpkm1
147            DO ji = 2, jpim1             
148               zwy(ji,jk) = tb(ji,jj,jk) + r2dt(jk) * ta(ji,jj,jk)
149            END DO
150         END DO
151#endif
152
153         ! Matrix inversion from the first level
154         ikst = 1
155
156#   include "zdf.matrixsolver.h90"
157
158         ! Save the masked temperature after in ta
159         ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done
160         !                  it will not be done in tranxt)
161         DO jk = 1, jpkm1
162            DO ji = 2, jpim1
163               ta(ji,jj,jk) = zwx(ji,jk) * tmask(ji,jj,jk)
164            END DO
165         END DO
166
167
168         ! 2. Vertical diffusion on salinity
169         ! -------------------------========
170
171#if defined key_zdfddm
172         ! Rebuild the Matrix as avt /= avs
173
174         ! Diagonal, inferior, superior
175         ! (including the bottom boundary condition via avs masked)
176         DO jk = 1, jpkm1
177            DO ji = 2, jpim1
178               zwi(ji,jk) = - r2dt(jk) * fsavs(ji,jj,jk  )   &
179                  /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk  ) )
180               zws(ji,jk) = - r2dt(jk) * fsavs(ji,jj,jk+1)   &
181                  /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) )
182               zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk)
183            END DO
184         END DO
185
186         ! Surface boudary conditions
187         DO ji = 2, jpim1
188            zwi(ji,1) = 0.e0
189            zwd(ji,1) = 1. - zws(ji,1)
190         END DO
191#endif
192         ! Second member construction
193#if defined key_zdfkpp
194         ! add non-local salinity flux ( in convective case only)
195         DO jk = 1, jpkm1
196            DO ji = 2, jpim1 
197               zwy(ji,jk) = sb(ji,jj,jk) + r2dt(jk) * sa(ji,jj,jk)  &
198                  &  - r2dt(jk) * ( ghats(ji,jj,jk) * fsavs(ji,jj,jk) - ghats(ji,jj,jk+1) * fsavs(ji,jj,jk+1) ) &
199                  &               * ws0(ji,jj) / fse3t(ji,jj,jk) 
200            END DO
201         END DO
202#else
203         DO jk = 1, jpkm1
204            DO ji = 2, jpim1             
205               zwy(ji,jk) = sb(ji,jj,jk) + r2dt(jk) * sa(ji,jj,jk)
206            END DO
207         END DO
208#endif
209 
210         ! Matrix inversion from the first level
211         ikst = 1
212
213#   include "zdf.matrixsolver.h90"
214
215
216         ! Save the masked salinity after in sa
217         ! (c a u t i o n:  salinity not its trend, Leap-frog scheme done
218         !                  it will not be done in tranxt)
219         DO jk = 1, jpkm1
220            DO ji = 2, jpim1
221               sa(ji,jj,jk) = zwx(ji,jk)  * tmask(ji,jj,jk)
222            END DO
223         END DO
224
225         !                                             ! ===============
226      END DO                                           !   End of slab
227      !                                                ! ===============
228
229      ! save the trends for diagnostic
230      ! Compute and save the vertical diffusive temperature & salinity trends
231      IF( l_trdtra )   THEN
232         ! compute the vertical diffusive trends in substracting the previous
233         ! trends ztdta()/ztdsa() to the new one computed (dT/dt or dS/dt)
234         ! with the new temperature/salinity ta/sa
235         DO jk = 1, jpkm1
236            ztdta(:,:,jk) = ( ( ta(:,:,jk) - tb(:,:,jk) ) / r2dt(jk) )   & ! new trend
237                &           - ztdta(:,:,jk)                                ! old trend
238            ztdsa(:,:,jk) = ( ( sa(:,:,jk) - sb(:,:,jk) ) / r2dt(jk) )   & ! new trend
239                &           - ztdsa(:,:,jk)                                ! old trend
240         END DO
241
242         CALL trd_mod(ztdta, ztdsa, jpttdzdf, 'TRA', kt)
243      ENDIF
244
245      IF(l_ctl) THEN         ! print mean trends (used for debugging)
246         zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
247         zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
248         WRITE(numout,*) ' zdf  - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
249         t_ctl = zta   ;   s_ctl = zsa
250      ENDIF
251
252   END SUBROUTINE tra_zdf_imp
253
254   !!==============================================================================
255END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.