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.
tranxt.F90 in branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA – NEMO

source: branches/DEV_r1837_MLF/NEMO/OPA_SRC/TRA/tranxt.F90 @ 1870

Last change on this file since 1870 was 1870, checked in by gm, 14 years ago

ticket: #663 step-1 : introduce the modified forcing term

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 23.4 KB
Line 
1MODULE tranxt
2   !!======================================================================
3   !!                       ***  MODULE  tranxt  ***
4   !! Ocean active tracers:  time stepping on temperature and salinity
5   !!======================================================================
6   !! History :  OPA  !  1991-11  (G. Madec)  Original code
7   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
8   !!            8.0  !  1996-02  (G. Madec & M. Imbard)  opa release 8.0
9   !!             -   !  1996-04  (A. Weaver)  Euler forward step
10   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
11   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
12   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
13   !!             -   !  2005-04  (C. Deltel) Add Asselin trend in the ML budget
14   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
15   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
16   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
17   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA
18   !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   tra_nxt       : time stepping on temperature and salinity
22   !!   tra_nxt_fix   : time stepping on temperature and salinity : fixed    volume case
23   !!   tra_nxt_vvl   : time stepping on temperature and salinity : variable volume case
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers variables
26   USE dom_oce         ! ocean space and time domain variables
27   USE zdf_oce         ! ???
28   USE domvvl          ! variable volume
29   USE dynspg_oce      ! surface     pressure gradient variables
30   USE dynhpg          ! hydrostatic pressure gradient
31   USE trdmod_oce      ! ocean variables trends
32   USE trdmod          ! ocean active tracers trends
33   USE phycst
34   USE obctra          ! open boundary condition (obc_tra routine)
35   USE bdytra          ! Unstructured open boundary condition (bdy_tra routine)
36   USE in_out_manager  ! I/O manager
37   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
38   USE prtctl          ! Print control
39   USE agrif_opa_update
40   USE agrif_opa_interp
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   tra_nxt    ! routine called by step.F90
46
47   REAL(wp)                 ::   rbcp, r1_2bcp   ! Brown & Campana parameters for semi-implicit hpg
48   REAL(wp), DIMENSION(jpk) ::   r2dt_t          ! vertical profile time step, =2*rdttra (leapfrog) or =rdttra (Euler)
49
50   !! * Substitutions
51#  include "domzgr_substitute.h90"
52   !!----------------------------------------------------------------------
53   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
54   !! $Id$
55   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
56   !!----------------------------------------------------------------------
57
58CONTAINS
59
60   SUBROUTINE tra_nxt( kt )
61      !!----------------------------------------------------------------------
62      !!                   ***  ROUTINE tranxt  ***
63      !!
64      !! ** Purpose :   Apply the boundary condition on the after temperature 
65      !!             and salinity fields, achieved the time stepping by adding
66      !!             the Asselin filter on now fields and swapping the fields.
67      !!
68      !! ** Method  :   At this stage of the computation, ta and sa are the
69      !!             after temperature and salinity as the time stepping has
70      !!             been performed in trazdf_imp or trazdf_exp module.
71      !!
72      !!              - Apply lateral boundary conditions on (ta,sa)
73      !!             at the local domain   boundaries through lbc_lnk call,
74      !!             at the radiative open boundaries (lk_obc=T),
75      !!             at the relaxed   open boundaries (lk_bdy=T), and
76      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
77      !!
78      !!              - Update lateral boundary conditions on AGRIF children
79      !!             domains (lk_agrif=T)
80      !!
81      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
82      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
83      !!
84      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
85      !!----------------------------------------------------------------------
86      USE oce, ONLY :    ztrdt => ua   ! use ua as 3D workspace   
87      USE oce, ONLY :    ztrds => va   ! use va as 3D workspace   
88      !!
89      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
90      !!
91      INTEGER  ::   jk      ! dummy loop indices
92      REAL(wp) ::   zfact   ! temporary scalars
93      !!----------------------------------------------------------------------
94
95      IF( kt == nit000 ) THEN          !==  initialisation  ==!
96         IF(lwp) WRITE(numout,*)
97         IF(lwp) WRITE(numout,*) 'tra_nxt : achieve the time stepping by Asselin filter and array swap'
98         IF(lwp) WRITE(numout,*) '~~~~~~~'
99         !
100         rbcp    = 0.25 * (1 + atfp) * (1 + atfp * atfp)       ! Brown & Campana parameter for semi-implicit hpg
101         r1_2bcp = 1.e0 - 2.e0 * rbcp
102      ENDIF
103      !                                      ! set time step size (Euler/Leapfrog)
104      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt_t(:) =     rdttra(:)      ! at nit000             (Euler)
105      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt_t(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
106      ENDIF
107
108      !                                !==  Update after tracer on domain lateral boundaries  ==!
109      !
110      CALL lbc_lnk( ta, 'T', 1. )            ! local domain boundaries  (T-point, unchanged sign)
111      CALL lbc_lnk( sa, 'T', 1. )
112      !
113#if defined key_obc
114      CALL obc_tra( kt )                     ! OBC open boundaries
115#endif
116#if defined key_bdy
117      CALL bdy_tra( kt )                     ! BDY open boundaries
118#endif
119#if defined key_agrif
120      CALL Agrif_tra                         ! AGRIF zoom boundaries
121#endif
122 
123      IF( l_trdtra ) THEN              ! trends computation: store now fields before applying the Asselin filter
124         ztrdt(:,:,:) = tn(:,:,:)     
125         ztrds(:,:,:) = sn(:,:,:)
126      ENDIF
127
128      !                                !==  modifed Leap-Frog + Asselin filter (modified LF-RA) scheme  ==!
129      IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl( kt )      ! variable volume level (vvl)
130      ELSE                  ;   CALL tra_nxt_fix( kt )      ! fixed    volume level
131      ENDIF
132
133#if defined key_agrif
134      IF( .NOT.Agrif_Root() )    CALL Agrif_Update_Tra( kt )   ! Update tracer at AGRIF zoom boundaries (children only)
135#endif     
136
137      IF( l_trdtra ) THEN              ! trends computation: trend of the Asselin filter (tb filtered - tb)/dt     
138         DO jk = 1, jpkm1
139            zfact = 1.e0 / r2dt_t(jk)             
140            ztrdt(:,:,jk) = ( tb(:,:,jk) - ztrdt(:,:,jk) ) * zfact
141            ztrds(:,:,jk) = ( sb(:,:,jk) - ztrds(:,:,jk) ) * zfact
142         END DO
143         CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt )
144      END IF
145
146      !                        ! control print
147      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt  - Tn: ', mask1=tmask,   &
148         &                       tab3d_2=sn, clinfo2=       ' Sn: ', mask2=tmask )
149      !
150   END SUBROUTINE tra_nxt
151
152
153   SUBROUTINE tra_nxt_fix( kt )
154      !!----------------------------------------------------------------------
155      !!                   ***  ROUTINE tra_nxt_fix  ***
156      !!
157      !! ** Purpose :   fixed volume: apply the Asselin time filter and
158      !!                swap the tracer fields.
159      !!
160      !! ** Method  : - Apply a Asselin time filter on now fields.
161      !!              - save in (ta,sa) an average over the three time levels
162      !!             which will be used to compute rdn and thus the semi-implicit
163      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
164      !!              - swap tracer fields to prepare the next time_step.
165      !!                This can be summurized for temperature as:
166      !!             ztm = (rbcp*ta+(2-rbcp)*tn+rbcp*tb)       ln_dynhpg_imp = T
167      !!             ztm = 0                                   otherwise
168      !!                   with rbcp=(1+atfp)*(1+atfp*atfp)/4
169      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
170      !!             tn  = ta
171      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
172      !!
173      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
174      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
175      !!----------------------------------------------------------------------
176      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
177      !!
178      INTEGER  ::   ji, jj, jk      ! dummy loop indices
179      REAL(wp) ::   ztm, ztf, zac   ! temporary scalars
180      REAL(wp) ::   zsm, zsf        !    -         -
181      !!----------------------------------------------------------------------
182
183      IF( kt == nit000 ) THEN
184         IF(lwp) WRITE(numout,*)
185         IF(lwp) WRITE(numout,*) 'tra_nxt_fix : time stepping'
186         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
187      ENDIF
188      !
189      !                                              ! ----------------------- !
190      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !
191         !                                           ! ----------------------- !
192         !
193         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
194            DO jk = 1, jpkm1                              ! (only swap)
195               tn(:,:,jk) = ta(:,:,jk)                                                          ! tb <-- tn
196               sn(:,:,jk) = sa(:,:,jk)
197            END DO
198         ELSE                                             ! general case (Leapfrog + Asselin filter)
199            DO jk = 1, jpkm1
200               DO jj = 1, jpj
201                  DO ji = 1, jpi
202                     ztm = rbcp * ( ta(ji,jj,jk) + tb(ji,jj,jk) ) + r1_2bcp * tn(ji,jj,jk)      ! mean t
203                     zsm = rbcp * ( sa(ji,jj,jk) + sb(ji,jj,jk) ) + r1_2bcp * sn(ji,jj,jk) 
204                     ztf = atfp * ( ta(ji,jj,jk) + tb(ji,jj,jk)   - 2.      * tn(ji,jj,jk) )    ! Asselin filter on t
205                     zsf = atfp * ( sa(ji,jj,jk) + sb(ji,jj,jk)   - 2.      * sn(ji,jj,jk) )
206                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                          ! tb <-- filtered tn
207                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
208                     tn(ji,jj,jk) = ta(ji,jj,jk)                                                ! tn <-- ta
209                     sn(ji,jj,jk) = sa(ji,jj,jk)
210                     ta(ji,jj,jk) = ztm                                                         ! ta <-- mean t
211                     sa(ji,jj,jk) = zsm
212                  END DO
213               END DO
214            END DO
215         ENDIF
216         !                                           ! ----------------------- !
217      ELSE                                           !    explicit hpg case    !
218         !                                           ! ----------------------- !
219         !
220         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
221            DO jk = 1, jpkm1
222               tn(:,:,jk) = ta(:,:,jk)                                                    ! tn <-- ta
223               sn(:,:,jk) = sa(:,:,jk)
224            END DO
225         ELSE                                             ! general case (Leapfrog + Asselin filter)
226            DO jk = 1, jpkm1
227               DO jj = 1, jpj
228                  DO ji = 1, jpi
229                     ztf = atfp * ( ta(ji,jj,jk) - 2.* tn(ji,jj,jk) + tb(ji,jj,jk) )      ! Asselin filter on t
230                     zsf = atfp * ( sa(ji,jj,jk) - 2.* sn(ji,jj,jk) + sb(ji,jj,jk) )
231                     tb(ji,jj,jk) = tn(ji,jj,jk) + ztf                                    ! tb <-- filtered tn
232                     sb(ji,jj,jk) = sn(ji,jj,jk) + zsf
233                     tn(ji,jj,jk) = ta(ji,jj,jk)                                          ! tn <-- ta
234                     sn(ji,jj,jk) = sa(ji,jj,jk)
235                  END DO
236               END DO
237            END DO
238         ENDIF
239      ENDIF
240      !
241!!gm from Matthieu : unchecked
242      IF( neuler /= 0 .OR. kt /= nit000 ) THEN           ! remove the forcing from the Asselin filter
243         zac = atfp * rdttra(1)
244         tb(:,:,1) = tb(:,:,1) - zac * ( qns (:,:) - qns_b (:,:) )         ! non solar surface flux
245         sb(:,:,1) = sn(:,:,1) - zac * ( emps(:,:) - emps_b(:,:) )         ! surface salt flux
246         !
247         IF( ln_traqsr )                                                   ! solar penetrating flux
248            DO jk = 1, nksr
249                  DO jj = 1, jpj
250                     DO ji = 1, jpi
251                        IF( ln_sco ) THEN
252                           z_cor_qsr = etot3(ji,jj,jk) * ( qsr(ji,jj) - qsrb(ji,jj) )
253                        ELSEIF( ln_zps .OR. ln_zco ) THEN
254                           z_cor_qsr = ( qsr(ji,jj) - qsrb(ji,jj) ) *   &
255                                &   ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) )
256                        ENDIF
257                        zt = zt - zfact1 * z_cor_qsr
258                     END DO
259                  END DO
260      ENDIF
261      !
262   END SUBROUTINE tra_nxt_fix
263
264
265   SUBROUTINE tra_nxt_vvl( kt )
266      !!----------------------------------------------------------------------
267      !!                   ***  ROUTINE tra_nxt_vvl  ***
268      !!
269      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
270      !!                and swap the tracer fields.
271      !!
272      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
273      !!              - save in (ta,sa) a thickness weighted average over the three
274      !!             time levels which will be used to compute rdn and thus the semi-
275      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
276      !!              - swap tracer fields to prepare the next time_step.
277      !!                This can be summurized for temperature as:
278      !!             ztm = ( rbcp*e3t_a*ta + (2-rbtp)*e3t_n*tn + rbcp*e3t_b*tb )   ln_dynhpg_imp = T
279      !!                  /( rbcp*e3t_a    + (2-rbcp)*e3t_n    + rbcp*e3t_b    )   
280      !!             ztm = 0                                                       otherwise
281      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
282      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
283      !!             tn  = ta
284      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
285      !!
286      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
287      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
288      !!----------------------------------------------------------------------
289      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
290      !!     
291      INTEGER  ::   ji, jj, jk             ! dummy loop indices
292      REAL(wp) ::   ztm , ztc_f , ztf , ztca, ztcn, ztcb   ! temporary scalar
293      REAL(wp) ::   zsm , zsc_f , zsf , zsca, zscn, zscb   !    -         -
294      REAL(wp) ::   ze3mr, ze3fr                           !    -         -
295      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f         !    -         -
296      !!----------------------------------------------------------------------
297
298      IF( kt == nit000 ) THEN
299         IF(lwp) WRITE(numout,*)
300         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl : time stepping'
301         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
302      ENDIF
303
304      !                                              ! ----------------------- !
305      IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg case !  (mean tracer computation)
306         !                                           ! ----------------------- !
307         !
308         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! Euler time-stepping at first time-step
309            DO jk = 1, jpkm1                                   ! (only swap)
310               tn(:,:,jk) = ta(:,:,jk)                         ! tn <-- ta
311               sn(:,:,jk) = sa(:,:,jk)
312            END DO
313            !                                             ! general case (Leapfrog + Asselin filter)
314         ELSE                                             ! apply filter on thickness weighted tracer and swap
315            DO jk = 1, jpkm1
316               DO jj = 1, jpj
317                  DO ji = 1, jpi
318                     ze3t_b = fse3t_b(ji,jj,jk)
319                     ze3t_n = fse3t_n(ji,jj,jk)
320                     ze3t_a = fse3t_a(ji,jj,jk)
321                     !                                         ! tracer content at Before, now and after
322                     ztcb = tb(ji,jj,jk) *  ze3t_b   ;   zscb = sb(ji,jj,jk) * ze3t_b
323                     ztcn = tn(ji,jj,jk) *  ze3t_n   ;   zscn = sn(ji,jj,jk) * ze3t_n
324                     ztca = ta(ji,jj,jk) *  ze3t_a   ;   zsca = sa(ji,jj,jk) * ze3t_a
325                     !
326                     !                                         ! Asselin filter on thickness and tracer content
327                     ze3t_f = atfp * ( ze3t_a + ze3t_b - 2.* ze3t_n )
328                     ztc_f  = atfp * ( ztca   + ztcb   - 2.* ztcn   ) 
329                     zsc_f  = atfp * ( zsca   + zscb   - 2.* zscn   ) 
330                     !
331                     !                                         ! filtered tracer including the correction
332                     ze3fr = 1.e0  / ( ze3t_n + ze3t_f )
333                     ztf   = ze3fr * ( ztcn   + ztc_f  )
334                     zsf   = ze3fr * ( zscn   + zsc_f  )
335                     !                                         ! mean thickness and tracer
336                     ze3mr = 1.e0  / (  rbcp * ( ze3t_a + ze3t_b ) + r1_2bcp * ze3t_n )
337                     ztm   = ze3mr * (  rbcp * ( ztca   + ztcb   ) + r1_2bcp * ztcn   )
338                     zsm   = ze3mr * (  rbcp * ( zsca   + zscb   ) + r1_2bcp * zscn   )
339!!gm mean e3t have to be saved and used in dynhpg  or it can be recomputed in dynhpg !!
340!!gm                 e3t_m(ji,jj,jk) = 1.e0 / ze3mr
341                     !                                         ! swap of arrays
342                     tb(ji,jj,jk) = ztf                             ! tb <-- tn + filter
343                     sb(ji,jj,jk) = zsf
344                     tn(ji,jj,jk) = ta(ji,jj,jk)                    ! tn <-- ta
345                     sn(ji,jj,jk) = sa(ji,jj,jk)
346                     ta(ji,jj,jk) = ztm                             ! ta <-- mean t
347                     sa(ji,jj,jk) = zsm
348                  END DO
349               END DO
350            END DO
351         ENDIF
352         !                                           ! ----------------------- !
353      ELSE                                           !    explicit hpg case    !  (NO mean tracer computation)
354         !                                           ! ----------------------- !
355         !
356         IF( neuler == 0 .AND. kt == nit000 ) THEN        ! case of Euler time-stepping at first time-step
357            DO jk = 1, jpkm1                                   ! (only swap)
358               tn(:,:,jk) = ta(:,:,jk)                         ! tn <-- ta
359               sn(:,:,jk) = sa(:,:,jk)
360            END DO
361            !                                             ! general case (Leapfrog + Asselin filter)
362         ELSE                                             ! apply filter on thickness weighted tracer and swap
363            DO jk = 1, jpkm1
364               DO jj = 1, jpj
365                  DO ji = 1, jpi
366                     ze3t_b = fse3t_b(ji,jj,jk)
367                     ze3t_n = fse3t_n(ji,jj,jk)
368                     ze3t_a = fse3t_a(ji,jj,jk)
369                     !                                         ! tracer content at Before, now and after
370                     ztcb = tb(ji,jj,jk) *  ze3t_b   ;   zscb = sb(ji,jj,jk) * ze3t_b
371                     ztcn = tn(ji,jj,jk) *  ze3t_n   ;   zscn = sn(ji,jj,jk) * ze3t_n
372                     ztca = ta(ji,jj,jk) *  ze3t_a   ;   zsca = sa(ji,jj,jk) * ze3t_a
373                     !
374                     !                                         ! Asselin filter on thickness and tracer content
375                     ze3t_f = atfp * ( ze3t_a - 2.* ze3t_n + ze3t_b )
376                     ztc_f  = atfp * ( ztca   - 2.* ztcn   + ztcb   ) 
377                     zsc_f  = atfp * ( zsca   - 2.* zscn   + zscb   ) 
378                     !
379                     !                                         ! filtered tracer including the correction
380                     ze3fr = 1.e0 / ( ze3t_n + ze3t_f )
381                     ztf   =  ( ztcn  + ztc_f ) * ze3fr
382                     zsf   =  ( zscn  + zsc_f ) * ze3fr
383                     !                                         ! swap of arrays
384                     tb(ji,jj,jk) = ztf                             ! tb <-- tn filtered
385                     sb(ji,jj,jk) = zsf
386                     tn(ji,jj,jk) = ta(ji,jj,jk)                    ! tn <-- ta
387                     sn(ji,jj,jk) = sa(ji,jj,jk)
388                  END DO
389               END DO
390            END DO
391         ENDIF
392      ENDIF
393!!gm from Matthieu : unchecked
394      !
395      IF( neuler /= 0 .OR. kt /= nit000 ) THEN           ! remove the forcing from the Asselin filter
396               IF( lk_vvl ) THEN                                        ! Varying levels
397                  DO jj = 1, jpj
398                     DO ji = 1, jpi
399                        ! - ML - modification for varaiance diagnostics
400                        zssh1 = tmask(ji,jj,jk) / ( fse3t(ji,jj,jk) + atfp * ( zfse3ta(ji,jj,jk) - 2*fse3t(ji,jj,jk) &
401                           &                                                 + fse3tb(ji,jj,jk) ) )
402                        zt = zssh1 * ( tn(ji,jj,jk) +  atfp * ( ta(ji,jj,jk) - 2 * tn(ji,jj,jk) + tb(ji,jj,jk) ) )
403                        zs = zssh1 * ( sn(ji,jj,jk) +  atfp * ( sa(ji,jj,jk) - 2 * sn(ji,jj,jk) + sb(ji,jj,jk) ) )
404                        ! filter correction for global conservation
405                        !------------------------------------------
406                        zfact1 = atfp * rdttra(1) * zssh1
407                        IF (jk == 1) THEN                        ! remove sbc trend from time filter
408                           zt = zt - zfact1 * ( qn(ji,jj) - qb(ji,jj) )
409!!???                           zs = zs - zfact1 * ( - emp( ji,jj) + empb(ji,jj) ) * zsrau * 35.5e0
410                        ENDIF
411                        ! remove qsr trend from time filter
412                        IF (jk <= nksr) THEN
413                           IF( ln_sco ) THEN
414                              z_cor_qsr = etot3(ji,jj,jk) * ( qsr(ji,jj) - qsrb(ji,jj) )
415                           ELSEIF( ln_zps .OR. ln_zco ) THEN
416                              z_cor_qsr = ( qsr(ji,jj) - qsrb(ji,jj) ) *   &
417                                   &   ( gdsr(jk) * tmask(ji,jj,jk) - gdsr(jk+1) * tmask(ji,jj,jk+1) )
418                           ENDIF
419                           zt = zt - zfact1 * z_cor_qsr
420                           IF (jk == nksr) qsrb(ji,jj) = qsr(ji,jj)
421                        ENDIF
422                        ! - ML - end of modification
423                        zssh1 = tmask(ji,jj,1) / zfse3ta(ji,jj,jk)
424                        tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1
425                        sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1
426                     END DO
427                  END DO
428      ENDIF
429!!gm end
430      !
431   END SUBROUTINE tra_nxt_vvl
432
433   !!======================================================================
434END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.