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 branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90 @ 2598

Last change on this file since 2598 was 2590, checked in by trackstand2, 13 years ago

Merge branch 'dynamic_memory' into master-svn-dyn
  • Property svn:keywords set to Id
File size: 15.3 KB
RevLine 
[3]1MODULE trazdf_imp
[1438]2   !!======================================================================
[457]3   !!                 ***  MODULE  trazdf_imp  ***
[2528]4   !! Ocean  tracers:  vertical component of the tracer mixing trend
[1438]5   !!======================================================================
6   !! History :  OPA  !  1990-10  (B. Blanke)  Original code
7   !!            7.0  !  1991-11  (G. Madec)
8   !!                 !  1992-06  (M. Imbard) correction on tracer trend loops
9   !!                 !  1996-01  (G. Madec) statement function for e3
10   !!                 !  1997-05  (G. Madec) vertical component of isopycnal
11   !!                 !  1997-07  (G. Madec) geopotential diffusion in s-coord
12   !!                 !  2000-08  (G. Madec) double diffusive mixing
13   !!   NEMO     1.0  !  2002-08  (G. Madec) F90: Free form and module
14   !!            2.0  !  2006-11  (G. Madec) New step reorganisation
15   !!            3.2  !  2009-03  (G. Madec)  heat and salt content trends
[2528]16   !!            3.3  !  2010-06  (C. Ethe, G. Madec) Merge TRA-TRC
[3]17   !!----------------------------------------------------------------------
[1438]18 
19   !!----------------------------------------------------------------------
[457]20   !!   tra_zdf_imp : Update the tracer trend with the diagonal vertical 
21   !!                 part of the mixing tensor.
[3]22   !!----------------------------------------------------------------------
[457]23   USE oce             ! ocean dynamics and tracers variables
24   USE dom_oce         ! ocean space and time domain variables
25   USE zdf_oce         ! ocean vertical physics variables
[2528]26   USE trc_oce         ! share passive tracers/ocean variables
27   USE domvvl          ! variable volume
[216]28   USE ldftra_oce      ! ocean active tracers: lateral physics
[2528]29   USE ldftra          ! lateral mixing type
[457]30   USE ldfslp          ! lateral physics: slope of diffusion
31   USE zdfddm          ! ocean vertical physics: double diffusion
[2528]32   USE traldf_iso_grif ! active tracers: Griffies operator
[3]33   USE in_out_manager  ! I/O manager
[457]34   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[3]35
36   IMPLICIT NONE
37   PRIVATE
38
[1438]39   PUBLIC   tra_zdf_imp   !  routine called by step.F90
[3]40
41   !! * Substitutions
42#  include "domzgr_substitute.h90"
[457]43#  include "ldftra_substitute.h90"
[3]44#  include "zdfddm_substitute.h90"
[457]45#  include "vectopt_loop_substitute.h90"
[3]46   !!----------------------------------------------------------------------
[2528]47   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1156]48   !! $Id$
[2528]49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[457]50   !!----------------------------------------------------------------------
[3]51CONTAINS
[2528]52 
53   SUBROUTINE tra_zdf_imp( kt, cdtype, p2dt, ptb, pta, kjpt ) 
[3]54      !!----------------------------------------------------------------------
55      !!                  ***  ROUTINE tra_zdf_imp  ***
56      !!
[457]57      !! ** Purpose :   Compute the trend due to the vertical tracer diffusion
58      !!     including the vertical component of lateral mixing (only for 2nd
59      !!     order operator, for fourth order it is already computed and add
60      !!     to the general trend in traldf.F) and add it to the general trend
61      !!     of the tracer equations.
[3]62      !!
[457]63      !! ** Method  :   The vertical component of the lateral diffusive trends
[503]64      !!      is provided by a 2nd order operator rotated along neutral or geo-
[457]65      !!      potential surfaces to which an eddy induced advection can be
66      !!      added. It is computed using before fields (forward in time) and
67      !!      isopycnal or geopotential slopes computed in routine ldfslp.
68      !!
69      !!      Second part: vertical trend associated with the vertical physics
70      !!      ===========  (including the vertical flux proportional to dk[t]
71      !!                  associated with the lateral mixing, through the
72      !!                  update of avt)
[2528]73      !!      The vertical diffusion of the tracer t  is given by:
[457]74      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
75      !!      It is computed using a backward time scheme (t=ta).
[3]76      !!      Surface and bottom boundary conditions: no diffusive flux on
77      !!      both tracers (bottom, applied through the masked field avt).
78      !!      Add this trend to the general trend ta,sa :
[457]79      !!         ta = ta + dz( avt dz(t) )
[2528]80      !!         if lk_zdfddm=T, use avs for salinity or for passive tracers
81      !!         (sa = sa + dz( avs dz(t) )
[3]82      !!
[457]83      !!      Third part: recover avt resulting from the vertical physics
84      !!      ==========  alone, for further diagnostics (for example to
85      !!                  compute the turbocline depth in zdfmxl.F90).
[2528]86      !!         if lk_zdfddm=T, use avt = zavt
[457]87      !!         (avs = zavs if lk_zdfddm=T )
[3]88      !!
[2528]89      !! ** Action  : - Update (ta) with before vertical diffusion trend
[3]90      !!---------------------------------------------------------------------
[1438]91      USE oce    , ONLY :   zwd   => ua   ! ua used as workspace
92      USE oce    , ONLY :   zws   => va   ! va  -          -
[2590]93      USE wrk_nemo, ONLY: wrk_use, wrk_release
94      USE wrk_nemo, ONLY: zwi => wrk_3d_1, zwt => wrk_3d_2  ! workspace arrays
[2528]95      !!
96      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
97      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
98      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
99      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt     ! vertical profile of tracer time-step
100      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb      ! before and now tracer fields
101      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta      ! tracer trend
[1438]102      !!
[2528]103      INTEGER  ::  ji, jj, jk, jn        ! dummy loop indices
104      REAL(wp) ::  zavi, zrhs, znvvl     ! local scalars
105      REAL(wp) ::  ze3tb, ze3tn, ze3ta   ! variable vertical scale factors
[3]106      !!---------------------------------------------------------------------
107
[2590]108      IF(.NOT. wrk_use(3, 1,2))THEN
109         CALL ctl_stop('tra_zdf_imp : requested workspace arrays unavailable.')
110         RETURN
111      END IF
112
[2528]113      IF( kt == nit000 )  THEN
[457]114         IF(lwp)WRITE(numout,*)
[2528]115         IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype
[457]116         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
[2528]117         zavi = 0._wp      ! avoid warning at compilation phase when lk_ldfslp=F
[457]118      ENDIF
[2528]119      !
[457]120      ! I. Local initialization
121      ! -----------------------
[2528]122      zwd(1,:, : ) = 0._wp     ;     zwd(jpi,:,:) = 0._wp
123      zws(1,:, : ) = 0._wp     ;     zws(jpi,:,:) = 0._wp
124      zwi(1,:, : ) = 0._wp     ;     zwi(jpi,:,:) = 0._wp
125      zwt(1,:, : ) = 0._wp     ;     zwt(jpi,:,:) = 0._wp
126      zwt(:,:,jpk) = 0._wp     ;     zwt( : ,:,1) = 0._wp
[3]127
[592]128      ! I.1 Variable volume : to take into account vertical variable vertical scale factors
129      ! -------------------
[2528]130      IF( lk_vvl ) THEN   ;    znvvl = 1._wp
131      ELSE                ;    znvvl = 0._wp
[592]132      ENDIF
133
[457]134      ! II. Vertical trend associated with the vertical physics
135      ! =======================================================
136      !     (including the vertical flux proportional to dk[t] associated
137      !      with the lateral mixing, through the avt update)
138      !     dk[ avt dk[ (t,s) ] ] diffusive trends
[3]139
[2528]140      !
[457]141      ! II.0 Matrix construction
142      ! ------------------------
[2528]143      DO jn = 1, kjpt
144         !
145         !  Matrix construction
146         ! ------------------------
147         IF( cdtype == 'TRA' .AND. jn == jp_tem )  THEN 
[457]148#if defined key_ldfslp
[2528]149            IF( ln_traldf_grif ) THEN
150               DO jk = 2, jpkm1
151                  DO jj = 2, jpjm1
152                     DO ji = fs_2, fs_jpim1   ! vector opt.
153                        ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
154                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing
155                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
156                     END DO
157                  END DO
158               END DO
159            ! update and save of avt (and avs if double diffusive mixing)
160            ELSE IF( l_traldf_rot ) THEN
161               DO jk = 2, jpkm1
162                  DO jj = 2, jpjm1
163                     DO ji = fs_2, fs_jpim1   ! vector opt.
164                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
165                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
166                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
167                        zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
168                     END DO
169                  END DO
170               END DO
171            ELSE                         ! no rotation but key_ldfslp defined
172               zwt(:,:,:) = avt(:,:,:)
173            ENDIF
174#else
175            ! No isopycnal diffusion
176            zwt(:,:,:) = avt(:,:,:)           
177#endif
178            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
179            DO jk = 1, jpkm1
180               DO jj = 2, jpjm1
181                  DO ji = fs_2, fs_jpim1   ! vector opt.
182                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
183                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
184                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
185                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
186                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
187                 END DO
188               END DO
189            END DO
190            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
[902]191            DO jj = 2, jpjm1
[2528]192               DO ji = fs_2, fs_jpim1
193                  zwt(ji,jj,1) = zwd(ji,jj,1)
[902]194               END DO
[3]195            END DO
[2528]196            DO jk = 2, jpkm1
197               DO jj = 2, jpjm1
198                  DO ji = fs_2, fs_jpim1
199                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
200                  END DO
201               END DO
202            END DO
203            !
204         ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
205#if defined key_ldfslp
206            IF( ln_traldf_grif ) THEN
207               DO jk = 2, jpkm1
208                  DO jj = 2, jpjm1
209                     DO ji = fs_2, fs_jpim1   ! vector opt.
210                        zavi = ah_wslp2(ji,jj,jk)                ! vertical mixing coef. due to lateral mixing
211                        ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk)         ! vertical mixing coef. due to lateral mixing
212                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on temperature)
213                     END DO
214                  END DO
215               END DO
216            ELSE IF( l_traldf_rot ) THEN
217               DO jk = 2, jpkm1
218                  DO jj = 2, jpjm1
219                     DO ji = fs_2, fs_jpim1   ! vector opt.
220                        zavi = fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
221                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
222                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
223                        zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi              ! zwt=avt+zavi (total vertical mixing coef. on salinity)
224                     END DO
225                  END DO
226               END DO
227            ELSE                         ! no rotation but key_ldfslp defined
228               zwt(:,:,:) = fsavs(:,:,:)
229            ENDIF
[255]230#else
[2528]231            ! No isopycnal diffusion
232            zwt(:,:,:) = fsavs(:,:,:)           
[592]233#endif
[2528]234            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
235            DO jk = 1, jpkm1
236               DO jj = 2, jpjm1
237                  DO ji = fs_2, fs_jpim1   ! vector opt.
238                     ze3ta =  ( 1. - znvvl ) +        znvvl   * fse3t_a(ji,jj,jk)   ! after scale factor at T-point
239                     ze3tn =         znvvl   + ( 1. - znvvl ) * fse3t_n(ji,jj,jk)   ! now   scale factor at T-point
240                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
241                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
242                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
243                 END DO
244               END DO
[3]245            END DO
[2528]246            ! Surface boudary conditions
247            DO jj = 2, jpjm1
248               DO ji = fs_2, fs_jpim1   ! vector opt.
249                 ze3ta = ( 1. - znvvl ) +  znvvl * fse3t_a(ji,jj,1)    ! after scale factor at T-point
250                 zwi(ji,jj,1) = 0._wp
251                 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
252               END DO
[3]253            END DO
[2528]254            !
255            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
256            DO jj = 2, jpjm1
257               DO ji = fs_2, fs_jpim1
258                  zwt(ji,jj,1) = zwd(ji,jj,1)
259               END DO
[457]260            END DO
[2528]261            DO jk = 2, jpkm1
262               DO jj = 2, jpjm1
263                  DO ji = fs_2, fs_jpim1
264                    zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
265                  END DO
266               END DO
267            END DO
268            !
269         END IF 
[3]270
[2528]271         ! II.1. Vertical diffusion on tracer
272         ! ---------------------------
273         
274         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
[457]275         DO jj = 2, jpjm1
276            DO ji = fs_2, fs_jpim1
[2528]277               ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)
278               ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)
279               pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn)
[457]280            END DO
281         END DO
[2528]282         DO jk = 2, jpkm1
283            DO jj = 2, jpjm1
284               DO ji = fs_2, fs_jpim1
285                  ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)
286                  ze3tn = ( 1. - znvvl ) + znvvl * fse3t  (ji,jj,jk)
287                  zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn)   ! zrhs=right hand side
288                  pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)
289               END DO
[3]290            END DO
291         END DO
[457]292
[2528]293         ! third recurrence:    Xk = (Zk - Sk Xk+1 ) / Tk
[457]294         DO jj = 2, jpjm1
295            DO ji = fs_2, fs_jpim1
[2528]296               pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
[3]297            END DO
298         END DO
[2528]299         DO jk = jpk-2, 1, -1
300            DO jj = 2, jpjm1
301               DO ji = fs_2, fs_jpim1
302                  pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) )   &
303                     &             / zwt(ji,jj,jk) * tmask(ji,jj,jk)
304               END DO
[457]305            END DO
306         END DO
[2528]307         !
[457]308      END DO
[1438]309      !
[2590]310      IF(.NOT. wrk_release(3, 1,2))THEN
311         CALL ctl_stop('tra_zdf_imp : failed to release workspace arrays.')
312      END IF
313      !
[3]314   END SUBROUTINE tra_zdf_imp
315
316   !!==============================================================================
317END MODULE trazdf_imp
Note: See TracBrowser for help on using the repository browser.