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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRD – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRD/trdmod.F90 @ 4443

Last change on this file since 4443 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

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