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 trunk/NEMO/OPA_SRC/TRA – NEMO

source: trunk/NEMO/OPA_SRC/TRA/tranxt.F90 @ 888

Last change on this file since 888 was 888, checked in by ctlod, 16 years ago

merge dev_001_SBC branche with the trunk to include the New Surface Module package, see ticket: #113
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.6 KB
Line 
1MODULE tranxt
2   !!======================================================================
3   !!                       ***  MODULE  tranxt  ***
4   !! Ocean active tracers:  time stepping on temperature and salinity
5   !!======================================================================
6   !! History :  7.0  !  91-11  (G. Madec)  Original code
7   !!                 !  93-03  (M. Guyon)  symetrical conditions
8   !!                 !  96-02  (G. Madec & M. Imbard)  opa release 8.0
9   !!            8.0  !  96-04  (A. Weaver)  Euler forward step
10   !!            8.2  !  99-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
11   !!            8.5  !  02-08  (G. Madec)  F90: Free form and module
12   !!                 !  02-11  (C. Talandier, A-M Treguier) Open boundaries
13   !!                 !  05-04  (C. Deltel) Add Asselin trend in the ML budget
14   !!            9.0  !  06-02  (L. Debreu, C. Mazauric) Agrif implementation
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   tra_nxt     : time stepping on temperature and salinity
19   !!----------------------------------------------------------------------
20   USE oce             ! ocean dynamics and tracers variables
21   USE dom_oce         ! ocean space and time domain variables
22   USE zdf_oce         ! ???
23   USE dynspg_oce      ! surface pressure gradient variables
24   USE trdmod_oce      ! ocean variables trends
25   USE trdmod          ! ocean active tracers trends
26   USE phycst
27   USE domvvl          ! variable volume
28   USE obctra          ! open boundary condition (obc_tra routine)
29   USE in_out_manager  ! I/O manager
30   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
31   USE prtctl          ! Print control
32   USE agrif_opa_update
33   USE agrif_opa_interp
34
35
36   IMPLICIT NONE
37   PRIVATE
38
39   !! * Routine accessibility
40   PUBLIC   tra_nxt   ! routine called by step.F90
41
42   REAL(wp) ::   vemp ! total amount of volume added or removed by E-P forcing
43
44   !! * Substitutions
45#  include "domzgr_substitute.h90"
46   !!----------------------------------------------------------------------
47   !!   OPA 9.0 , LOCEAN-IPSL (2006)
48   !! $Id$
49   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51
52CONTAINS
53
54   SUBROUTINE tra_nxt( kt )
55      !!----------------------------------------------------------------------
56      !!                   ***  ROUTINE tranxt  ***
57      !!
58      !! ** Purpose :   Compute the temperature and salinity fields at the
59      !!      next time-step from their temporal trends and swap the fields.
60      !!
61      !! ** Method  :   Apply lateral boundary conditions on (ua,va) through
62      !!      call to lbc_lnk routine
63      !!      After t and s are compute using a leap-frog scheme environment:
64      !!         ta = tb + 2 rdttra(k) * ta
65      !!         sa = sb + 2 rdttra(k) * sa
66      !!      Compute and save in (ta,sa) an average over three time levels
67      !!      (before,now and after) of temperature and salinity which is
68      !!      used to compute rhd in eos routine and thus the hydrostatic
69      !!      pressure gradient (ln_dynhpg_imp = T)
70      !!      Apply an Asselin time filter on now tracers (tn,sn) to avoid
71      !!      the divergence of two consecutive time-steps and swap tracer
72      !!      arrays to prepare the next time_step:
73      !!         (zt,zs) = (ta+2tn+tb,sa+2sn+sb)/4       (ln_dynhpg_imp = T)
74      !!         (zt,zs) = (0,0)                            (default option)
75      !!         (tb,sb) = (tn,vn) + atfp [ (tb,sb) + (ta,sa) - 2 (tn,sn) ]
76      !!         (tn,sn) = (ta,sa)
77      !!         (ta,sa) = (zt,zs)  (NB: reset to 0 after use in eos.F)
78      !!
79      !! ** Action  : - update (tb,sb) and (tn,sn)
80      !!              - (ta,sa) time averaged (t,s)      (ln_dynhpg_imp = T)
81      !!----------------------------------------------------------------------
82      USE oce, ONLY :    ztrdt => ua   ! use ua as 3D workspace   
83      USE oce, ONLY :    ztrds => va   ! use va as 3D workspace   
84      !!
85      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
86      !!
87      INTEGER  ::   ji, jj, jk    ! dummy loop indices
88      REAL(wp) ::   zt, zs, zssh1 ! temporary scalars
89      REAL(wp) ::   zfact         ! temporary scalar
90      !! Variable volume
91      REAL(wp), DIMENSION(jpi,jpj) ::   zssh                         ! temporary scalars
92      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfse3tb, zfse3tn, zfse3ta  ! 3D workspace
93
94      !!----------------------------------------------------------------------
95
96      !! Explicit physics with thickness weighted updates
97      IF( lk_vvl .AND. ln_zdfexp ) THEN
98
99         ! Scale factors at before and after time step
100         ! -------------------------------------------
101         CALL dom_vvl_sf( sshb, 'T', zfse3tb ) ;    CALL dom_vvl_sf( ssha, 'T', zfse3ta )
102
103         ! Asselin filtered scale factor at now time step
104         ! ----------------------------------------------
105         IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN
106            CALL dom_vvl_sf_ini( 'T', zfse3tn ) 
107         ELSE
108            zssh(:,:) = atfp * ( sshb(:,:) + ssha(:,:) ) + atfp1 * sshn(:,:)
109            CALL dom_vvl_sf( zssh, 'T', zfse3tn ) 
110         ENDIF
111
112         ! Thickness weighting
113         ! -------------------
114         DO jk = 1, jpkm1
115            DO jj = 1, jpj
116               DO ji = 1, jpi
117                  ta(ji,jj,jk) = ta(ji,jj,jk) * fse3t(ji,jj,jk)
118                  sa(ji,jj,jk) = sa(ji,jj,jk) * fse3t(ji,jj,jk)
119
120                  tn(ji,jj,jk) = tn(ji,jj,jk) * fse3t(ji,jj,jk)
121                  sn(ji,jj,jk) = sn(ji,jj,jk) * fse3t(ji,jj,jk)
122
123                  tb(ji,jj,jk) = tb(ji,jj,jk) * zfse3tb(ji,jj,jk)
124                  sb(ji,jj,jk) = sb(ji,jj,jk) * zfse3tb(ji,jj,jk)
125               END DO
126            END DO
127         END DO
128
129      ENDIF
130
131      IF( l_trdtra ) THEN
132         ztrdt(:,:,jpk) = 0.e0
133         ztrds(:,:,jpk) = 0.e0
134      ENDIF
135
136      ! 0. Lateral boundary conditions on ( ta, sa )   (T-point, unchanged sign)
137      ! ---------------------------------============
138      CALL lbc_lnk( ta, 'T', 1. )   
139      CALL lbc_lnk( sa, 'T', 1. )
140
141      !                                                ! ===============
142      DO jk = 1, jpkm1                                 ! Horizontal slab
143         !                                             ! ===============
144
145         ! 1. Leap-frog scheme (only in explicit case, otherwise the
146         ! -------------------  time stepping is already done in trazdf)
147         IF( ln_zdfexp ) THEN
148            zfact = 2. * rdttra(jk)
149            IF( neuler == 0 .AND. kt == nit000 )   zfact = rdttra(jk)
150            ta(:,:,jk) = ( tb(:,:,jk) + zfact * ta(:,:,jk) ) * tmask(:,:,jk)
151            sa(:,:,jk) = ( sb(:,:,jk) + zfact * sa(:,:,jk) ) * tmask(:,:,jk)
152            IF(l_trdtra)   CALL ctl_stop( 'tranxt: Asselin ML trend not yet accounted for.' )
153         ENDIF
154
155#if defined key_obc
156         !                                             ! ===============
157      END DO                                           !   End of slab
158      !                                                ! ===============
159      ! Update tracers on open boundaries.
160      CALL obc_tra( kt )
161      !                                                ! ===============
162      DO jk = 1, jpkm1                                 ! Horizontal slab
163         !                                             ! ===============
164#endif
165#if defined key_agrif
166         !                                             ! ===============
167      END DO                                           !   End of slab
168      !                                                ! ===============
169      ! Update tracers on open boundaries.
170      CALL Agrif_tra
171      !                                                ! ===============
172      DO jk = 1, jpkm1                                 ! Horizontal slab
173         !                                             ! ===============
174#endif
175         ! 2. Time filter and swap of arrays
176         ! ---------------------------------
177
178         IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg
179            IF( neuler == 0 .AND. kt == nit000 ) THEN
180               DO jj = 1, jpj
181                  DO ji = 1, jpi
182                     zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25
183                     zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25
184                     tb(ji,jj,jk) = tn(ji,jj,jk)
185                     sb(ji,jj,jk) = sn(ji,jj,jk)
186                     tn(ji,jj,jk) = ta(ji,jj,jk)
187                     sn(ji,jj,jk) = sa(ji,jj,jk)
188                     ta(ji,jj,jk) = zt
189                     sa(ji,jj,jk) = zs
190                  END DO
191               END DO
192               IF( l_trdtra ) THEN
193                  ztrdt(:,:,jk) = 0.e0
194                  ztrds(:,:,jk) = 0.e0
195               END IF
196            ELSE
197               DO jj = 1, jpj
198                  DO ji = 1, jpi
199                     zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25
200                     zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25
201                     tb(ji,jj,jk) = atfp  * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk)
202                     sb(ji,jj,jk) = atfp  * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk)
203                     IF( l_trdtra ) THEN ! ChD ceci est a optimiser, mais ca marche
204                        ztrdt(ji,jj,jk) = tb(ji,jj,jk) - tn(ji,jj,jk)
205                        ztrds(ji,jj,jk) = sb(ji,jj,jk) - sn(ji,jj,jk)
206                     END IF
207                     tn(ji,jj,jk) = ta(ji,jj,jk)
208                     sn(ji,jj,jk) = sa(ji,jj,jk)
209                     ta(ji,jj,jk) = zt
210                     sa(ji,jj,jk) = zs
211                  END DO
212               END DO
213            ENDIF
214         ELSE                                          ! Default case
215            IF( neuler == 0 .AND. kt == nit000 ) THEN
216               IF( (lk_vvl .AND. ln_zdfexp) ) THEN                      ! Varying levels
217                  DO jj = 1, jpj
218                     DO ji = 1, jpi
219                        zssh1 = tmask(ji,jj,jk) / fse3t(ji,jj,jk)
220                        tb(ji,jj,jk) = tn(ji,jj,jk) * zssh1 * tmask(ji,jj,jk)
221                        sb(ji,jj,jk) = sn(ji,jj,jk) * zssh1 * tmask(ji,jj,jk)
222                        zssh1 = tmask(ji,jj,jk) / zfse3ta(ji,jj,jk)
223                        tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1 * tmask(ji,jj,jk)
224                        sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1 * tmask(ji,jj,jk)
225                     END DO
226                  END DO
227               ELSE                                                     ! Fixed levels
228                 DO jj = 1, jpj
229                     DO ji = 1, jpi
230                        tb(ji,jj,jk) = tn(ji,jj,jk)
231                        sb(ji,jj,jk) = sn(ji,jj,jk)
232                        tn(ji,jj,jk) = ta(ji,jj,jk)
233                        sn(ji,jj,jk) = sa(ji,jj,jk)
234                     END DO
235                  END DO
236               ENDIF
237               IF( l_trdtra ) THEN
238                  ztrdt(:,:,jk) = 0.e0
239                  ztrds(:,:,jk) = 0.e0
240               END IF
241            ELSE
242               IF( l_trdtra ) THEN
243                  DO jj = 1, jpj
244                     DO ji = 1, jpi
245                        ztrdt(ji,jj,jk) = atfp * ( tb(ji,jj,jk) - 2*tn(ji,jj,jk) + ta(ji,jj,jk) )
246                        ztrds(ji,jj,jk) = atfp * ( sb(ji,jj,jk) - 2*sn(ji,jj,jk) + sa(ji,jj,jk) )
247                     END DO
248                  END DO
249               END IF
250               IF( (lk_vvl .AND. ln_zdfexp) ) THEN                      ! Varying levels
251                  DO jj = 1, jpj
252                     DO ji = 1, jpi
253                        zssh1 = tmask(ji,jj,jk) / zfse3tn(ji,jj,jk)
254                        tb(ji,jj,jk) = ( atfp  * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) &
255                          &            + atfp1 * tn(ji,jj,jk) ) * zssh1
256                        sb(ji,jj,jk) = ( atfp  * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) &
257                          &            + atfp1 * sn(ji,jj,jk) ) * zssh1
258                        zssh1 = tmask(ji,jj,1) / zfse3ta(ji,jj,jk)
259                        tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1
260                        sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1
261                     END DO
262                  END DO
263               ELSE                                                     ! Fixed levels or first varying level
264                  DO jj = 1, jpj
265                     DO ji = 1, jpi
266                        tb(ji,jj,jk) = atfp  * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk)
267                        sb(ji,jj,jk) = atfp  * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk)
268                        tn(ji,jj,jk) = ta(ji,jj,jk)
269                        sn(ji,jj,jk) = sa(ji,jj,jk)
270                     END DO
271                  END DO
272               ENDIF
273            ENDIF
274         ENDIF
275         !                                             ! ===============
276      END DO                                           !   End of slab
277      !                                                ! ===============
278
279      IF( l_trdtra )   THEN      ! Take the Asselin trend into account
280         ztrdt(:,:,:) = ztrdt(:,:,:) / ( 2.*rdt )
281         ztrds(:,:,:) = ztrds(:,:,:) / ( 2.*rdt )         
282         CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt )
283      END IF
284
285      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt  - Tn: ', mask1=tmask,   &
286         &                       tab3d_2=sn, clinfo2=       ' Sn: ', mask2=tmask )
287#if defined key_agrif
288      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Tra( kt )
289#endif     
290      !
291   END SUBROUTINE tra_nxt
292
293   !!======================================================================
294END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.