source: NEMO/branches/UKMO/NEMO_4.0_momentum_trends/src/OCE/TRD/trddyn.F90 @ 11903

Last change on this file since 11903 was 11903, checked in by davestorkey, 9 months ago

UKMO/NEMO_4.0_momentum_trends : modification to bottom friction trend diagnostic.

File size: 17.6 KB
Line 
1MODULE trddyn
2   !!======================================================================
3   !!                       ***  MODULE  trddyn  ***
4   !! Ocean diagnostics:  ocean dynamic trends
5   !!=====================================================================
6   !! History :  3.5  !  2012-02  (G. Madec) creation from trdmod: split DYN and TRA trends
7   !!                                        and manage  3D trends output for U, V, and KE
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   trd_dyn       : manage the type of momentum trend diagnostics (3D I/O, domain averaged, KE)
12   !!   trd_dyn_iom   : output 3D momentum and/or tracer trends using IOM
13   !!   trd_dyn_init  : initialization step
14   !!----------------------------------------------------------------------
15   USE oce            ! ocean dynamics and tracers variables
16   USE dom_oce        ! ocean space and time domain variables
17   USE phycst         ! physical constants
18   USE sbc_oce        ! surface boundary condition: ocean
19   USE zdf_oce        ! ocean vertical physics: variables
20   USE trd_oce        ! trends: ocean variables
21   USE trdken         ! trends: Kinetic ENergy
22   USE trdglo         ! trends: global domain averaged
23   USE trdvor         ! trends: vertical averaged vorticity
24   USE trdmxl         ! trends: mixed layer averaged
25   !
26   USE in_out_manager ! I/O manager
27   USE lbclnk         ! lateral boundary condition
28   USE iom            ! I/O manager library
29   USE lib_mpp        ! MPP library
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC trd_dyn        ! called by all dynXXX modules
35
36   INTERFACE trd_dyn
37      module procedure trd_dyn_3d, trd_dyn_2d
38   END INTERFACE
39
40   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zutrd_hpg, zvtrd_hpg
41   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zutrd_pvo, zvtrd_pvo
42   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: zutrd_bfr, zvtrd_bfr
43
44   !! * Substitutions
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
48   !! $Id$
49   !! Software governed by the CeCILL license (see ./LICENSE)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   SUBROUTINE trd_dyn_3d( putrd, pvtrd, ktrd, kt )
54      !!---------------------------------------------------------------------
55      !!                  ***  ROUTINE trd_dyn_3d  ***
56      !!
57      !! ** Purpose :   Dispatch momentum trend computation, e.g. 3D output,
58      !!              integral constraints, barotropic vorticity, kinetic enrgy,
59      !!              and/or mixed layer budget.
60      !!----------------------------------------------------------------------
61      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends
62      INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
63      INTEGER                   , INTENT(in   ) ::   kt             ! time step
64      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zue, zve       ! temporary 2D arrays
65      INTEGER                                   ::   jk
66      !!----------------------------------------------------------------------
67      !
68      putrd(:,:,:) = putrd(:,:,:) * umask(:,:,:)                       ! mask the trends
69      pvtrd(:,:,:) = pvtrd(:,:,:) * vmask(:,:,:)
70      !
71
72!!gm NB : here a lbc_lnk should probably be added
73
74      SELECT CASE( ktrd )
75      CASE( jpdyn_hpg_save ) 
76         !
77         ! save 3D HPG trends to possibly have barotropic part corrected later before writing out
78         ALLOCATE( zutrd_hpg(jpi,jpj,jpk), zvtrd_hpg(jpi,jpj,jpk) )
79         zutrd_hpg(:,:,:) = putrd(:,:,:)
80         zvtrd_hpg(:,:,:) = pvtrd(:,:,:)
81
82      CASE( jpdyn_pvo_save ) 
83         !
84         ! save 3D coriolis trends to possibly have barotropic part corrected later before writing out
85         ALLOCATE( zutrd_pvo(jpi,jpj,jpk), zvtrd_pvo(jpi,jpj,jpk) )
86         zutrd_pvo(:,:,:) = putrd(:,:,:)
87         zvtrd_pvo(:,:,:) = pvtrd(:,:,:)
88
89      CASE DEFAULT
90
91         IF( ktrd == jpdyn_spg ) THEN
92            ! For explicit scheme SPG trends come here as 3D fields
93            ! Add SPG trend to 3D HPG trend and also output as 2D diagnostic in own right.
94            CALL trd_dyn_iom_2d( putrd(:,:,1), pvtrd(:,:,1), jpdyn_spg, kt ) 
95            putrd(:,:,:) = putrd(:,:,:) + zutrd_hpg(:,:,:) 
96            pvtrd(:,:,:) = pvtrd(:,:,:) + zvtrd_hpg(:,:,:) 
97            DEALLOCATE( zutrd_hpg, zvtrd_hpg )
98         ENDIF
99
100         IF( ktrd == jpdyn_bfr ) THEN
101            !
102            ! Add 3D part of BFR trend minus its depth-mean part to depth-mean trend already saved.
103            ALLOCATE( zue(jpi,jpj), zve(jpi,jpj) )
104            zue(:,:) = e3u_a(:,:,1) * putrd(:,:,1) * umask(:,:,1)
105            zve(:,:) = e3v_a(:,:,1) * pvtrd(:,:,1) * vmask(:,:,1)
106            DO jk = 2, jpkm1
107               zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * putrd(:,:,jk) * umask(:,:,jk)
108               zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * pvtrd(:,:,jk) * vmask(:,:,jk)
109            END DO
110            DO jk = 1, jpkm1
111               putrd(:,:,jk) = zutrd_bfr(:,:,jk) + putrd(:,:,jk) - zue(:,:) * r1_hu_a(:,:)
112               pvtrd(:,:,jk) = zvtrd_bfr(:,:,jk) + pvtrd(:,:,jk) - zve(:,:) * r1_hv_a(:,:)
113            END DO
114            DEALLOCATE( zue, zve, zutrd_bfr, zvtrd_bfr )
115         ENDIF
116
117         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
118         !   3D output of momentum and/or tracers trends using IOM interface
119         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
120         IF( ln_dyn_trd )   CALL trd_dyn_iom_3d( putrd, pvtrd, ktrd, kt )
121
122         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
123         !  Integral Constraints Properties for momentum and/or tracers trends
124         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
125         IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt )
126
127         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
128         !  Kinetic Energy trends
129         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
130         IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt )
131
132         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
133         !  Vorticity trends
134         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
135         IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt )
136
137         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
138         !  Mixed layer trends for active tracers
139         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
140         !!gm      IF( ln_dyn_mxl )   CALL trd_mxl_dyn   
141         !
142      END SELECT
143      !
144   END SUBROUTINE trd_dyn_3d
145
146
147   SUBROUTINE trd_dyn_2d( putrd, pvtrd, ktrd, kt )
148      !!---------------------------------------------------------------------
149      !!                  ***  ROUTINE trd_mod  ***
150      !!
151      !! ** Purpose :   Dispatch momentum trend computation, e.g. 2D output,
152      !!              integral constraints, barotropic vorticity, kinetic enrgy,
153      !!              and/or mixed layer budget.
154      !!----------------------------------------------------------------------
155      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends
156      INTEGER                 , INTENT(in   ) ::   ktrd           ! trend index
157      INTEGER                 , INTENT(in   ) ::   kt             ! time step
158      INTEGER                                 ::   jk
159      !!----------------------------------------------------------------------
160      !
161      putrd(:,:) = putrd(:,:) * umask(:,:,1)                       ! mask the trends
162      pvtrd(:,:) = pvtrd(:,:) * vmask(:,:,1)
163      !
164
165!!gm NB : here a lbc_lnk should probably be added
166
167      SELECT CASE(ktrd)
168
169      CASE ( jpdyn_hpg_corr )
170         !
171         ! Remove "first-guess" SPG trend from 3D HPG trend.
172         DO jk = 1, jpkm1
173            zutrd_hpg(:,:,jk) = zutrd_hpg(:,:,jk) - putrd(:,:)
174            zvtrd_hpg(:,:,jk) = zvtrd_hpg(:,:,jk) - pvtrd(:,:)
175         ENDDO
176         CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_hpg_corr, kt ) 
177
178      CASE( jpdyn_pvo_corr )
179         !
180         ! Remove "first-guess" barotropic coriolis trend from 3D PVO trend.
181         DO jk = 1, jpkm1
182            zutrd_pvo(:,:,jk) = zutrd_pvo(:,:,jk) - putrd(:,:)
183            zvtrd_pvo(:,:,jk) = zvtrd_pvo(:,:,jk) - pvtrd(:,:)
184         ENDDO
185         CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_pvo_corr, kt ) 
186
187      CASE( jpdyn_spg )
188          !
189          ! For split-explicit scheme SPG trends come here as 2D fields
190          ! Add SPG trend to 3D HPG trend and also output as 2D diagnostic in own right.
191          DO jk = 1, jpkm1
192             zutrd_hpg(:,:,jk) = zutrd_hpg(:,:,jk) + putrd(:,:)
193             zvtrd_hpg(:,:,jk) = zvtrd_hpg(:,:,jk) + pvtrd(:,:)
194          ENDDO
195          CALL trd_dyn_iom_2d( putrd(:,:), pvtrd(:,:), jpdyn_spg, kt ) 
196          CALL trd_dyn_3d( zutrd_hpg, zvtrd_hpg, jpdyn_hpg, kt )
197          DEALLOCATE( zutrd_hpg, zvtrd_hpg )
198
199      CASE( jpdyn_pvo )
200          !
201          ! Add 2D PVO trend to 3D PVO trend and also output as diagnostic in own right.
202          DO jk = 1, jpkm1
203             zutrd_pvo(:,:,jk) = zutrd_pvo(:,:,jk) + putrd(:,:)
204             zvtrd_pvo(:,:,jk) = zvtrd_pvo(:,:,jk) + pvtrd(:,:)
205          ENDDO
206          CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_pvo, kt ) 
207          CALL trd_dyn_3d( zutrd_pvo, zvtrd_pvo, jpdyn_pvo, kt )
208          DEALLOCATE( zutrd_pvo, zvtrd_pvo )
209
210      CASE( jpdyn_bfr )
211          !
212          ! Create 3D BFR trend from 2D field and also output 2D field as diagnostic in own right.
213          ALLOCATE( zutrd_bfr(jpi,jpj,jpk), zvtrd_bfr(jpi,jpj,jpk) )
214          zutrd_bfr(:,:,:) = 0.0
215          zvtrd_bfr(:,:,:) = 0.0
216          DO jk = 1, jpkm1
217             zutrd_bfr(:,:,jk) = putrd(:,:)
218             zvtrd_bfr(:,:,jk) = pvtrd(:,:)
219          ENDDO
220          CALL trd_dyn_iom_2d( putrd, pvtrd, jpdyn_bfr, kt ) 
221
222      CASE DEFAULT 
223
224         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
225         !   2D output of momentum and/or tracers trends using IOM interface
226         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
227         IF( ln_dyn_trd )   CALL trd_dyn_iom_2d( putrd, pvtrd, ktrd, kt )
228         
229      END SELECT
230
231!!$   CALLS TO THESE ROUTINES FOR 2D DIAGOSTICS NOT CODED YET
232!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
233!!$      !  Integral Constraints Properties for momentum and/or tracers trends
234!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
235!!$      IF( ln_glo_trd )   CALL trd_glo( putrd, pvtrd, ktrd, 'DYN', kt )
236!!$
237!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
238!!$      !  Kinetic Energy trends
239!!$      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
240!!$      IF( ln_KE_trd  )   CALL trd_ken( putrd, pvtrd, ktrd, kt )
241!!$
242!!$      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
243!!$      !  Vorticity trends
244!!$      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
245!!$      IF( ln_vor_trd )   CALL trd_vor( putrd, pvtrd, ktrd, kt )
246
247      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
248      !  Mixed layer trends for active tracers
249      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
250!!gm      IF( ln_dyn_mxl )   CALL trd_mxl_dyn   
251      !
252   END SUBROUTINE trd_dyn_2d
253
254
255   SUBROUTINE trd_dyn_iom_3d( putrd, pvtrd, ktrd, kt )
256      !!---------------------------------------------------------------------
257      !!                  ***  ROUTINE trd_dyn_iom  ***
258      !!
259      !! ** Purpose :   output 3D trends using IOM
260      !!----------------------------------------------------------------------
261      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends
262      INTEGER                   , INTENT(in   ) ::   ktrd           ! trend index
263      INTEGER                   , INTENT(in   ) ::   kt             ! time step
264      !
265      INTEGER ::   ji, jj, jk   ! dummy loop indices
266      INTEGER ::   ikbu, ikbv   ! local integers
267      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace
268      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z3dx, z3dy   ! 3D workspace
269      !!----------------------------------------------------------------------
270      !
271      SELECT CASE( ktrd )
272      CASE( jpdyn_hpg )   ;   CALL iom_put( "utrd_hpg", putrd )    ! hydrostatic pressure gradient
273                              CALL iom_put( "vtrd_hpg", pvtrd )
274      CASE( jpdyn_pvo )   ;   CALL iom_put( "utrd_pvo", putrd )    ! planetary vorticity
275                              CALL iom_put( "vtrd_pvo", pvtrd )
276      CASE( jpdyn_rvo )   ;   CALL iom_put( "utrd_rvo", putrd )    ! relative  vorticity     (or metric term)
277                              CALL iom_put( "vtrd_rvo", pvtrd )
278      CASE( jpdyn_keg )   ;   CALL iom_put( "utrd_keg", putrd )    ! Kinetic Energy gradient (or had)
279                              CALL iom_put( "vtrd_keg", pvtrd )
280                              ALLOCATE( z3dx(jpi,jpj,jpk) , z3dy(jpi,jpj,jpk) )
281                              z3dx(:,:,:) = 0._wp                  ! U.dxU & V.dyV (approximation)
282                              z3dy(:,:,:) = 0._wp
283                              DO jk = 1, jpkm1   ! no mask as un,vn are masked
284                                 DO jj = 2, jpjm1
285                                    DO ji = 2, jpim1
286                                       z3dx(ji,jj,jk) = un(ji,jj,jk) * ( un(ji+1,jj,jk) - un(ji-1,jj,jk) ) / ( 2._wp * e1u(ji,jj) )
287                                       z3dy(ji,jj,jk) = vn(ji,jj,jk) * ( vn(ji,jj+1,jk) - vn(ji,jj-1,jk) ) / ( 2._wp * e2v(ji,jj) )
288                                    END DO
289                                 END DO
290                              END DO
291                              CALL lbc_lnk_multi( 'trddyn', z3dx, 'U', -1., z3dy, 'V', -1. )
292                              CALL iom_put( "utrd_udx", z3dx  )
293                              CALL iom_put( "vtrd_vdy", z3dy  )
294                              DEALLOCATE( z3dx , z3dy )
295      CASE( jpdyn_zad )   ;   CALL iom_put( "utrd_zad", putrd )    ! vertical advection
296                              CALL iom_put( "vtrd_zad", pvtrd )
297      CASE( jpdyn_ldf )   ;   CALL iom_put( "utrd_ldf", putrd )    ! lateral  diffusion
298                              CALL iom_put( "vtrd_ldf", pvtrd )
299      CASE( jpdyn_zdf )   ;   CALL iom_put( "utrd_zdf", putrd )    ! vertical diffusion
300                              CALL iom_put( "vtrd_zdf", pvtrd )
301                              !
302                              !                                    ! wind stress trends
303                              ALLOCATE( z2dx(jpi,jpj) , z2dy(jpi,jpj) )
304                              z2dx(:,:) = ( utau_b(:,:) + utau(:,:) ) / ( e3u_n(:,:,1) * rau0 )
305                              z2dy(:,:) = ( vtau_b(:,:) + vtau(:,:) ) / ( e3v_n(:,:,1) * rau0 )
306                              CALL iom_put( "utrd_tau", z2dx )
307                              CALL iom_put( "vtrd_tau", z2dy )
308                              DEALLOCATE( z2dx , z2dy )
309      CASE( jpdyn_bfr )   ;   CALL iom_put( "utrd_bfr", putrd )    ! bottom friction (explicit case)
310                              CALL iom_put( "vtrd_bfr", pvtrd )
311      CASE( jpdyn_bfri)   ;   CALL iom_put( "utrd_bfri", putrd )    ! bottom friction (implicit case)
312                              CALL iom_put( "vtrd_bfri", pvtrd )
313      CASE( jpdyn_atf )   ;   CALL iom_put( "utrd_atf", putrd )        ! asselin filter trends
314                              CALL iom_put( "vtrd_atf", pvtrd )
315      END SELECT
316      !
317   END SUBROUTINE trd_dyn_iom_3d
318
319
320   SUBROUTINE trd_dyn_iom_2d( putrd, pvtrd, ktrd, kt )
321      !!---------------------------------------------------------------------
322      !!                  ***  ROUTINE trd_dyn_iom  ***
323      !!
324      !! ** Purpose :   output 2D trends using IOM
325      !!----------------------------------------------------------------------
326      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   putrd, pvtrd   ! U and V trends
327      INTEGER                 , INTENT(in   ) ::   ktrd           ! trend index
328      INTEGER                 , INTENT(in   ) ::   kt             ! time step
329      !
330      INTEGER ::   ji, jj, jk   ! dummy loop indices
331      INTEGER ::   ikbu, ikbv   ! local integers
332      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   z2dx, z2dy   ! 2D workspace
333      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z3dx, z3dy   ! 3D workspace
334      !!----------------------------------------------------------------------
335      !
336      SELECT CASE( ktrd )
337      CASE( jpdyn_spg )      ;   CALL iom_put( "utrd_spg2d", putrd )      ! surface pressure gradient
338                                 CALL iom_put( "vtrd_spg2d", pvtrd )
339      CASE( jpdyn_pvo )      ;   CALL iom_put( "utrd_pvo2d", putrd )      ! planetary vorticity (barotropic part)
340                                 CALL iom_put( "vtrd_pvo2d", pvtrd )
341      CASE( jpdyn_hpg_corr ) ;   CALL iom_put( "utrd_hpg_corr", putrd )   ! horizontal pressure gradient correction
342                                 CALL iom_put( "vtrd_hpg_corr", pvtrd )
343      CASE( jpdyn_pvo_corr ) ;   CALL iom_put( "utrd_pvo_corr", putrd )   ! planetary vorticity correction
344                                 CALL iom_put( "vtrd_pvo_corr", pvtrd )
345      CASE( jpdyn_bfr )      ;   CALL iom_put( "utrd_bfr2d", putrd )      ! bottom friction due to barotropic currents
346                                 CALL iom_put( "vtrd_bfr2d", pvtrd )
347      END SELECT
348      !
349   END SUBROUTINE trd_dyn_iom_2d
350
351   !!======================================================================
352END MODULE trddyn
Note: See TracBrowser for help on using the repository browser.