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 trunk/NEMO/OPA_SRC/TRA – NEMO

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