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

source: tags/nemo_v2_2/NEMO/OPA_SRC/TRA/tranxt.F90 @ 8115

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

nemo_v2_update_001 : CT : - add non linear free surface (variable volume) with new cpp key key_vvl
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.9 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       ! temporary scalars
90      REAL(wp) ::   zfact        ! temporary scalar
91      !! Variable volume
92      REAL(wp) ::   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         DO jk = 1, jpkm1
103            zfse3tb(:,:,jk)  = fsve3t(:,:,jk) * ( 1 + sshb(:,:) * mut(:,:,jk) )
104            zfse3ta(:,:,jk)  = fsve3t(:,:,jk) * ( 1 + ssha(:,:) * mut(:,:,jk) )
105         END DO
106
107         ! Asselin filtered scale factor at now time step
108         ! ----------------------------------------------
109         IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN
110            zfse3tn(:,:,:) = fse3t(:,:,:)
111         ELSE
112            DO jk = 1, jpkm1
113               DO jj = 1, jpj
114                  DO ji = 1, jpi
115                     zssh = atfp * ( sshb(ji,jj) + ssha(ji,jj) ) + atfp1 * sshn(ji,jj)
116                     zfse3tn(ji,jj,jk) = fsve3t(ji,jj,jk) * ( 1 + zssh * mut(ji,jj,jk) )
117                  END DO
118               END DO
119            END DO
120         ENDIF
121
122         ! Thickness weighting
123         ! -------------------
124         ta(:,:,1:jpkm1) = ta(:,:,1:jpkm1) * fse3t (:,:,1:jpkm1)
125         sa(:,:,1:jpkm1) = sa(:,:,1:jpkm1) * fse3t (:,:,1:jpkm1)
126
127         tn(:,:,1:jpkm1) = tn(:,:,1:jpkm1) * fse3t (:,:,1:jpkm1)
128         sn(:,:,1:jpkm1) = sn(:,:,1:jpkm1) * fse3t (:,:,1:jpkm1)
129
130         tb(:,:,1:jpkm1) = tb(:,:,1:jpkm1) * zfse3tb(:,:,1:jpkm1)
131         sb(:,:,1:jpkm1) = sb(:,:,1:jpkm1) * zfse3tb(:,:,1:jpkm1)
132
133      ENDIF
134
135      IF( l_trdtra ) THEN
136         ztrdt(:,:,jpk) = 0.e0
137         ztrds(:,:,jpk) = 0.e0
138      ENDIF
139
140      ! 0. Lateral boundary conditions on ( ta, sa )   (T-point, unchanged sign)
141      ! ---------------------------------============
142      CALL lbc_lnk( ta, 'T', 1. )   
143      CALL lbc_lnk( sa, 'T', 1. )
144
145      !                                                ! ===============
146      DO jk = 1, jpkm1                                 ! Horizontal slab
147         !                                             ! ===============
148
149         ! 1. Leap-frog scheme (only in explicit case, otherwise the
150         ! -------------------  time stepping is already done in trazdf)
151         IF( ln_zdfexp ) THEN
152            zfact = 2. * rdttra(jk)
153            IF( neuler == 0 .AND. kt == nit000 )   zfact = rdttra(jk)
154            ta(:,:,jk) = ( tb(:,:,jk) + zfact * ta(:,:,jk) ) * tmask(:,:,jk)
155            sa(:,:,jk) = ( sb(:,:,jk) + zfact * sa(:,:,jk) ) * tmask(:,:,jk)
156            IF(l_trdtra)   CALL ctl_stop( 'tranxt: Asselin ML trend not yet accounted for.' )
157         ENDIF
158
159#if defined key_obc
160         !                                             ! ===============
161      END DO                                           !   End of slab
162      !                                                ! ===============
163      ! Update tracers on open boundaries.
164      CALL obc_tra( kt )
165      !                                                ! ===============
166      DO jk = 1, jpkm1                                 ! Horizontal slab
167         !                                             ! ===============
168#endif
169#if defined key_agrif
170         !                                             ! ===============
171      END DO                                           !   End of slab
172      !                                                ! ===============
173      ! Update tracers on open boundaries.
174      CALL Agrif_tra( kt )
175      !                                                ! ===============
176      DO jk = 1, jpkm1                                 ! Horizontal slab
177         !                                             ! ===============
178#endif
179         ! 2. Time filter and swap of arrays
180         ! ---------------------------------
181
182         IF( ln_dynhpg_imp ) THEN                       ! semi-implicite hpg
183            IF( neuler == 0 .AND. kt == nit000 ) THEN
184               DO jj = 1, jpj
185                  DO ji = 1, jpi
186                     zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25
187                     zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25
188                     tb(ji,jj,jk) = tn(ji,jj,jk)
189                     sb(ji,jj,jk) = sn(ji,jj,jk)
190                     tn(ji,jj,jk) = ta(ji,jj,jk)
191                     sn(ji,jj,jk) = sa(ji,jj,jk)
192                     ta(ji,jj,jk) = zt
193                     sa(ji,jj,jk) = zs
194                  END DO
195               END DO
196               IF( l_trdtra ) THEN
197                  ztrdt(:,:,jk) = 0.e0
198                  ztrds(:,:,jk) = 0.e0
199               END IF
200            ELSE
201               DO jj = 1, jpj
202                  DO ji = 1, jpi
203                     zt = ( ta(ji,jj,jk) + 2. * tn(ji,jj,jk) + tb(ji,jj,jk) ) * 0.25
204                     zs = ( sa(ji,jj,jk) + 2. * sn(ji,jj,jk) + sb(ji,jj,jk) ) * 0.25
205                     tb(ji,jj,jk) = atfp  * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk)
206                     sb(ji,jj,jk) = atfp  * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk)
207                     IF( l_trdtra ) THEN ! ChD ceci est a optimiser, mais ca marche
208                        ztrdt(ji,jj,jk) = tb(ji,jj,jk) - tn(ji,jj,jk)
209                        ztrds(ji,jj,jk) = sb(ji,jj,jk) - sn(ji,jj,jk)
210                     END IF
211                     tn(ji,jj,jk) = ta(ji,jj,jk)
212                     sn(ji,jj,jk) = sa(ji,jj,jk)
213                     ta(ji,jj,jk) = zt
214                     sa(ji,jj,jk) = zs
215                  END DO
216               END DO
217            ENDIF
218         ELSE                                          ! Default case
219            IF( neuler == 0 .AND. kt == nit000 ) THEN
220               IF( (lk_vvl .AND. ln_zdfexp) ) THEN                      ! Varying levels
221                  DO jj = 1, jpj
222                     DO ji = 1, jpi
223                        zssh = tmask(ji,jj,jk) / fse3t(ji,jj,jk)
224                        tb(ji,jj,jk) = tn(ji,jj,jk) * zssh * tmask(ji,jj,jk)
225                        sb(ji,jj,jk) = sn(ji,jj,jk) * zssh * tmask(ji,jj,jk)
226                        zssh = tmask(ji,jj,jk) / zfse3ta(ji,jj,jk)
227                        tn(ji,jj,jk) = ta(ji,jj,jk) * zssh * tmask(ji,jj,jk)
228                        sn(ji,jj,jk) = sa(ji,jj,jk) * zssh * tmask(ji,jj,jk)
229                     END DO
230                  END DO
231               ELSE                                                     ! Fixed levels
232                 DO jj = 1, jpj
233                     DO ji = 1, jpi
234                        tb(ji,jj,jk) = tn(ji,jj,jk)
235                        sb(ji,jj,jk) = sn(ji,jj,jk)
236                        tn(ji,jj,jk) = ta(ji,jj,jk)
237                        sn(ji,jj,jk) = sa(ji,jj,jk)
238                     END DO
239                  END DO
240               ENDIF
241               IF( l_trdtra ) THEN
242                  ztrdt(:,:,jk) = 0.e0
243                  ztrds(:,:,jk) = 0.e0
244               END IF
245            ELSE
246               IF( l_trdtra ) THEN
247                  DO jj = 1, jpj
248                     DO ji = 1, jpi
249                        ztrdt(ji,jj,jk) = atfp * ( tb(ji,jj,jk) - 2*tn(ji,jj,jk) + ta(ji,jj,jk) )
250                        ztrds(ji,jj,jk) = atfp * ( sb(ji,jj,jk) - 2*sn(ji,jj,jk) + sa(ji,jj,jk) )
251                     END DO
252                  END DO
253               END IF
254               IF( (lk_vvl .AND. ln_zdfexp) ) THEN                      ! Varying levels
255                  DO jj = 1, jpj
256                     DO ji = 1, jpi
257                        zssh = tmask(ji,jj,jk) / zfse3tn(ji,jj,jk)
258                        tb(ji,jj,jk) = ( atfp  * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) &
259                          &            + atfp1 * tn(ji,jj,jk) ) * zssh
260                        sb(ji,jj,jk) = ( atfp  * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) &
261                          &            + atfp1 * sn(ji,jj,jk) ) * zssh
262                        zssh = tmask(ji,jj,1) / zfse3ta(ji,jj,jk)
263                        tn(ji,jj,jk) = ta(ji,jj,jk) * zssh
264                        sn(ji,jj,jk) = sa(ji,jj,jk) * zssh
265                     END DO
266                  END DO
267               ELSE                                                     ! Fixed levels or first varying level
268                  DO jj = 1, jpj
269                     DO ji = 1, jpi
270                        tb(ji,jj,jk) = atfp  * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk)
271                        sb(ji,jj,jk) = atfp  * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk)
272                        tn(ji,jj,jk) = ta(ji,jj,jk)
273                        sn(ji,jj,jk) = sa(ji,jj,jk)
274                     END DO
275                  END DO
276               ENDIF
277            ENDIF
278         ENDIF
279         !                                             ! ===============
280      END DO                                           !   End of slab
281      !                                                ! ===============
282
283      IF( l_trdtra )   THEN      ! Take the Asselin trend into account
284         ztrdt(:,:,:) = ztrdt(:,:,:) / ( 2.*rdt )
285         ztrds(:,:,:) = ztrds(:,:,:) / ( 2.*rdt )         
286         CALL trd_mod( ztrdt, ztrds, jptra_trd_atf, 'TRA', kt )
287      END IF
288
289      IF(ln_ctl)   CALL prt_ctl( tab3d_1=tn, clinfo1=' nxt  - Tn: ', mask1=tmask,   &
290         &                       tab3d_2=sn, clinfo2=       ' Sn: ', mask2=tmask )
291#if defined key_agrif
292      IF (.NOT.Agrif_Root())    CALL Agrif_Update_Tra( kt )
293#endif     
294      !
295   END SUBROUTINE tra_nxt
296
297   !!======================================================================
298END MODULE tranxt
Note: See TracBrowser for help on using the repository browser.