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

Last change on this file since 642 was 642, checked in by opalod, 17 years ago

nemo_v2_bugfix_028 : CT : correction to be able to compile without compilation stop when using key_zco cpp key
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.7 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 in_out_manager  ! I/O manager
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25   USE obctra          ! open boundary condition (obc_tra routine)
26   USE trdmod          ! ocean active tracers trends
27   USE trdmod_oce      ! ocean variables trends
28   USE prtctl          ! Print control
29   USE agrif_opa_update
30   USE agrif_opa_interp
31
32   USE ocesbc          ! ocean surface boundary condition
33   USE domvvl          ! variable volume
34   USE dynspg_oce      ! surface pressure gradient variables
35   USE phycst
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   !! $Header$
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         zfse3tb(:,:,:) = sfe3( sshb, 'T' )    ;    zfse3ta(:,:,:) = sfe3( ssha, 'T' )
103
104         ! Asselin filtered scale factor at now time step
105         ! ----------------------------------------------
106         IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN
107            zfse3tn(:,:,:) = sfe3ini( 'T' )
108         ELSE
109            zssh(:,:) = atfp * ( sshb(:,:) + ssha(:,:) ) + atfp1 * sshn(:,:)
110            zfse3tn(:,:,:) = sfe3( zssh, 'T' )
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      DO jk = 1, jpkm1                                 ! Horizontal slab
164         !                                             ! ===============
165#endif
166#if defined key_agrif
167         !                                             ! ===============
168      END DO                                           !   End of slab
169      !                                                ! ===============
170      ! Update tracers on open boundaries.
171      CALL Agrif_tra( kt )
172      !                                                ! ===============
173      DO jk = 1, jpkm1                                 ! Horizontal slab
174         !                                             ! ===============
175#endif
176         ! 2. Time filter and swap of arrays
177         ! ---------------------------------
178
179         IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg
180            IF( neuler == 0 .AND. kt == nit000 ) THEN
181               DO jj = 1, jpj
182                  DO ji = 1, jpi
183                     zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25
184                     zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25
185                     tb(ji,jj,jk) = tn(ji,jj,jk)
186                     sb(ji,jj,jk) = sn(ji,jj,jk)
187                     tn(ji,jj,jk) = ta(ji,jj,jk)
188                     sn(ji,jj,jk) = sa(ji,jj,jk)
189                     ta(ji,jj,jk) = zt
190                     sa(ji,jj,jk) = zs
191                  END DO
192               END DO
193               IF( l_trdtra ) THEN
194                  ztrdt(:,:,jk) = 0.e0
195                  ztrds(:,:,jk) = 0.e0
196               END IF
197            ELSE
198               DO jj = 1, jpj
199                  DO ji = 1, jpi
200                     zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25
201                     zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25
202                     tb(ji,jj,jk) = atfp  * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk)
203                     sb(ji,jj,jk) = atfp  * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk)
204                     IF( l_trdtra ) THEN ! ChD ceci est a optimiser, mais ca marche
205                        ztrdt(ji,jj,jk) = tb(ji,jj,jk) - tn(ji,jj,jk)
206                        ztrds(ji,jj,jk) = sb(ji,jj,jk) - sn(ji,jj,jk)
207                     END IF
208                     tn(ji,jj,jk) = ta(ji,jj,jk)
209                     sn(ji,jj,jk) = sa(ji,jj,jk)
210                     ta(ji,jj,jk) = zt
211                     sa(ji,jj,jk) = zs
212                  END DO
213               END DO
214            ENDIF
215         ELSE                                          ! Default case
216            IF( neuler == 0 .AND. kt == nit000 ) THEN
217               IF( (lk_vvl .AND. ln_zdfexp) ) THEN                      ! Varying levels
218                  DO jj = 1, jpj
219                     DO ji = 1, jpi
220                        zssh1 = tmask(ji,jj,jk) / fse3t(ji,jj,jk)
221                        tb(ji,jj,jk) = tn(ji,jj,jk) * zssh1 * tmask(ji,jj,jk)
222                        sb(ji,jj,jk) = sn(ji,jj,jk) * zssh1 * tmask(ji,jj,jk)
223                        zssh1 = tmask(ji,jj,jk) / zfse3ta(ji,jj,jk)
224                        tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1 * tmask(ji,jj,jk)
225                        sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1 * tmask(ji,jj,jk)
226                     END DO
227                  END DO
228               ELSE                                                     ! Fixed levels
229                 DO jj = 1, jpj
230                     DO ji = 1, jpi
231                        tb(ji,jj,jk) = tn(ji,jj,jk)
232                        sb(ji,jj,jk) = sn(ji,jj,jk)
233                        tn(ji,jj,jk) = ta(ji,jj,jk)
234                        sn(ji,jj,jk) = sa(ji,jj,jk)
235                     END DO
236                  END DO
237               ENDIF
238               IF( l_trdtra ) THEN
239                  ztrdt(:,:,jk) = 0.e0
240                  ztrds(:,:,jk) = 0.e0
241               END IF
242            ELSE
243               IF( l_trdtra ) THEN
244                  DO jj = 1, jpj
245                     DO ji = 1, jpi
246                        ztrdt(ji,jj,jk) = atfp * ( tb(ji,jj,jk) - 2*tn(ji,jj,jk) + ta(ji,jj,jk) )
247                        ztrds(ji,jj,jk) = atfp * ( sb(ji,jj,jk) - 2*sn(ji,jj,jk) + sa(ji,jj,jk) )
248                     END DO
249                  END DO
250               END IF
251               IF( (lk_vvl .AND. ln_zdfexp) ) THEN                      ! Varying levels
252                  DO jj = 1, jpj
253                     DO ji = 1, jpi
254                        zssh1 = tmask(ji,jj,jk) / zfse3tn(ji,jj,jk)
255                        tb(ji,jj,jk) = ( atfp  * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) &
256                          &            + atfp1 * tn(ji,jj,jk) ) * zssh1
257                        sb(ji,jj,jk) = ( atfp  * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) &
258                          &            + atfp1 * sn(ji,jj,jk) ) * zssh1
259                        zssh1 = tmask(ji,jj,1) / zfse3ta(ji,jj,jk)
260                        tn(ji,jj,jk) = ta(ji,jj,jk) * zssh1
261                        sn(ji,jj,jk) = sa(ji,jj,jk) * zssh1
262                     END DO
263                  END DO
264               ELSE                                                     ! Fixed levels or first varying level
265                  DO jj = 1, jpj
266                     DO ji = 1, jpi
267                        tb(ji,jj,jk) = atfp  * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk)
268                        sb(ji,jj,jk) = atfp  * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk)
269                        tn(ji,jj,jk) = ta(ji,jj,jk)
270                        sn(ji,jj,jk) = sa(ji,jj,jk)
271                     END DO
272                  END DO
273               ENDIF
274            ENDIF
275         ENDIF
276         !                                             ! ===============
277      END DO                                           !   End of slab
278      !                                                ! ===============
279
280      IF( l_trdtra )   THEN      ! Take the Asselin trend into account
281         ztrdt(:,:,:) = ztrdt(:,:,:) / ( 2.*rdt )
282         ztrds(:,:,:) = ztrds(:,:,:) / ( 2.*rdt )         
283         CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt )
284      END IF
285
286      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt  - Tn: ', mask1=tmask,   &
287         &                       tab3d_2=sn, clinfo2=       ' Sn: ', mask2=tmask )
288#if defined key_agrif
289      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Tra( kt )
290#endif     
291      !
292   END SUBROUTINE tra_nxt
293
294   !!======================================================================
295END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.