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.
trdtra.F90 in branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: branches/2012/dev_r3309_LOCEAN12_Ediag/NEMOGCM/NEMO/OPA_SRC/TRD/trdtra.F90 @ 3325

Last change on this file since 3325 was 3318, checked in by gm, 12 years ago

Ediag branche: #927 split TRA/DYN trd computation
  • Property svn:keywords set to Id
File size: 19.4 KB
Line 
1MODULE trdtra
2   !!======================================================================
3   !!                       ***  MODULE  trdtra  ***
4   !! Ocean diagnostics:  ocean tracers trends pre-processing
5   !!=====================================================================
6   !! History :  3.3  !  2010-06  (C. Ethe) creation for the TRA/TRC merge
7   !!            3.5  !  2012-02  (G. Madec) update the comments
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   trd_tra       : pre-process the tracer trends
12   !!   trd_tra_adv   : transform a div(U.T) trend into a U.grad(T) trend
13   !!   trd_tra_mng   : tracer trend manager: dispatch to the diagnostic modules
14   !!   trd_tra_iom   : output 3D tracer trends using IOM
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and tracers variables
17   USE dom_oce        ! ocean domain
18   USE sbc_oce        ! surface boundary condition: ocean
19   USE zdf_oce        ! ocean vertical physics
20   USE trd_oce        ! trends: ocean variables
21   USE trdmod_trc     ! ocean passive mixed layer tracers trends
22   USE trdglo         ! trends:global domain averaged
23   USE trdmld         ! ocean active mixed layer tracers trends
24   USE ldftra_oce     ! ocean active tracers lateral physics
25   USE zdfddm         ! vertical physics: double diffusion
26   USE phycst         ! physical constants
27   USE in_out_manager ! I/O manager
28   USE iom            ! I/O manager library
29   USE lib_mpp        ! MPP library
30   USE wrk_nemo       ! Memory allocation
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   trd_tra   ! called by all tra_... modules
36
37   REAL(wp) ::   r2dt   ! time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=0
38
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   trdtx, trdty, trdt   ! use to store the temperature trends
40
41   !! * Substitutions
42#  include "domzgr_substitute.h90"
43#  include "zdfddm_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   INTEGER FUNCTION trd_tra_alloc()
53      !!---------------------------------------------------------------------
54      !!                  ***  FUNCTION trd_tra_alloc  ***
55      !!---------------------------------------------------------------------
56      ALLOCATE( trdtx(jpi,jpj,jpk) , trdty(jpi,jpj,jpk) , trdt(jpi,jpj,jpk) , STAT= trd_tra_alloc )
57      !
58      IF( lk_mpp             )   CALL mpp_sum ( trd_tra_alloc )
59      IF( trd_tra_alloc /= 0 )   CALL ctl_warn('trd_tra_alloc: failed to allocate arrays')
60   END FUNCTION trd_tra_alloc
61
62
63   SUBROUTINE trd_tra( kt, ctype, ktra, ktrd, ptrd, pun, ptra )
64      !!---------------------------------------------------------------------
65      !!                  ***  ROUTINE trd_tra  ***
66      !!
67      !! ** Purpose : pre-process tracer trends
68      !!
69      !! ** Method  : - mask the trend
70      !!              - advection (ptra present) converte the incoming flux (U.T)
71      !!              into trend (U.T => -U.grat(T)=div(U.T)-T.div(U)) through a
72      !!              call to trd_tra_adv
73      !!              - 'TRA' case : regroup T & S trends
74      !!              - send the trends to trd_tra_mng (trdmod_trc) for further processing
75      !!----------------------------------------------------------------------
76      INTEGER                         , INTENT(in)           ::   kt      ! time step
77      CHARACTER(len=3)                , INTENT(in)           ::   ctype   ! tracers trends type 'TRA'/'TRC'
78      INTEGER                         , INTENT(in)           ::   ktra    ! tracer index
79      INTEGER                         , INTENT(in)           ::   ktrd    ! tracer trend index
80      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in)           ::   ptrd    ! tracer trend  or flux
81      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   pun     ! now velocity
82      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in), OPTIONAL ::   ptra    ! now tracer variable
83      !
84      INTEGER  ::   jk   ! loop indices
85      REAL(wp), POINTER, DIMENSION(:,:,:)  ::   zwt, zws, ztrdt, ztrds   ! 3D workspace
86      !!----------------------------------------------------------------------
87      !
88      CALL wrk_alloc( jpi, jpj, jpk, ztrds )
89      !     
90      IF( .NOT. ALLOCATED( trdtx ) ) THEN      ! allocate trdtra arrays
91         IF( trd_tra_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trd_tra : unable to allocate arrays' )
92      ENDIF
93
94      IF( ctype == 'TRA' .AND. ktra == jp_tem ) THEN   !==  Temperature trend  ==!
95         !
96         IF( PRESENT( ptra ) ) THEN                       ! advection: transform flux into trend
97            SELECT CASE( ktrd )     
98            CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd, pun, ptra, 'X', trdtx ) 
99            CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd, pun, ptra, 'Y', trdty ) 
100            CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd, pun, ptra, 'Z', trdt  ) 
101            END SELECT
102         ELSE                                             ! other trends:
103            trdt(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)              ! mask & store
104            IF( ktrd == jptra_bbc .OR. ktrd == jptra_qsr ) THEN   ! qsr, bbc: on temperature only
105               ztrds(:,:,:) = 0._wp
106               CALL trd_tra_mng( trdt, ztrds, ktrd, kt )               ! send to trd_tra_mng
107            ENDIF
108         ENDIF
109         !
110      ENDIF
111
112      IF( ctype == 'TRA' .AND. ktra == jp_sal ) THEN      !==  Salinity trends  ==!
113         !
114         IF( PRESENT( ptra ) ) THEN      ! advection: transform the advective flux into a trend
115            SELECT CASE( ktrd )          !            and send T & S trends to trd_tra_mng
116            CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd , pun  , ptra, 'X'  , ztrds ) 
117                                    CALL trd_tra_mng( trdtx, ztrds, ktrd, kt   )
118            CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd , pun  , ptra, 'Y'  , ztrds ) 
119                                ;   CALL trd_tra_mng( trdty, ztrds, ktrd, kt   )
120            CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd , pun  , ptra, 'Z'  , ztrds ) 
121                                    CALL trd_tra_mng( trdt , ztrds, ktrd, kt   )
122            END SELECT
123         ELSE                            ! other trends: mask and send T & S trends to trd_tra_mng
124            ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
125            CALL trd_tra_mng( trdt, ztrds, ktrd, kt ) 
126         ENDIF
127         !
128         IF( ktrd == jptra_zdfp ) THEN     ! diagnose the "PURE" Kz trend (here: just before the swap)
129            !
130            IF( ln_traldf_iso ) THEN      ! iso-neutral diffusion only otherwise jptra_zdf is "PURE"
131               !
132               CALL wrk_alloc( jpi, jpj, jpk, zwt, zws, ztrdt )
133               !
134               zwt(:,:, 1 ) = 0._wp   ;   zws(:,:, 1 ) = 0._wp            ! vertical diffusive fluxes
135               zwt(:,:,jpk) = 0._wp   ;   zws(:,:,jpk) = 0._wp
136               DO jk = 2, jpk
137                  zwt(:,:,jk) =   avt(:,:,jk) * ( tsa(:,:,jk-1,jp_tem) - tsa(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk)
138                  zws(:,:,jk) = fsavs(:,:,jk) * ( tsa(:,:,jk-1,jp_sal) - tsa(:,:,jk,jp_sal) ) / fse3w(:,:,jk) * tmask(:,:,jk)
139               END DO
140               !
141               ztrdt(:,:,jpk) = 0._wp   ;   ztrds(:,:,jpk) = 0._wp
142               DO jk = 1, jpkm1
143                  ztrdt(:,:,jk) = ( zwt(:,:,jk) - zwt(:,:,jk+1) ) / fse3t(:,:,jk)
144                  ztrds(:,:,jk) = ( zws(:,:,jk) - zws(:,:,jk+1) ) / fse3t(:,:,jk) 
145               END DO
146               CALL trd_tra_mng( ztrdt, ztrds, jptra_zdfp, kt ) 
147               !
148               CALL wrk_dealloc( jpi, jpj, jpk, zwt, zws, ztrdt )
149               !
150            ENDIF
151            !
152         ENDIF
153      ENDIF
154
155      IF( ctype == 'TRC' ) THEN                           !==  passive tracer trend  ==!
156         !
157         IF( PRESENT( ptra ) ) THEN                          ! advection: transform flux into a masked trend
158            SELECT CASE( ktrd )
159            CASE( jptra_xad )   ;   CALL trd_tra_adv( ptrd , pun , ptra, 'X', ztrds ) 
160            CASE( jptra_yad )   ;   CALL trd_tra_adv( ptrd , pun , ptra, 'Y', ztrds ) 
161            CASE( jptra_zad )   ;   CALL trd_tra_adv( ptrd , pun , ptra, 'Z', ztrds ) 
162            END SELECT
163         ELSE                                                ! other trends: masked trend
164            ztrds(:,:,:) = ptrd(:,:,:) * tmask(:,:,:)
165         END IF
166         !                                 
167         CALL trd_mod_trc( ztrds, ktra, ktrd, kt )           ! send trend to trd_mod_trc
168         !
169      ENDIF
170      !
171      CALL wrk_dealloc( jpi, jpj, jpk, ztrds )
172      !
173   END SUBROUTINE trd_tra
174
175
176   SUBROUTINE trd_tra_adv( pf, pun, ptn, cdir, ptrd )
177      !!---------------------------------------------------------------------
178      !!                  ***  ROUTINE trd_tra_adv  ***
179      !!
180      !! ** Purpose :   transformed a advective flux into a masked advective trends
181      !!
182      !! ** Method  :   use the following transformation: -div(U.T) = - U grad(T) + T.div(U)
183      !!       i-advective trends = -un. di-1[T] = -( di-1[fi] - tn di-1[un] )
184      !!       j-advective trends = -un. di-1[T] = -( dj-1[fi] - tn dj-1[un] )
185      !!       k-advective trends = -un. di+1[T] = -( dk+1[fi] - tn dk+1[un] )
186      !!                where fi is the incoming advective flux.
187      !!----------------------------------------------------------------------
188      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pf      ! advective flux in one direction
189      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pun     ! now velocity   in one direction
190      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   ptn     ! now or before tracer
191      CHARACTER(len=1)                , INTENT(in   ) ::   cdir    ! X/Y/Z direction
192      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   ptrd    ! advective trend in one direction
193      !
194      INTEGER  ::   ji, jj, jk   ! dummy loop indices
195      INTEGER  ::   ii, ij, ik   ! index shift as function of the direction
196      !!----------------------------------------------------------------------
197      !
198      SELECT CASE( cdir )      ! shift depending on the direction
199      CASE( 'X' )   ;   ii = 1   ;   ij = 0   ;   ik = 0      ! i-trend
200      CASE( 'Y' )   ;   ii = 0   ;   ij = 1   ;   ik = 0      ! j-trend
201      CASE( 'Z' )   ;   ii = 0   ;   ij = 0   ;   ik =-1      ! k-trend
202      END SELECT
203      !
204      !                        ! set to zero uncomputed values
205      ptrd(jpi,:,:) = 0._wp   ;   ptrd(1,:,:) = 0._wp
206      ptrd(:,jpj,:) = 0._wp   ;   ptrd(:,1,:) = 0._wp
207      ptrd(:,:,jpk) = 0._wp
208      !
209      DO jk = 1, jpkm1         ! advective trend
210         DO jj = 2, jpjm1
211            DO ji = fs_2, fs_jpim1   ! vector opt.
212               ptrd(ji,jj,jk) = - (     pf (ji,jj,jk) - pf (ji-ii,jj-ij,jk-ik)                        &
213                 &                  - ( pun(ji,jj,jk) - pun(ji-ii,jj-ij,jk-ik) ) * ptn(ji,jj,jk)  )   &
214                 &              / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )  * tmask(ji,jj,jk)
215            END DO
216         END DO
217      END DO
218      !
219   END SUBROUTINE trd_tra_adv
220
221
222   SUBROUTINE trd_tra_mng( ptrdx, ptrdy, ktrd, kt )
223      !!---------------------------------------------------------------------
224      !!                  ***  ROUTINE trd_tra_mng  ***
225      !!
226      !! ** Purpose :   Dispatch all trends computation, e.g. 3D output, integral
227      !!                constraints, barotropic vorticity, kinetic enrgy,
228      !!                potential energy, and/or mixed layer budget.
229      !!----------------------------------------------------------------------
230      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend
231      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
232      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
233      INTEGER                   , INTENT(in   ) ::   kt      ! time step
234      !!----------------------------------------------------------------------
235
236      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restart with Euler time stepping)
237      ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog)
238      ENDIF
239
240      !                   ! 3D output of tracers trends using IOM interface
241      IF( ln_tra_trd )   CALL trd_tra_iom ( ptrdx, ptrdy, ktrd, kt )
242
243      !                   ! Integral Constraints Properties for tracers trends                                       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
244      IF( ln_glo_trd )   CALL trd_glo( ptrdx, ptrdy, ktrd, 'TRA', kt )
245
246      !                   ! Mixed layer trends for active tracers
247      IF( ln_tra_mld )   THEN   
248         !-----------------------------------------------------------------------------------------------
249         ! W.A.R.N.I.N.G :
250         ! jptra_ldf : called by traldf.F90
251         !                 at this stage we store:
252         !                  - the lateral geopotential diffusion (here, lateral = horizontal)
253         !                  - and the iso-neutral diffusion if activated
254         ! jptra_zdf : called by trazdf.F90
255         !                 * in case of iso-neutral diffusion we store the vertical diffusion component in the
256         !                   lateral trend including the K_z contrib, which will be removed later (see trd_mld)
257         !-----------------------------------------------------------------------------------------------
258
259         SELECT CASE ( ktrd )
260         CASE ( jptra_xad )        ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_xad, '3D' )   ! zonal    advection
261         CASE ( jptra_yad )        ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_yad, '3D' )   ! merid.   advection
262         CASE ( jptra_zad )        ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_zad, '3D' )   ! vertical advection
263         CASE ( jptra_ldf )        ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_ldf, '3D' )   ! lateral  diffusion
264         CASE ( jptra_bbl )        ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_bbl, '3D' )   ! bottom boundary layer
265         CASE ( jptra_zdf )
266            IF( ln_traldf_iso ) THEN ; CALL trd_mld_zint( ptrdx, ptrdy, jpmld_ldf, '3D' )   ! lateral  diffusion (K_z)
267            ELSE                   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_zdf, '3D' )   ! vertical diffusion (K_z)
268            ENDIF
269         CASE ( jptra_dmp )        ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_dmp, '3D' )   ! internal 3D restoring (tradmp)
270         CASE ( jptra_qsr )        ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_for, '3D' )   ! air-sea : penetrative sol radiat
271         CASE ( jptra_nsr )        ;   ptrdx(:,:,2:jpk) = 0._wp   ;   ptrdy(:,:,2:jpk) = 0._wp
272                                   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_for, '2D' )   ! air-sea : non penetr sol radiation
273         CASE ( jptra_bbc )        ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_bbc, '3D' )   ! bottom bound cond (geoth flux)
274         CASE ( jptra_npc )        ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_npc, '3D' )   ! non penetr convect adjustment
275         CASE ( jptra_atf )        ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_atf, '3D' )   ! asselin time filter (last trend)
276                                   !
277                                       CALL trd_mld( kt )                                   ! trends: Mixed-layer (output)
278         END SELECT
279         !
280      ENDIF
281      !
282   END SUBROUTINE trd_tra_mng
283
284
285   SUBROUTINE trd_tra_iom( ptrdx, ptrdy, ktrd, kt )
286      !!---------------------------------------------------------------------
287      !!                  ***  ROUTINE trd_tra_iom  ***
288      !!
289      !! ** Purpose :   output 3D tracer trends using IOM
290      !!----------------------------------------------------------------------
291      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdx   ! Temperature or U trend
292      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptrdy   ! Salinity    or V trend
293      INTEGER                   , INTENT(in   ) ::   ktrd    ! tracer trend index
294      INTEGER                   , INTENT(in   ) ::   kt      ! time step
295      !!
296      INTEGER ::   ji, jj, jk   ! dummy loop indices
297      INTEGER ::   ikbu, ikbv   ! local integers
298      REAL(wp), POINTER, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace
299      !!----------------------------------------------------------------------
300      !
301      CALL wrk_alloc( jpi, jpj, z2dx, z2dy )
302      !
303!!gm Rq: mask the trends already masked in trd_tra, but lbc_lnk should probably be added
304      !
305      SELECT CASE( ktrd )
306      CASE( jptra_xad  )   ;   CALL iom_put( "ttrd_xad" , ptrdx )        ! x- horizontal advection
307                                   CALL iom_put( "strd_xad" , ptrdy )
308      CASE( jptra_yad  )   ;   CALL iom_put( "ttrd_yad" , ptrdx )        ! y- horizontal advection
309                                   CALL iom_put( "strd_yad" , ptrdy )
310      CASE( jptra_zad  )   ;   CALL iom_put( "ttrd_zad" , ptrdx )        ! z- vertical   advection
311                               CALL iom_put( "strd_zad" , ptrdy )
312                               IF( .NOT.lk_vvl ) THEN                   ! cst volume : adv flux through z=0 surface
313                                  z2dx(:,:) = wn(:,:,1) * tsn(:,:,1,jp_tem) / fse3t(:,:,1)
314                                  z2dy(:,:) = wn(:,:,1) * tsn(:,:,1,jp_sal) / fse3t(:,:,1)
315                                  CALL iom_put( "ttrd_sad", z2dx )
316                                  CALL iom_put( "strd_sad", z2dy )
317                               ENDIF
318      CASE( jptra_ldf  )   ;   CALL iom_put( "ttrd_ldf" , ptrdx )        ! lateral diffusion
319                               CALL iom_put( "strd_ldf" , ptrdy )
320      CASE( jptra_zdf  )   ;   CALL iom_put( "ttrd_zdf" , ptrdx )        ! vertical diffusion (including Kz contribution)
321                               CALL iom_put( "strd_zdf" , ptrdy )
322      CASE( jptra_zdfp )   ;   CALL iom_put( "ttrd_zdfp", ptrdx )        ! PURE vertical diffusion (no isoneutral contribution)
323                               CALL iom_put( "strd_zdfp", ptrdy )
324      CASE( jptra_dmp  )   ;   CALL iom_put( "ttrd_dmp" , ptrdx )        ! internal restoring (damping)
325                               CALL iom_put( "strd_dmp" , ptrdy )
326      CASE( jptra_bbl  )   ;   CALL iom_put( "ttrd_bbl" , ptrdx )        ! bottom boundary layer
327                               CALL iom_put( "strd_bbl" , ptrdy )
328      CASE( jptra_npc  )   ;   CALL iom_put( "ttrd_npc" , ptrdx )        ! static instability mixing
329                               CALL iom_put( "strd_npc" , ptrdy )
330      CASE( jptra_nsr  )   ;   CALL iom_put( "ttrd_qns" , ptrdx )        ! surface forcing + runoff (ln_rnf=T)
331                               CALL iom_put( "strd_cdt" , ptrdy )
332      CASE( jptra_qsr  )   ;   CALL iom_put( "ttrd_qsr" , ptrdx )        ! penetrative solar radiat. (only on temperature)
333      CASE( jptra_bbc  )   ;   CALL iom_put( "ttrd_bbc" , ptrdx )        ! geothermal heating   (only on temperature)
334      CASE( jptra_atf  )   ;   CALL iom_put( "ttrd_atf" , ptrdx )        ! asselin time Filter
335                               CALL iom_put( "strd_atf" , ptrdy )
336      END SELECT
337      !
338      CALL wrk_dealloc( jpi, jpj, z2dx, z2dy )
339      !
340   END SUBROUTINE trd_tra_iom
341
342   !!======================================================================
343END MODULE trdtra
Note: See TracBrowser for help on using the repository browser.