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.
trdmod.F90 in trunk/NEMO/OPA_SRC/TRD – NEMO

source: trunk/NEMO/OPA_SRC/TRD/trdmod.F90 @ 1129

Last change on this file since 1129 was 1129, checked in by ctlod, 16 years ago

trunk: compilation error when key_trddyn is active related to the dynamics advection flux formulation, see ticket: #211
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 18.8 KB
Line 
1MODULE trdmod
2   !!======================================================================
3   !!                       ***  MODULE  trdmod  ***
4   !! Ocean diagnostics:  ocean tracers and dynamic trends
5   !!=====================================================================
6   !! History :  9.0  !  04-08  (C. Talandier) Original code
7   !!                 !  05-04  (C. Deltel)    Add Asselin trend in the ML budget
8   !!----------------------------------------------------------------------
9#if  defined key_trdtra || defined key_trddyn || defined key_trdmld || defined key_trdvor || defined key_esopa
10   !!----------------------------------------------------------------------
11   !!   trd_mod          : Call the trend to be computed
12   !!   trd_mod_init     : Initialization step
13   !!----------------------------------------------------------------------
14   USE oce                     ! ocean dynamics and tracers variables
15   USE dom_oce                 ! ocean space and time domain variables
16   USE zdf_oce                 ! ocean vertical physics variables
17   USE trdmod_oce              ! ocean variables trends
18   USE ldftra_oce              ! ocean active tracers lateral physics
19   USE sbc_oce                 ! surface boundary condition: ocean
20   USE phycst                  ! physical constants
21   USE trdvor                  ! ocean vorticity trends
22   USE trdicp                  ! ocean bassin integral constraints properties
23   USE trdmld                  ! ocean active mixed layer tracers trends
24   USE in_out_manager          ! I/O manager
25
26   IMPLICIT NONE
27   PRIVATE
28
29   REAL(wp) ::   r2dt          ! time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=0
30
31   PUBLIC trd_mod              ! called by all dynXX or traXX modules
32   PUBLIC trd_mod_init         ! called by opa.F90 module
33
34   !! * Substitutions
35#  include "domzgr_substitute.h90"
36#  include "vectopt_loop_substitute.h90"
37   !!----------------------------------------------------------------------
38   !!   OPA 9.0 , LOCEAN-IPSL (2005)
39   !! $Id$
40   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42
43CONTAINS
44
45   SUBROUTINE trd_mod( ptrdx, ptrdy, ktrd, ctype, kt, cnbpas )
46      !!---------------------------------------------------------------------
47      !!                  ***  ROUTINE trd_mod  ***
48      !!
49      !! ** Purpose : Dispatch all trends computation, e.g. vorticity, mld or
50      !!              integral constraints
51      !!----------------------------------------------------------------------
52      INTEGER, INTENT( in ) ::   kt                                ! time step
53      INTEGER, INTENT( in ) ::   ktrd                              ! tracer trend index
54      CHARACTER(len=3), INTENT( in ) ::   ctype                    ! momentum or tracers trends type 'DYN'/'TRA'
55      CHARACTER(len=3), INTENT( in ), OPTIONAL ::   cnbpas         ! number of passage
56      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   ptrdx ! Temperature or U trend
57      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   ptrdy ! Salinity    or V trend
58      !!
59      INTEGER ::   ji, ikbu, ikbum1
60      INTEGER ::   jj, ikbv, ikbvm1
61      CHARACTER(len=3) ::   ccpas                                  ! number of passage
62      REAL(wp) ::   zua, zva                                       ! scalars
63      REAL(wp), DIMENSION(jpi,jpj) ::   ztswu, ztswv               ! 2D workspace
64      REAL(wp), DIMENSION(jpi,jpj) ::   ztbfu, ztbfv               ! 2D workspace
65      REAL(wp), DIMENSION(jpi,jpj) ::   z2dx, z2dy                 ! workspace arrays
66      !!----------------------------------------------------------------------
67
68      z2dx(:,:)   = 0.e0   ;   z2dy(:,:)   = 0.e0                  ! initialization of workspace arrays
69
70      ! Control of optional arguments
71      ccpas = 'fst'
72      IF( PRESENT(cnbpas) )  ccpas = cnbpas
73
74      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restarting with Euler time stepping)
75      ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog)
76      ENDIF
77
78      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
79      ! I. Integral Constraints Properties for momentum and/or tracers trends
80      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
81
82      IF( ( mod(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend) )   THEN
83         !
84         IF( lk_trdtra .AND. ctype == 'TRA' )   THEN       ! active tracer trends
85            SELECT CASE ( ktrd )
86            CASE ( jptra_trd_ldf )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_ldf, ctype )   ! lateral diff
87            CASE ( jptra_trd_zdf )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_zdf, ctype )   ! vertical diff (Kz)
88            CASE ( jptra_trd_bbc )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_bbc, ctype )   ! bottom boundary cond
89            CASE ( jptra_trd_bbl )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_bbl, ctype )   ! bottom boundary layer
90            CASE ( jptra_trd_npc )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_npc, ctype )   ! static instability mixing
91            CASE ( jptra_trd_dmp )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_dmp, ctype )   ! damping
92            CASE ( jptra_trd_qsr )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_qsr, ctype )   ! penetrative solar radiat.
93            CASE ( jptra_trd_nsr )   
94               z2dx(:,:) = ptrdx(:,:,1)   ;   z2dy(:,:) = ptrdy(:,:,1)
95               CALL trd_icp( z2dx, z2dy, jpicpt_nsr, ctype )                               ! non solar radiation
96            CASE ( jptra_trd_xad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_xad, ctype )   ! x- horiz adv
97            CASE ( jptra_trd_yad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_yad, ctype )   ! y- horiz adv
98            CASE ( jptra_trd_zad )                                                         ! z- vertical adv
99               CALL trd_icp( ptrdx, ptrdy, jpicpt_zad, ctype, clpas=ccpas )   
100               ! compute the surface flux condition wn(:,:,1)*tn(:,:,1)
101               z2dx(:,:) = wn(:,:,1)*tn(:,:,1)/fse3t(:,:,1)
102               z2dy(:,:) = wn(:,:,1)*sn(:,:,1)/fse3t(:,:,1)
103               CALL trd_icp( z2dx , z2dy , jpicpt_zl1, ctype )                             ! 1st z- vertical adv
104            END SELECT
105         END IF
106
107         IF( lk_trddyn .AND. ctype == 'DYN' )   THEN       ! momentum trends
108            !
109            SELECT CASE ( ktrd )
110            CASE ( jpdyn_trd_hpg )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_hpg, ctype )   ! hydrost. pressure grad
111            CASE ( jpdyn_trd_keg )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_keg, ctype )   ! KE gradient
112            CASE ( jpdyn_trd_rvo )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_rvo, ctype )   ! relative vorticity
113            CASE ( jpdyn_trd_pvo )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_pvo, ctype )   ! planetary vorticity
114            CASE ( jpdyn_trd_ldf )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_ldf, ctype )   ! lateral diffusion
115            CASE ( jpdyn_trd_had )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_had, ctype )   ! horizontal advection
116            CASE ( jpdyn_trd_zad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_zad, ctype )   ! vertical advection
117            CASE ( jpdyn_trd_spg )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_spg, ctype )   ! surface pressure grad.
118            CASE ( jpdyn_trd_dat )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_dat, ctype )   ! damping term
119            CASE ( jpdyn_trd_zdf )                                                         ! vertical diffusion
120               ! subtract surface forcing/bottom friction trends
121               ! from vertical diffusive momentum trends
122               ztswu(:,:) = 0.e0   ;   ztswv(:,:) = 0.e0
123               ztbfu(:,:) = 0.e0   ;   ztbfv(:,:) = 0.e0 
124               DO jj = 2, jpjm1   
125                  DO ji = fs_2, fs_jpim1   ! vector opt.
126                     ! save the surface forcing momentum fluxes
127                     ztswu(ji,jj) = utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 )
128                     ztswv(ji,jj) = vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 )
129                     ! save bottom friction momentum fluxes
130                     ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
131                     ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) )
132                     ikbum1 = MAX( ikbu-1, 1 )
133                     ikbvm1 = MAX( ikbv-1, 1 )
134                     zua = ua(ji,jj,ikbum1) * r2dt + ub(ji,jj,ikbum1)
135                     zva = va(ji,jj,ikbvm1) * r2dt + vb(ji,jj,ikbvm1)
136                     ztbfu(ji,jj) = - avmu(ji,jj,ikbu) * zua / ( fse3u(ji,jj,ikbum1)*fse3uw(ji,jj,ikbu) )
137                     ztbfv(ji,jj) = - avmv(ji,jj,ikbv) * zva / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) )
138                     !
139                     ptrdx(ji,jj,1     ) = ptrdx(ji,jj,1     ) - ztswu(ji,jj)
140                     ptrdy(ji,jj,1     ) = ptrdy(ji,jj,1     ) - ztswv(ji,jj)
141                     ptrdx(ji,jj,ikbum1) = ptrdx(ji,jj,ikbum1) - ztbfu(ji,jj)
142                     ptrdy(ji,jj,ikbvm1) = ptrdy(ji,jj,ikbvm1) - ztbfv(ji,jj)
143                  END DO
144               END DO
145               !
146               CALL trd_icp( ptrdx, ptrdy, jpicpd_zdf, ctype )   
147               CALL trd_icp( ztswu, ztswv, jpicpd_swf, ctype )                               ! wind stress forcing term
148               CALL trd_icp( ztbfu, ztbfv, jpicpd_bfr, ctype )                               ! bottom friction term
149            END SELECT
150            !
151         END IF
152         !
153      END IF
154
155      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
156      ! II. Vorticity trends
157      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
158
159      IF( lk_trdvor .AND. ctype == 'DYN' )   THEN
160         !
161         SELECT CASE ( ktrd ) 
162         CASE ( jpdyn_trd_hpg )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_prg )   ! Hydrostatique Pressure Gradient
163         CASE ( jpdyn_trd_keg )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_keg )   ! KE Gradient
164         CASE ( jpdyn_trd_rvo )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_rvo )   ! Relative Vorticity
165         CASE ( jpdyn_trd_pvo )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_pvo )   ! Planetary Vorticity Term
166         CASE ( jpdyn_trd_ldf )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_ldf )   ! Horizontal Diffusion
167         CASE ( jpdyn_trd_had )   ;   CALL ctl_warn('Vorticity for horizontal advection trend never checked')   
168         CASE ( jpdyn_trd_zad )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_zad )   ! Vertical Advection
169         CASE ( jpdyn_trd_spg )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_spg )   ! Surface Pressure Grad.
170         CASE ( jpdyn_trd_dat )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_bev )   ! Beta V 
171         CASE ( jpdyn_trd_zdf )                                                      ! Vertical Diffusion
172            ! subtract surface forcing/bottom friction trends
173            ! from vertical diffusive momentum trends
174            ztswu(:,:) = 0.e0   ;   ztswv(:,:) = 0.e0
175            ztbfu(:,:) = 0.e0   ;   ztbfv(:,:) = 0.e0 
176            DO jj = 2, jpjm1   
177               DO ji = fs_2, fs_jpim1   ! vector opt.
178                  ! save the surface forcing momentum fluxes
179                  ztswu(ji,jj) = utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 )
180                  ztswv(ji,jj) = vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 )
181                  ! save bottom friction momentum fluxes
182                  ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
183                  ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) )
184                  ikbum1 = MAX( ikbu-1, 1 )
185                  ikbvm1 = MAX( ikbv-1, 1 )
186                  zua = ua(ji,jj,ikbum1) * r2dt + ub(ji,jj,ikbum1)
187                  zva = va(ji,jj,ikbvm1) * r2dt + vb(ji,jj,ikbvm1)
188                  ztbfu(ji,jj) = - avmu(ji,jj,ikbu) * zua / ( fse3u(ji,jj,ikbum1)*fse3uw(ji,jj,ikbu) )
189                  ztbfv(ji,jj) = - avmv(ji,jj,ikbv) * zva / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) )
190                  !
191                  ptrdx(ji,jj,1     ) = ptrdx(ji,jj,1     ) - ztswu(ji,jj)
192                  ptrdx(ji,jj,ikbum1) = ptrdx(ji,jj,ikbum1) - ztbfu(ji,jj)
193                  ptrdy(ji,jj,1     ) = ptrdy(ji,jj,1     ) - ztswv(ji,jj)
194                  ptrdy(ji,jj,ikbvm1) = ptrdy(ji,jj,ikbvm1) - ztbfv(ji,jj)
195               END DO
196            END DO
197            !
198            CALL trd_vor_zint( ptrdx, ptrdy, jpvor_zdf )   
199            CALL trd_vor_zint( ztswu, ztswv, jpvor_swf )                               ! Wind stress forcing term
200            CALL trd_vor_zint( ztbfu, ztbfv, jpvor_bfr )                               ! Bottom friction term
201         END SELECT
202         !
203      ENDIF
204
205      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
206      ! III. Mixed layer trends for active tracers
207      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
208
209      IF( lk_trdmld .AND. ctype == 'TRA' )   THEN
210         
211         !-----------------------------------------------------------------------------------------------
212         ! W.A.R.N.I.N.G :
213         ! jptra_trd_ldf : called by traldf.F90
214         !                 at this stage we store:
215         !                  - the lateral geopotential diffusion (here, lateral = horizontal)
216         !                  - and the iso-neutral diffusion if activated
217         ! jptra_trd_zdf : called by trazdf.F90
218         !                 * in case of iso-neutral diffusion we store the vertical diffusion component in the
219         !                   lateral trend including the K_z contrib, which will be removed later (see trd_mld)
220         !-----------------------------------------------------------------------------------------------
221
222         SELECT CASE ( ktrd )
223         CASE ( jptra_trd_xad )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_xad, '3D' )   ! merid. advection
224         CASE ( jptra_trd_yad )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_yad, '3D' )   ! zonal  advection
225         CASE ( jptra_trd_zad )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_zad, '3D' )   ! vertical advection
226         CASE ( jptra_trd_ldf )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_ldf, '3D' )   ! lateral diffusive
227         CASE ( jptra_trd_bbl )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_bbl, '3D' )   ! bottom boundary layer
228         CASE ( jptra_trd_zdf )
229            IF( ln_traldf_iso )   THEN
230               CALL trd_mld_zint( ptrdx, ptrdy, jpmld_ldf, '3D' )   ! vertical diffusion (K_z)
231            ELSE
232               CALL trd_mld_zint( ptrdx, ptrdy, jpmld_zdf, '3D' )   ! vertical diffusion (K_z)
233            ENDIF
234         CASE ( jptra_trd_dmp )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_dmp, '3D' )   ! internal 3D restoring (tradmp)
235         CASE ( jptra_trd_qsr )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_for, '3D' )   ! air-sea : penetrative sol radiat
236         CASE ( jptra_trd_nsr )
237            ptrdx(:,:,2:jpk) = 0.e0   ;   ptrdy(:,:,2:jpk) = 0.e0
238            CALL trd_mld_zint( ptrdx, ptrdy, jpmld_for, '2D' )                             ! air-sea : non penetr sol radiat
239         CASE ( jptra_trd_bbc )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_bbc, '3D' )   ! bottom bound cond (geoth flux)
240         CASE ( jptra_trd_atf )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_atf, '3D' )   ! asselin numerical
241         CASE ( jptra_trd_npc )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_npc, '3D' )   ! non penetr convect adjustment
242         END SELECT
243
244      ENDIF
245
246   END SUBROUTINE trd_mod
247
248#   else
249   !!----------------------------------------------------------------------
250   !!   Default case :                                         Empty module
251   !!----------------------------------------------------------------------
252   USE trdmod_oce      ! ocean variables trends
253   USE trdvor          ! ocean vorticity trends
254   USE trdicp          ! ocean bassin integral constraints properties
255   USE trdmld          ! ocean active mixed layer tracers trends
256
257CONTAINS
258   SUBROUTINE trd_mod(ptrd3dx, ptrd3dy, ktrd , ctype, kt, cnbpas)   ! Empty routine
259      REAL, DIMENSION(:,:,:), INTENT( in ) ::   &
260          ptrd3dx,                     &                           ! Temperature or U trend
261          ptrd3dy                                                  ! Salinity    or V trend
262      INTEGER, INTENT( in ) ::   ktrd                              ! momentum or tracer trend index
263      INTEGER, INTENT( in ) ::   kt                                ! Time step
264      CHARACTER(len=3), INTENT( in ) ::  ctype                     ! momentum or tracers trends type
265      CHARACTER(len=3), INTENT( in ), OPTIONAL ::   cnbpas         ! number of passage
266      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', ptrd3dx(1,1,1)
267      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ptrd3dy(1,1,1)
268      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ktrd
269      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ctype
270      WRITE(*,*) ' "   ": You should not have seen this print! error ?', kt
271      WRITE(*,*) ' "   ": You should not have seen this print! error ?', cnbpas
272   END SUBROUTINE trd_mod
273#   endif
274
275   SUBROUTINE trd_mod_init
276      !!----------------------------------------------------------------------
277      !!                  ***  ROUTINE trd_mod_init  ***
278      !!
279      !! ** Purpose :   Initialization of activated trends
280      !!----------------------------------------------------------------------
281      USE in_out_manager          ! I/O manager
282
283      NAMELIST/namtrd/ ntrd, nctls, ln_trdmld_restart, ucf, ln_trdmld_instant
284      !!----------------------------------------------------------------------
285
286      IF( l_trdtra .OR. l_trddyn )   THEN
287         REWIND( numnam )
288         READ  ( numnam, namtrd )      ! namelist namtrd : trends diagnostic
289
290         IF(lwp) THEN
291            WRITE(numout,*)
292            WRITE(numout,*) ' trd_mod_init : Momentum/Tracers trends'
293            WRITE(numout,*) ' ~~~~~~~~~~~~~'
294            WRITE(numout,*) '       Namelist namtrd : set trends parameters'
295            WRITE(numout,*) '           * frequency of trends diagnostics   ntrd               = ', ntrd
296            WRITE(numout,*) '           * control surface type              nctls              = ', nctls
297            WRITE(numout,*) '           * restart for ML diagnostics        ln_trdmld_restart  = ', ln_trdmld_restart
298            WRITE(numout,*) '           * instantaneous or mean ML T/S      ln_trdmld_instant  = ', ln_trdmld_instant
299            WRITE(numout,*) '           * unit conversion factor            ucf                = ', ucf
300        ENDIF
301      ENDIF
302      !
303      IF( lk_trddyn .OR. lk_trdtra )    CALL trd_icp_init       ! integral constraints trends
304      IF( lk_trdmld                )    CALL trd_mld_init       ! mixed-layer trends (active  tracers) 
305      IF( lk_trdvor                )    CALL trd_vor_init       ! vorticity trends       
306      !
307   END SUBROUTINE trd_mod_init
308
309   !!======================================================================
310END MODULE trdmod
Note: See TracBrowser for help on using the repository browser.