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.
dynatf.F90 in NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/DYN/dynatf.F90 @ 13710

Last change on this file since 13710 was 13710, checked in by emanuelaclementi, 4 years ago

branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves: merge with trunk@13708, see #2155 and #2339

  • Property svn:keywords set to Id
File size: 19.0 KB
RevLine 
[11050]1MODULE dynatf
[1502]2   !!=========================================================================
[11050]3   !!                       ***  MODULE  dynatf  ***
4   !! Ocean dynamics: time filtering
[1502]5   !!=========================================================================
[1438]6   !! History :  OPA  !  1987-02  (P. Andrich, D. L Hostis)  Original code
7   !!                 !  1990-10  (C. Levy, G. Madec)
8   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
9   !!            8.0  !  1997-02  (G. Madec & M. Imbard)  opa, release 8.0
10   !!            8.2  !  1997-04  (A. Weaver)  Euler forward step
11   !!             -   !  1997-06  (G. Madec)  lateral boudary cond., lbc routine
12   !!    NEMO    1.0  !  2002-08  (G. Madec)  F90: Free form and module
13   !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond.
14   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
15   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.
[1502]16   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option
[2528]17   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module
[2723]18   !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
[4292]19   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
[6140]20   !!            3.6  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
[5930]21   !!            3.7  !  2015-11  (J. Chanut) Free surface simplification
[11475]22   !!            4.1  !  2019-08  (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatf.F90. Now just does time filtering.
[1502]23   !!-------------------------------------------------------------------------
[1438]24 
[11050]25   !!----------------------------------------------------------------------------------------------
26   !!   dyn_atf       : apply Asselin time filtering to "now" velocities and vertical scale factors
27   !!----------------------------------------------------------------------------------------------
[6140]28   USE oce            ! ocean dynamics and tracers
29   USE dom_oce        ! ocean space and time domain
30   USE sbc_oce        ! Surface boundary condition: ocean fields
[9023]31   USE sbcrnf         ! river runoffs
[6140]32   USE phycst         ! physical constants
33   USE dynadv         ! dynamics: vector invariant versus flux form
34   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
35   USE domvvl         ! variable volume
[13710]36   USE bdy_oce , ONLY : ln_bdy
[6140]37   USE bdydta         ! ocean open boundary conditions
38   USE bdydyn         ! ocean open boundary conditions
39   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
40   USE trd_oce        ! trends: ocean variables
41   USE trddyn         ! trend manager: dynamics
42   USE trdken         ! trend manager: kinetic energy
[12150]43   USE isf_oce   , ONLY: ln_isf     ! ice shelf
44   USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine
[4990]45   !
[6140]46   USE in_out_manager ! I/O manager
47   USE iom            ! I/O manager library
48   USE lbclnk         ! lateral boundary condition (or mpp link)
49   USE lib_mpp        ! MPP library
50   USE prtctl         ! Print control
51   USE timing         ! Timing
[13710]52   USE zdfdrg ,  ONLY : ln_drgice_imp, rCdU_top
[2528]53#if defined key_agrif
[9570]54   USE agrif_oce_interp
[11050]55#endif
[3]56
57   IMPLICIT NONE
58   PRIVATE
59
[11050]60   PUBLIC    dyn_atf   ! routine called by step.F90
[1438]61
[13710]62#if defined key_qco
63   !!----------------------------------------------------------------------
64   !!   'key_qco'      EMPTY ROUTINE     Quasi-Eulerian vertical coordonate
65   !!----------------------------------------------------------------------
66CONTAINS
67
68   SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v )
69      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index
70      INTEGER                             , INTENT(in   ) :: Kbb, Kmm, Kaa    ! before and after time level indices
71      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities to be time filtered
72      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: pe3t, pe3u, pe3v ! scale factors to be time filtered
73
74      WRITE(*,*) 'dyn_atf: You should not have seen this print! error?', kt
75   END SUBROUTINE dyn_atf
76
77#else
78
[12340]79   !! * Substitutions
80#  include "do_loop_substitute.h90"
[2715]81   !!----------------------------------------------------------------------
[9598]82   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1438]83   !! $Id$
[10068]84   !! Software governed by the CeCILL license (see ./LICENSE)
[2715]85   !!----------------------------------------------------------------------
[3]86CONTAINS
87
[11050]88   SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v )
[3]89      !!----------------------------------------------------------------------
[11050]90      !!                  ***  ROUTINE dyn_atf  ***
[3]91      !!                   
[5930]92      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
[11475]93      !!             condition on the after velocity and apply the Asselin time
94      !!             filter to the now fields.
[3]95      !!
[5930]96      !! ** Method  : * Ensure after velocities transport matches time splitting
97      !!             estimate (ln_dynspg_ts=T)
[3]98      !!
[1502]99      !!              * Apply lateral boundary conditions on after velocity
100      !!             at the local domain boundaries through lbc_lnk call,
[7646]101      !!             at the one-way open boundaries (ln_bdy=T),
[4990]102      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
[3]103      !!
[11475]104      !!              * Apply the Asselin time filter to the now fields
[1502]105      !!             arrays to start the next time step:
[11475]106      !!                (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm))
[12489]107      !!                                    + rn_atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ]
[6140]108      !!             Note that with flux form advection and non linear free surface,
109      !!             the time filter is applied on thickness weighted velocity.
[11475]110      !!             As a result, dyn_atf MUST be called after tra_atf.
[1502]111      !!
[11475]112      !! ** Action :   puu(Kmm),pvv(Kmm)   filtered now horizontal velocity
[3]113      !!----------------------------------------------------------------------
[11050]114      INTEGER                             , INTENT(in   ) :: kt               ! ocean time-step index
115      INTEGER                             , INTENT(in   ) :: Kbb, Kmm, Kaa    ! before and after time level indices
116      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv         ! velocities to be time filtered
117      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: pe3t, pe3u, pe3v ! scale factors to be time filtered
[2715]118      !
[3]119      INTEGER  ::   ji, jj, jk   ! dummy loop indices
[11050]120      REAL(wp) ::   zue3a, zue3n, zue3b, zcoef    ! local scalars
121      REAL(wp) ::   zve3a, zve3n, zve3b, z1_2dt   !   -      -
[12372]122      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve, zwfld
[13710]123      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zutau, zvtau
[11050]124      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3t_f, ze3u_f, ze3v_f, zua, zva 
[1502]125      !!----------------------------------------------------------------------
[3294]126      !
[11050]127      IF( ln_timing    )   CALL timing_start('dyn_atf')
[9019]128      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
129      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
[3294]130      !
[3]131      IF( kt == nit000 ) THEN
132         IF(lwp) WRITE(numout,*)
[11050]133         IF(lwp) WRITE(numout,*) 'dyn_atf : Asselin time filtering'
[3]134         IF(lwp) WRITE(numout,*) '~~~~~~~'
135      ENDIF
136
[5930]137      IF ( ln_dynspg_ts ) THEN
138         ! Ensure below that barotropic velocities match time splitting estimate
139         ! Compute actual transport and replace it with ts estimate at "after" time step
[11050]140         zue(:,:) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
141         zve(:,:) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
[5930]142         DO jk = 2, jpkm1
[11050]143            zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
144            zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
[1502]145         END DO
146         DO jk = 1, jpkm1
[11050]147            puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
148            pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
[592]149         END DO
[6140]150         !
151         IF( .NOT.ln_bt_fw ) THEN
[5930]152            ! Remove advective velocity from "now velocities"
153            ! prior to asselin filtering     
154            ! In the forward case, this is done below after asselin filtering   
155            ! so that asselin contribution is removed at the same time
156            DO jk = 1, jpkm1
[11050]157               puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)
158               pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)
[7753]159            END DO 
[5930]160         ENDIF
[4292]161      ENDIF
162
[1502]163      ! Update after velocity on domain lateral boundaries
164      ! --------------------------------------------------     
[5930]165# if defined key_agrif
166      CALL Agrif_dyn( kt )             !* AGRIF zoom boundaries
167# endif
168      !
[13710]169      CALL lbc_lnk_multi( 'dynatf', puu(:,:,:,Kaa), 'U', -1.0_wp, pvv(:,:,:,Kaa), 'V', -1.0_wp )     !* local domain boundaries
[1502]170      !
171      !                                !* BDY open boundaries
[11050]172      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )
173      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )
[3294]174
175!!$   Do we need a call to bdy_vol here??
176      !
[4990]177      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
178         !
179         !                                  ! Kinetic energy and Conversion
[11050]180         IF( ln_KE_trd  )   CALL trd_dyn( puu(:,:,:,Kaa), pvv(:,:,:,Kaa), jpdyn_ken, kt, Kmm )
[4990]181         !
182         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
[12489]183            zua(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) * r1_Dt
184            zva(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) * r1_Dt
[4990]185            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
186            CALL iom_put( "vtrd_tot", zva )
187         ENDIF
188         !
[11050]189         zua(:,:,:) = puu(:,:,:,Kmm)             ! save the now velocity before the asselin filter
190         zva(:,:,:) = pvv(:,:,:,Kmm)             ! (caution: there will be a shift by 1 timestep in the
[7753]191         !                                  !  computation of the asselin filter trends)
[4990]192      ENDIF
193
[1438]194      ! Time filter and swap of dynamics arrays
195      ! ------------------------------------------
[9226]196         
[12489]197      IF( .NOT. l_1st_euler ) THEN    !* Leap-Frog : Asselin time filter
[2528]198         !                                ! =============!
[6140]199         IF( ln_linssh ) THEN             ! Fixed volume !
[2528]200            !                             ! =============!
[13710]201            DO_3D( 1, 1, 1, 1, 1, jpkm1 )
[12489]202               puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
203               pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
[12340]204            END_3D
[2528]205            !                             ! ================!
206         ELSE                             ! Variable volume !
207            !                             ! ================!
[11050]208            ! Time-filtered scale factor at t-points
[4292]209            ! ----------------------------------------------------
[12372]210            ALLOCATE( ze3t_f(jpi,jpj,jpk), zwfld(jpi,jpj) )
[9023]211            DO jk = 1, jpkm1
[12489]212               ze3t_f(:,:,jk) = pe3t(:,:,jk,Kmm) + rn_atfp * ( pe3t(:,:,jk,Kbb) - 2._wp * pe3t(:,:,jk,Kmm) + pe3t(:,:,jk,Kaa) )
[9023]213            END DO
214            ! Add volume filter correction: compatibility with tracer advection scheme
[12372]215            ! => time filter + conservation correction
[12489]216            zcoef = rn_atfp * rn_Dt * r1_rho0
[12372]217            zwfld(:,:) = emp_b(:,:) - emp(:,:)
218            IF ( ln_rnf ) zwfld(:,:) =  zwfld(:,:) - ( rnf_b(:,:) - rnf(:,:) )
[13710]219
[12372]220            DO jk = 1, jpkm1
221               ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) &
222                              &                        * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 
223            END DO
[2528]224            !
[12150]225            ! ice shelf melting (deal separately as it can be in depth)
226            ! PM: we could probably define a generic subroutine to do the in depth correction
227            !     to manage rnf, isf and possibly in the futur icb, tide water glacier (...)
228            !     ...(kt, coef, ktop, kbot, hz, fwf_b, fwf)
[12489]229            IF ( ln_isf ) CALL isf_dynatf( kt, Kmm, ze3t_f, rn_atfp * rn_Dt )
[12150]230            !
[11050]231            pe3t(:,:,1:jpkm1,Kmm) = ze3t_f(:,:,1:jpkm1)        ! filtered scale factor at T-points
232            !
[6140]233            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
234               ! Before filtered scale factor at (u/v)-points
[11050]235               CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3u(:,:,:,Kmm), 'U' )
236               CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' )
[13710]237               DO_3D( 1, 1, 1, 1, 1, jpkm1 )
[12489]238                  puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )
239                  pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )
[12340]240               END_3D
[2528]241               !
[6140]242            ELSE                          ! Asselin filter applied on thickness weighted velocity
243               !
[9019]244               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
[11050]245               ! Now filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
246               CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3u_f, 'U' )
247               CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), ze3v_f, 'V' )
[13710]248               DO_3D( 1, 1, 1, 1, 1, jpkm1 )
[12340]249                  zue3a = pe3u(ji,jj,jk,Kaa) * puu(ji,jj,jk,Kaa)
250                  zve3a = pe3v(ji,jj,jk,Kaa) * pvv(ji,jj,jk,Kaa)
251                  zue3n = pe3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
252                  zve3n = pe3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
253                  zue3b = pe3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb)
254                  zve3b = pe3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb)
255                  !
[12489]256                  puu(ji,jj,jk,Kmm) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
257                  pvv(ji,jj,jk,Kmm) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
[12340]258               END_3D
[11050]259               pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) 
260               pe3v(:,:,1:jpkm1,Kmm) = ze3v_f(:,:,1:jpkm1)
[6140]261               !
[9019]262               DEALLOCATE( ze3u_f , ze3v_f )
[2528]263            ENDIF
264            !
[12372]265            DEALLOCATE( ze3t_f, zwfld )
[3]266         ENDIF
[2528]267         !
[6140]268         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
[11050]269            ! Revert filtered "now" velocities to time split estimate
[4312]270            ! Doing it here also means that asselin filter contribution is removed 
[11050]271            zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
272            zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)   
[4990]273            DO jk = 2, jpkm1
[11050]274               zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
275               zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)   
[4370]276            END DO
277            DO jk = 1, jpkm1
[11050]278               puu(:,:,jk,Kmm) = puu(:,:,jk,Kmm) - (zue(:,:) * r1_hu(:,:,Kmm) - uu_b(:,:,Kmm)) * umask(:,:,jk)
279               pvv(:,:,jk,Kmm) = pvv(:,:,jk,Kmm) - (zve(:,:) * r1_hv(:,:,Kmm) - vv_b(:,:,Kmm)) * vmask(:,:,jk)
[4292]280            END DO
281         ENDIF
282         !
[12489]283      ENDIF ! .NOT. l_1st_euler
[4354]284      !
285      ! Set "now" and "before" barotropic velocities for next time step:
286      ! JC: Would be more clever to swap variables than to make a full vertical
287      ! integration
288      !
[6140]289      IF(.NOT.ln_linssh ) THEN
[11050]290         hu(:,:,Kmm) = pe3u(:,:,1,Kmm ) * umask(:,:,1)
291         hv(:,:,Kmm) = pe3v(:,:,1,Kmm ) * vmask(:,:,1)
[6140]292         DO jk = 2, jpkm1
[11050]293            hu(:,:,Kmm) = hu(:,:,Kmm) + pe3u(:,:,jk,Kmm ) * umask(:,:,jk)
294            hv(:,:,Kmm) = hv(:,:,Kmm) + pe3v(:,:,jk,Kmm ) * vmask(:,:,jk)
[4354]295         END DO
[11050]296         r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
297         r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
[4354]298      ENDIF
299      !
[11050]300      uu_b(:,:,Kaa) = pe3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
301      uu_b(:,:,Kmm) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1)
302      vv_b(:,:,Kaa) = pe3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
303      vv_b(:,:,Kmm) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1)
[6140]304      DO jk = 2, jpkm1
[11050]305         uu_b(:,:,Kaa) = uu_b(:,:,Kaa) + pe3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
306         uu_b(:,:,Kmm) = uu_b(:,:,Kmm) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk)
307         vv_b(:,:,Kaa) = vv_b(:,:,Kaa) + pe3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
308         vv_b(:,:,Kmm) = vv_b(:,:,Kmm) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk)
[4354]309      END DO
[11050]310      uu_b(:,:,Kaa) = uu_b(:,:,Kaa) * r1_hu(:,:,Kaa)
311      vv_b(:,:,Kaa) = vv_b(:,:,Kaa) * r1_hv(:,:,Kaa)
312      uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
313      vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
[4354]314      !
[6140]315      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
[11050]316         CALL iom_put(  "ubar", uu_b(:,:,Kmm) )
317         CALL iom_put(  "vbar", vv_b(:,:,Kmm) )
[6140]318      ENDIF
[4990]319      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
[11050]320         zua(:,:,:) = ( puu(:,:,:,Kmm) - zua(:,:,:) ) * z1_2dt
321         zva(:,:,:) = ( pvv(:,:,:,Kmm) - zva(:,:,:) ) * z1_2dt
[10946]322         CALL trd_dyn( zua, zva, jpdyn_atf, kt, Kmm )
[4990]323      ENDIF
324      !
[13710]325      IF ( iom_use("utau") ) THEN
326         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
327            ALLOCATE(zutau(jpi,jpj)) 
328            DO_2D( 0, 0, 0, 0 )
329               jk = miku(ji,jj) 
330               zutau(ji,jj) = utau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * puu(ji,jj,jk,Kaa)
331            END_2D
332            CALL iom_put(  "utau", zutau(:,:) )
333            DEALLOCATE(zutau)
334         ELSE
335            CALL iom_put(  "utau", utau(:,:) )
336         ENDIF
337      ENDIF
338      !
339      IF ( iom_use("vtau") ) THEN
340         IF ( ln_drgice_imp.OR.ln_isfcav ) THEN
341            ALLOCATE(zvtau(jpi,jpj))
342            DO_2D( 0, 0, 0, 0 )
343               jk = mikv(ji,jj)
344               zvtau(ji,jj) = vtau(ji,jj) + 0.5_wp * rho0 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * pvv(ji,jj,jk,Kaa)
345            END_2D
346            CALL iom_put(  "vtau", zvtau(:,:) )
347            DEALLOCATE(zvtau)
348         ELSE
349            CALL iom_put(  "vtau", vtau(:,:) )
350         ENDIF
351      ENDIF
352      !
[12236]353      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt  - puu(:,:,:,Kaa): ', mask1=umask,   &
354         &                                  tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): '       , mask2=vmask )
[6140]355      !
[9019]356      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
357      IF( l_trddyn     )   DEALLOCATE( zua, zva )
[11050]358      IF( ln_timing    )   CALL timing_stop('dyn_atf')
[2715]359      !
[11050]360   END SUBROUTINE dyn_atf
[3]361
[13710]362#endif
363
[1502]364   !!=========================================================================
[11050]365END MODULE dynatf
Note: See TracBrowser for help on using the repository browser.