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 branches/dev_001_GM/NEMO/OPA_SRC/TRD – NEMO

source: branches/dev_001_GM/NEMO/OPA_SRC/TRD/trdmod.F90 @ 790

Last change on this file since 790 was 790, checked in by gm, 16 years ago

dev_001_GM - complete theprevious comit with omitted routines

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.7 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_tra          : active tracer trend manager
12   !!   trd_tra_adv      : pre-treatment of the tracer advection trends
13   !!   trd_mod          : momentum trend manager
14   !!   trd_mod_init     : Initialization step
15   !!----------------------------------------------------------------------
16   USE phycst                  ! physical constants
17   USE oce                     ! ocean dynamics and tracers variables
18   USE dom_oce                 ! ocean space and time domain variables
19   USE zdf_oce                 ! ocean vertical physics variables
20   USE trdmod_oce              ! ocean variables trends
21   USE ldftra_oce              ! ocean active tracers lateral physics
22   USE trdvor                  ! ocean vorticity trends
23   USE trdicp                  ! ocean bassin integral constraints properties
24   USE trdmld                  ! ocean active mixed layer tracers trends
25   USE in_out_manager          ! I/O manager
26   USE taumod                  ! surface ocean stress
27
28   IMPLICIT NONE
29   PRIVATE
30
31   REAL(wp) ::   r2dt          ! time-step, = 2 rdttra except at nit000 (=rdttra) if neuler=0
32
33   PUBLIC trd_tra              ! called by all traXXX.F90 modules
34   PUBLIC trd_tra_adv          ! called by all traadv_XXX.F90 modules
35   PUBLIC trd_mod              ! called by all dynXXX.F90 modules
36   PUBLIC trd_mod_init         ! called by opa.F90 module
37
38   !! * Substitutions
39#  include "domzgr_substitute.h90"
40#  include "vectopt_loop_substitute.h90"
41   !!----------------------------------------------------------------------
42   !!   OPA 9.0 , LOCEAN-IPSL (2005)
43   !! $Header$
44   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46
47CONTAINS
48
49   SUBROUTINE trd_tra( kt, ktra, ktrd, ctype, ptrd2d, ptrd3d, cnbpas )
50      !!---------------------------------------------------------------------
51      !!                  ***  ROUTINE trd_mod  ***
52      !!
53      !! ** Purpose : Dispatch all trends computation, e.g. vorticity, mld or
54      !!              integral constraints
55      !!----------------------------------------------------------------------
56      INTEGER         , INTENT(in   )                                   ::   kt      ! time step
57      INTEGER         , INTENT(in   )                                   ::   ktra    ! tracer index
58      INTEGER         , INTENT(in   )                                   ::   ktrd    ! tracer trend index
59      CHARACTER(len=3), INTENT(in   )                                   ::   ctype   ! tracers type 'TRA' or 'TRC'
60      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj)    , OPTIONAL ::   ptrd2d  ! Temperature or U trend
61      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk), OPTIONAL ::   ptrd3d  ! Temperature or U trend
62      CHARACTER(len=3), INTENT(in   )                        , OPTIONAL ::   cnbpas  ! number of passage
63      !!
64      CHARACTER(len=3) ::   ccpas                                  ! number of passage
65!     REAL(wp), DIMENSION(jpi,jpj) ::   z2dx, z2dy                 ! workspace arrays
66      !!----------------------------------------------------------------------
67
68      ccpas = 'fst'                            ! Control of optional arguments
69      IF( PRESENT(cnbpas) )  ccpas = cnbpas
70
71      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (restarting with Euler time stepping)
72      ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog)
73      ENDIF
74
75      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
76      ! I. Integral Constraints Properties for momentum and/or tracers trends
77      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
78
79      IF( ( mod(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend) )   THEN
80         !
81         IF( lk_trdtra .AND. ctype == 'TRA' )   THEN       ! active tracer trends
82            IF( PRESENT(ptrd2d) ) THEN   ;   CALL trd_icp( ptrd2d, ktra, ktrd, ctype, clpas=ccpas )
83            ELSE                         ;   CALL trd_icp( ptrd3d, ktra, ktrd, ctype, clpas=ccpas )
84            ENDIF
85         ENDIF
86!           SELECT CASE ( ktrd )
87!           CASE ( jptra_trd_ldf )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_ldf, ctype )   ! lateral diff
88!           CASE ( jptra_trd_zdf )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_zdf, ctype )   ! vertical diff (Kz)
89!           CASE ( jptra_trd_bbc )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_bbc, ctype )   ! bottom boundary cond
90!           CASE ( jptra_trd_bbl )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_bbl, ctype )   ! bottom boundary layer
91!           CASE ( jptra_trd_npc )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_npc, ctype )   ! static instability mixing
92!           CASE ( jptra_trd_dmp )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_dmp, ctype )   ! damping
93!           CASE ( jptra_trd_qsr )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_qsr, ctype )   ! penetrative solar radiat.
94!           CASE ( jptra_trd_nsr )   
95!              z2dx(:,:) = ptrdx(:,:,1)   ;   z2dy(:,:) = ptrdy(:,:,1)
96!              CALL trd_icp( z2dx, z2dy, jpicpt_nsr, ctype )                               ! non solar radiation
97!           CASE ( jptra_trd_xad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_xad, ctype )   ! x- horiz adv
98!           CASE ( jptra_trd_yad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_yad, ctype )   ! y- horiz adv
99!           CASE ( jptra_trd_zad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpt_zad, ctype, clpas=ccpas )    ! z- adv
100!              CALL trd_icp( ptrdx, ptrdy, jpicpt_zad, ctype, clpas=ccpas )   
101!              ! compute the surface flux condition wn(:,:,1)*tn(:,:,1)
102!              z2dx(:,:) = wn(:,:,1)*tn(:,:,1)/fse3t(:,:,1)
103!              z2dy(:,:) = wn(:,:,1)*sn(:,:,1)/fse3t(:,:,1)
104!              CALL trd_icp( z2dx , z2dy , jpicpt_zl1, ctype )                             ! 1st z- vertical adv
105!           END SELECT
106!        END IF
107         !
108      END IF
109
110      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
111      ! III. Mixed layer trends for active tracers
112      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
113
114      IF( lk_trdmld .AND. ctype == 'TRA' )   THEN
115         
116         !-----------------------------------------------------------------------------------------------
117         ! W.A.R.N.I.N.G :
118         ! jptra_trd_ldf : called by traldf.F90
119         !                 at this stage we store:
120         !                  - the lateral geopotential diffusion (here, lateral = horizontal)
121         !                  - and the iso-neutral diffusion if activated
122         ! jptra_trd_zdf : called by trazdf.F90
123         !                 * in case of iso-neutral diffusion we store the vertical diffusion component in the
124         !                   lateral trend including the K_z contrib, which will be removed later (see trd_mld)
125         !-----------------------------------------------------------------------------------------------
126
127!        SELECT CASE ( ktrd )
128!        CASE ( jpt_trd_zdf )
129!           IF( ln_traldf_iso )   THEN
130!              CALL trd_mld_zint( ptrdx, ptrdy, jpmld_ldf, '3D' )   ! vertical diffusion (K_z)
131!           ELSE
132!              CALL trd_mld_zint( ptrdx, ptrdy, jpmld_zdf, '3D' )   ! vertical diffusion (K_z)
133!           ENDIF
134!        CASE ( jpt_trd_nsr )
135!           ptrdx(:,:,2:jpk) = 0.e0   ;   ptrdy(:,:,2:jpk) = 0.e0
136!           CALL trd_mld_zint( ptrdx, ptrdy, jpmld_for, '2D' )                             ! air-sea : non solar flux
137!        CASE ( jpt_trd_bbc )   ;   CALL trd_mld_zint( ptrdx, ptrdy, jpmld_bbc, '3D' )   ! bbc (geothermal flux)
138!        CASE DEFAULT
139!                                     CALL trd_mld_zint( ptrdx, ptrdy, ktra     , '3D' )   ! other 3D trends
140!        END SELECT
141
142      ENDIF
143
144   END SUBROUTINE trd_tra
145
146
147   SUBROUTINE trd_tra_adv( kt, ktra, ktrd, ctype, pf, pun, ptn, cnbpas )
148      !!---------------------------------------------------------------------
149      !!                  ***  ROUTINE trd_mod  ***
150      !!
151      !! ** Purpose :   transformed the i-advective flux into the i-advective trends
152      !! ** Method  :   i-advective trends = un. di[T] = di[fi] - tn di[un]
153      !!----------------------------------------------------------------------
154      INTEGER         , INTENT(in   )                         ::   kt      ! time step
155      INTEGER         , INTENT(in   )                         ::   ktra    ! tracer index
156      INTEGER         , INTENT(in   )                         ::   ktrd    ! tracer trend index
157      CHARACTER(len=3), INTENT(in   )                         ::   ctype   ! tracers type 'TRA' or 'TRC'
158      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pf      ! advective flux in one direction
159      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pun     ! now velocity  in one direction
160      REAL(wp)        , INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   ptn     ! now or before tracer
161      CHARACTER(len=3), INTENT(in   )                        , OPTIONAL ::   cnbpas  ! number of passage
162      !!
163      INTEGER                          ::   ji, jj, jk      ! dummy loop indices
164      CHARACTER(len=3) ::   ccpas                           ! number of passage
165      REAL(wp)                         ::   zbtr, z_hdivn   !
166      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrdt           ! 3D workspace
167      !!----------------------------------------------------------------------
168
169      ccpas = 'fst'                            ! Control of optional arguments
170      IF( PRESENT(cnbpas) )  ccpas = cnbpas
171
172      ztrdt(:,:,:) = 0.e0
173      !
174      IF( ccpas == 'fst' ) THEN      ! first treatment : remove the divergence
175         SELECT CASE( ktrd ) 
176         CASE( jpt_trd_xad )      ! i-advective trend
177            DO jk = 1, jpkm1
178               DO jj = 2, jpjm1
179                  DO ji = fs_2, fs_jpim1   ! vector opt.
180#if defined key_zco
181                     zbtr    = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) )
182                     z_hdivn = (  e2u(ji  ,jj) * pun(ji  ,jj,jk)                              &
183                        &       - e2u(ji-1,jj) * pun(ji-1,jj,jk) )
184#else
185                     zbtr    = 1.e0/ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
186                     z_hdivn = (  e2u(ji  ,jj) * fse3u(ji  ,jj,jk) * pun(ji  ,jj,jk)          &
187                        &       - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) )
188#endif
189                     ztrdt(ji,jj,jk) = - zbtr * (  pf(ji,jj,jk) - pf(ji-1,jj,jk)  -  ptn(ji,jj,jk) * z_hdivn  )
190                  END DO
191               END DO
192            END DO
193            !
194         CASE( jpt_trd_yad )      ! j-advective trend
195            DO jk = 1, jpkm1
196               DO jj = 2, jpjm1
197                  DO ji = fs_2, fs_jpim1   ! vector opt.
198#if defined key_zco
199                     zbtr    = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) )
200                     z_hdivn = (  e1v(ji,jj  ) * pun(ji,jj  ,jk)                              &
201                        &       - e1v(ji,jj-1) * pun(ji,jj-1,jk) )
202#else
203                     zbtr    = 1.e0/ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
204                     z_hdivn = (  e1v(ji  ,jj) * fse3v(ji,jj  ,jk) * pun(ji,jj  ,jk)          &
205                        &       - e1v(ji-1,jj) * fse3v(ji,jj-1,jk) * pun(ji,jj-1,jk) )
206#endif
207                     ztrdt(ji,jj,jk) = - zbtr * (  pf(ji,jj,jk) - pf(ji,jj-1,jk)  -  ptn(ji,jj,jk) * z_hdivn  )
208                  END DO
209               END DO
210            END DO
211            !
212         CASE( jpt_trd_zad )      ! z-advective trend
213            DO jk = 1, jpkm1
214               DO jj = 2, jpjm1
215                  DO ji = fs_2, fs_jpim1   ! vector opt.
216                     zbtr    = 1.e0 / fse3t(ji,jj,jk)
217                     z_hdivn = pun(ji,jj,jk) - pun(ji,jj,jk+1)
218                     ztrdt(ji,jj,jk) = - zbtr * (  pf(ji,jj,jk) - pf(ji,jj,jk+1)  -  ptn(ji,jj,jk) * z_hdivn  )
219                  END DO
220               END DO
221            END DO
222            !
223         END SELECT
224         !
225      ELSE                ! second call : just compute the trend   (TVD scheme)
226         SELECT CASE( ktrd ) 
227         CASE( jpt_trd_xad )      ! i-advective trend
228            DO jk = 1, jpkm1     
229               DO jj = 2, jpjm1
230                  DO ji = fs_2, fs_jpim1   ! vector opt.
231                     zbtr    = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
232                     ztrdt(ji,jj,jk) = - zbtr * (  pf(ji,jj,jk) - pf(ji-1,jj,jk)  )
233                  END DO
234               END DO
235            END DO
236            !
237         CASE( jpt_trd_yad )      ! j-advective trend
238            DO jk = 1, jpkm1
239               DO jj = 2, jpjm1
240                  DO ji = fs_2, fs_jpim1   ! vector opt.
241                     zbtr    = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
242                     ztrdt(ji,jj,jk) = - zbtr * (  pf(ji,jj,jk) - pf(ji,jj-1,jk)  )
243                  END DO
244               END DO
245            END DO
246            !
247         CASE( jpt_trd_zad )      ! z-advective trend
248            DO jk = 1, jpkm1
249               DO jj = 2, jpjm1
250                  DO ji = fs_2, fs_jpim1   ! vector opt.
251                     zbtr    = 1.e0 / fse3t(ji,jj,jk)
252                     ztrdt(ji,jj,jk) = - zbtr * (  pf(ji,jj,jk) - pf(ji,jj,jk+1)  )
253                  END DO
254               END DO
255            END DO
256            !
257         END SELECT
258         !
259      ENDIF
260      !
261      CALL trd_tra( kt, ktra, ktrd, ctype, ptrd3d=ztrdt)       ! trend diagnostics
262      !
263   END SUBROUTINE trd_tra_adv
264
265
266   SUBROUTINE trd_mod( ptrdx, ptrdy, ktrd, ctype, kt, cnbpas )
267      !!---------------------------------------------------------------------
268      !!                  ***  ROUTINE trd_mod  ***
269      !!
270      !! ** Purpose : Dispatch all trends computation, e.g. vorticity, mld or
271      !!              integral constraints
272      !!----------------------------------------------------------------------
273      INTEGER         , INTENT(in   )                                   ::   kt      ! time step
274      INTEGER         , INTENT(in   )                                   ::   ktrd    ! tracer trend index
275      CHARACTER(len=3), INTENT(in   )                                   ::   ctype   ! momentum trends type ='DYN'
276      CHARACTER(len=3), INTENT(in   )                        , OPTIONAL ::   cnbpas  ! number of passage
277      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk)           ::   ptrdx   ! Temperature or U trend
278      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk)           ::   ptrdy   ! Salinity    or V trend
279      !!
280      INTEGER ::   ji, ikbu, ikbum1
281      INTEGER ::   jj, ikbv, ikbvm1
282      CHARACTER(len=3) ::   ccpas                                  ! number of passage
283      REAL(wp) ::   zua, zva                                       ! scalars
284      REAL(wp), DIMENSION(jpi,jpj) ::   ztswu, ztswv               ! 2D workspace
285      REAL(wp), DIMENSION(jpi,jpj) ::   ztbfu, ztbfv               ! 2D workspace
286      REAL(wp), DIMENSION(jpi,jpj) ::   z2dx, z2dy                 ! workspace arrays
287      !!----------------------------------------------------------------------
288
289      z2dx(:,:)   = 0.e0   ;   z2dy(:,:)   = 0.e0                  ! initialization of workspace arrays
290
291      ccpas = 'fst'                            ! Control of optional arguments
292      IF( PRESENT(cnbpas) )  ccpas = cnbpas
293
294      IF( neuler == 0 .AND. kt == nit000    ) THEN   ;   r2dt =      rdt      ! = rdtra (start with Euler time stepping)
295      ELSEIF(               kt <= nit000 + 1) THEN   ;   r2dt = 2. * rdt      ! = 2 rdttra (leapfrog)
296      ENDIF
297
298      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
299      ! I. Integral Constraints Properties for momentum and/or tracers trends
300      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
301
302      IF( ( mod(kt,ntrd) == 0 .OR. kt == nit000 .OR. kt == nitend) )   THEN
303         !
304         IF( lk_trddyn .AND. ctype == 'DYN' )   THEN       ! momentum trends
305            !
306            SELECT CASE ( ktrd )
307            CASE ( jpdyn_trd_hpg )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_hpg, ctype )   ! hydrost. pressure grad
308            CASE ( jpdyn_trd_keg )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_keg, ctype )   ! KE gradient
309            CASE ( jpdyn_trd_rvo )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_rvo, ctype )   ! relative vorticity
310            CASE ( jpdyn_trd_pvo )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_pvo, ctype )   ! planetary vorticity
311            CASE ( jpdyn_trd_ldf )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_ldf, ctype )   ! lateral diffusion
312            CASE ( jpdyn_trd_zad )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_zad, ctype )   ! vertical advection
313            CASE ( jpdyn_trd_spg )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_spg, ctype )   ! surface pressure grad.
314            CASE ( jpdyn_trd_dat )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_dat, ctype )   ! damping term
315            CASE ( jpdyn_trd_zdf )                                                         ! vertical diffusion
316               ! subtract surface forcing/bottom friction trends
317               ! from vertical diffusive momentum trends
318               ztswu(:,:) = 0.e0   ;   ztswv(:,:) = 0.e0
319               ztbfu(:,:) = 0.e0   ;   ztbfv(:,:) = 0.e0 
320               DO jj = 2, jpjm1   
321                  DO ji = fs_2, fs_jpim1   ! vector opt.
322                     ! save the surface forcing momentum fluxes
323                     ztswu(ji,jj) = taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 )
324                     ztswv(ji,jj) = tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 )
325                     ! save bottom friction momentum fluxes
326                     ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
327                     ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) )
328                     ikbum1 = MAX( ikbu-1, 1 )
329                     ikbvm1 = MAX( ikbv-1, 1 )
330                     zua = ua(ji,jj,ikbum1) * r2dt + ub(ji,jj,ikbum1)
331                     zva = va(ji,jj,ikbvm1) * r2dt + vb(ji,jj,ikbvm1)
332                     ztbfu(ji,jj) = - avmu(ji,jj,ikbu) * zua / ( fse3u(ji,jj,ikbum1)*fse3uw(ji,jj,ikbu) )
333                     ztbfv(ji,jj) = - avmv(ji,jj,ikbv) * zva / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) )
334                     !
335                     ptrdx(ji,jj,1     ) = ptrdx(ji,jj,1     ) - ztswu(ji,jj)
336                     ptrdy(ji,jj,1     ) = ptrdy(ji,jj,1     ) - ztswv(ji,jj)
337                     ptrdx(ji,jj,ikbum1) = ptrdx(ji,jj,ikbum1) - ztbfu(ji,jj)
338                     ptrdy(ji,jj,ikbvm1) = ptrdy(ji,jj,ikbvm1) - ztbfv(ji,jj)
339                  END DO
340               END DO
341               !
342               CALL trd_icp( ptrdx, ptrdy, jpicpd_zdf, ctype )   
343               CALL trd_icp( ztswu, ztswv, jpicpd_swf, ctype )                               ! wind stress forcing term
344               CALL trd_icp( ztbfu, ztbfv, jpicpd_bfr, ctype )                               ! bottom friction term
345            END SELECT
346            !
347         END IF
348         !
349      END IF
350
351      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
352      ! II. Vorticity trends
353      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
354
355      IF( lk_trdvor .AND. ctype == 'DYN' )   THEN
356         !
357         SELECT CASE ( ktrd ) 
358         CASE ( jpdyn_trd_hpg )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_prg )   ! Hydrostatique Pressure Gradient
359         CASE ( jpdyn_trd_keg )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_keg )   ! KE Gradient
360         CASE ( jpdyn_trd_rvo )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_rvo )   ! Relative Vorticity
361         CASE ( jpdyn_trd_pvo )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_pvo )   ! Planetary Vorticity Term
362         CASE ( jpdyn_trd_ldf )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_ldf )   ! Horizontal Diffusion
363         CASE ( jpdyn_trd_zad )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_zad )   ! Vertical Advection
364         CASE ( jpdyn_trd_spg )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_spg )   ! Surface Pressure Grad.
365         CASE ( jpdyn_trd_dat )   ;   CALL trd_vor_zint( ptrdx, ptrdy, jpvor_bev )   ! Beta V 
366         CASE ( jpdyn_trd_zdf )                                                      ! Vertical Diffusion
367            ! subtract surface forcing/bottom friction trends
368            ! from vertical diffusive momentum trends
369            ztswu(:,:) = 0.e0   ;   ztswv(:,:) = 0.e0
370            ztbfu(:,:) = 0.e0   ;   ztbfv(:,:) = 0.e0 
371            DO jj = 2, jpjm1   
372               DO ji = fs_2, fs_jpim1   ! vector opt.
373                  ! save the surface forcing momentum fluxes
374                  ztswu(ji,jj) = taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 )
375                  ztswv(ji,jj) = tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 )
376                  ! save bottom friction momentum fluxes
377                  ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) )
378                  ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) )
379                  ikbum1 = MAX( ikbu-1, 1 )
380                  ikbvm1 = MAX( ikbv-1, 1 )
381                  zua = ua(ji,jj,ikbum1) * r2dt + ub(ji,jj,ikbum1)
382                  zva = va(ji,jj,ikbvm1) * r2dt + vb(ji,jj,ikbvm1)
383                  ztbfu(ji,jj) = - avmu(ji,jj,ikbu) * zua / ( fse3u(ji,jj,ikbum1)*fse3uw(ji,jj,ikbu) )
384                  ztbfv(ji,jj) = - avmv(ji,jj,ikbv) * zva / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) )
385                  !
386                  ptrdx(ji,jj,1     ) = ptrdx(ji,jj,1     ) - ztswu(ji,jj)
387                  ptrdx(ji,jj,ikbum1) = ptrdx(ji,jj,ikbum1) - ztbfu(ji,jj)
388                  ptrdy(ji,jj,1     ) = ptrdy(ji,jj,1     ) - ztswv(ji,jj)
389                  ptrdy(ji,jj,ikbvm1) = ptrdy(ji,jj,ikbvm1) - ztbfv(ji,jj)
390               END DO
391            END DO
392            !
393            CALL trd_vor_zint( ptrdx, ptrdy, jpvor_zdf )   
394            CALL trd_vor_zint( ztswu, ztswv, jpvor_swf )                               ! Wind stress forcing term
395            CALL trd_vor_zint( ztbfu, ztbfv, jpvor_bfr )                               ! Bottom friction term
396         END SELECT
397         !
398      ENDIF
399      !
400   END SUBROUTINE trd_mod
401
402#   else
403   !!----------------------------------------------------------------------
404   !!   Default case :                                         Empty module
405   !!----------------------------------------------------------------------
406   USE trdmod_oce      ! ocean variables trends
407   USE trdvor          ! ocean vorticity trends
408   USE trdicp          ! ocean bassin integral constraints properties
409   USE trdmld          ! ocean active mixed layer tracers trends
410
411CONTAINS
412
413   SUBROUTINE trd_tra( kt, ktra, ktrd, ctype, ptrd2d, ptrd3d, cnbpas )
414      INTEGER ::   kt, ktra, ktrd
415      CHARACTER(len=3) ::   ctype   ! tracers type 'TRA' or 'TRC'
416      REAL, DIMENSION(:,:)  , OPTIONAL ::   ptrd2d  ! Temperature or U trend
417      REAL, DIMENSION(:,:,:), OPTIONAL ::   ptrd3d  ! Temperature or U trend
418      CHARACTER(len=3)      , OPTIONAL ::   cnbpas  ! number of passage
419      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', kt, ktra, ktrd, ctype,   &
420         &                                                               ptrd2d(1,1), ptrd3d(1,1,1), cnbpas
421   END SUBROUTINE trd_tra
422
423   SUBROUTINE trd_tra_adv( kt, ktra, ktrd, ctype, pf, pun, ptn, cnbpas )
424      INTEGER ::   kt, ktra, ktrd 
425      CHARACTER(len=3) ::   ctype
426      REAL, DIMENSION(:,:,:) ::   pf, pun, ptn
427      CHARACTER(len=3), OPTIONAL ::   cnbpas
428      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', kt, ktra, ktrd, ctype, pf(1,1,1),   &
429         &                                                               pun(1,1,1), ptn(1,1,1), cnbpas
430   END SUBROUTINE trd_tra_adv
431
432   SUBROUTINE trd_mod(ptrd3dx, ptrd3dy, ktrd , ctype, kt, cnbpas)   ! Empty routine
433      REAL, DIMENSION(:,:,:) ::   ptrd3dx, ptrd3dy   
434      INTEGER ::   ktrd, kt
435      CHARACTER(len=3), INTENT( in ) ::  ctype                     ! momentum or tracers trends type
436      CHARACTER(len=3), INTENT( in ), OPTIONAL ::   cnbpas         ! number of passage
437      WRITE(*,*) 'trd_3d: You should not have seen this print! error ?', ptrd3dx(1,1,1)
438      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ptrd3dy(1,1,1)
439      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ktrd
440      WRITE(*,*) ' "   ": You should not have seen this print! error ?', ctype
441      WRITE(*,*) ' "   ": You should not have seen this print! error ?', kt
442      WRITE(*,*) ' "   ": You should not have seen this print! error ?', cnbpas
443   END SUBROUTINE trd_mod
444#   endif
445
446   SUBROUTINE trd_mod_init
447      !!----------------------------------------------------------------------
448      !!                  ***  ROUTINE trd_mod_init  ***
449      !!
450      !! ** Purpose :   Initialization of activated trends
451      !!----------------------------------------------------------------------
452      USE in_out_manager          ! I/O manager
453
454      NAMELIST/namtrd/ ntrd, nctls, ln_trdmld_restart, ucf, ln_trdmld_instant
455      !!----------------------------------------------------------------------
456
457      IF( l_trdtra .OR. l_trddyn )   THEN
458         REWIND( numnam )
459         READ  ( numnam, namtrd )      ! namelist namtrd : trends diagnostic
460
461         IF(lwp) THEN
462            WRITE(numout,*)
463            WRITE(numout,*) ' trd_mod_init : Momentum/Tracers trends'
464            WRITE(numout,*) ' ~~~~~~~~~~~~~'
465            WRITE(numout,*) '       Namelist namtrd : set trends parameters'
466            WRITE(numout,*) '           * frequency of trends diagnostics   ntrd               = ', ntrd
467            WRITE(numout,*) '           * control surface type              nctls              = ', nctls
468            WRITE(numout,*) '           * restart for ML diagnostics        ln_trdmld_restart  = ', ln_trdmld_restart
469            WRITE(numout,*) '           * instantaneous or mean ML T/S      ln_trdmld_instant  = ', ln_trdmld_instant
470            WRITE(numout,*) '           * unit conversion factor            ucf                = ', ucf
471        ENDIF
472      ENDIF
473      !
474      IF( lk_trddyn .OR. lk_trdtra )    CALL trd_icp_init       ! integral constraints trends
475      IF( lk_trdmld                )    CALL trd_mld_init       ! mixed-layer trends (active  tracers) 
476      IF( lk_trdvor                )    CALL trd_vor_init       ! vorticity trends       
477      !
478   END SUBROUTINE trd_mod_init
479
480   !!======================================================================
481END MODULE trdmod
Note: See TracBrowser for help on using the repository browser.