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 @ 911

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

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