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.
traadv_cen2.F90 in branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 2636

Last change on this file since 2636 was 2636, checked in by gm, 13 years ago

dynamic mem: #785 ; move ctl_stop & warn in lib_mpp to avoid a circular dependency + ctl_stop improvment
  • Property svn:keywords set to Id
File size: 18.5 KB
Line 
1MODULE traadv_cen2
2   !!======================================================================
3   !!                     ***  MODULE  traadv_cen2  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend
5   !!======================================================================
6   !! History :  8.2  ! 2001-08  (G. Madec, E. Durand)  trahad+trazad=traadv
7   !!            1.0  ! 2002-06  (G. Madec)  F90: Free form and module
8   !!            9.0  ! 2004-08  (C. Talandier) New trends organization
9   !!             -   ! 2005-11  (V. Garnier) Surface pressure gradient organization
10   !!            2.0  ! 2006-04  (R. Benshila, G. Madec) Step reorganization
11   !!             -   ! 2006-07  (G. madec)  add ups_orca_set routine
12   !!            3.2  ! 2009-07  (G. Madec) add avmb, avtb in restart for cen2 advection
13   !!            3.3  ! 2010-05  (C. Ethe, G. Madec)  merge TRC-TRA + switch from velocity to transport
14   !!----------------------------------------------------------------------
15
16   !!----------------------------------------------------------------------
17   !!   tra_adv_cen2 : update the tracer trend with the advection trends using a 2nd order centered scheme
18   !!   ups_orca_set : allow mixed upstream/centered scheme in specific area (set for orca 2 and 4 only)
19   !!----------------------------------------------------------------------
20   USE oce, ONLY: tsn  ! now ocean temperature and salinity
21   USE dom_oce         ! ocean space and time domain
22   USE eosbn2          ! equation of state
23   USE trdmod_oce      ! tracers trends
24   USE trdtra          ! tracers trends
25   USE closea          ! closed sea
26   USE sbcrnf          ! river runoffs
27   USE in_out_manager  ! I/O manager
28   USE iom             ! IOM library
29   USE diaptr          ! poleward transport diagnostics
30   USE zdf_oce         ! ocean vertical physics
31   USE restart         ! ocean restart
32   USE trc_oce         ! share passive tracers/Ocean variables
33   USE lib_mpp           ! MPP library
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   tra_adv_cen2       ! routine called by step.F90
39   PUBLIC   ups_orca_set       ! routine used by traadv_cen2_jki.F90
40
41   LOGICAL  :: l_trd       ! flag to compute trends
42
43   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits
44   !                                                             !  and in closed seas (orca 2 and 4 configurations)
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
50   !! $Id$
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE tra_adv_cen2( kt, cdtype, pun, pvn, pwn,        &
56      &                                 ptb, ptn, pta, kjpt   ) 
57      !!----------------------------------------------------------------------
58      !!                  ***  ROUTINE tra_adv_cen2  ***
59      !!                 
60      !! ** Purpose :   Compute the now trend due to the advection of tracers
61      !!      and add it to the general trend of passive tracer equations.
62      !!
63      !! ** Method  :   The advection is evaluated by a second order centered
64      !!      scheme using now fields (leap-frog scheme). In specific areas
65      !!      (vicinity of major river mouths, some straits, or where tn is
66      !!      approaching the freezing point) it is mixed with an upstream
67      !!      scheme for stability reasons.
68      !!         Part 0 : compute the upstream / centered flag
69      !!                  (3D array, zind, defined at T-point (0<zind<1))
70      !!         Part I : horizontal advection
71      !!       * centered flux:
72      !!               zcenu = e2u*e3u  un  mi(ptn)
73      !!               zcenv = e1v*e3v  vn  mj(ptn)
74      !!       * upstream flux:
75      !!               zupsu = e2u*e3u  un  (ptb(i) or ptb(i-1) ) [un>0 or <0]
76      !!               zupsv = e1v*e3v  vn  (ptb(j) or ptb(j-1) ) [vn>0 or <0]
77      !!       * mixed upstream / centered horizontal advection scheme
78      !!               zcofi = max(zind(i+1), zind(i))
79      !!               zcofj = max(zind(j+1), zind(j))
80      !!               zwx = zcofi * zupsu + (1-zcofi) * zcenu
81      !!               zwy = zcofj * zupsv + (1-zcofj) * zcenv
82      !!       * horizontal advective trend (divergence of the fluxes)
83      !!               ztra = 1/(e1t*e2t*e3t) { di-1[zwx] + dj-1[zwy] }
84      !!       * Add this trend now to the general trend of tracer (ta,sa):
85      !!               pta = pta + ztra
86      !!       * trend diagnostic ('key_trdtra' defined): the trend is
87      !!      saved for diagnostics. The trends saved is expressed as
88      !!      Uh.gradh(T), i.e.
89      !!                     save trend = ztra + ptn divn
90      !!
91      !!         Part II : vertical advection
92      !!      For temperature (idem for salinity) the advective trend is com-
93      !!      puted as follows :
94      !!            ztra = 1/e3t dk+1[ zwz ]
95      !!      where the vertical advective flux, zwz, is given by :
96      !!            zwz = zcofk * zupst + (1-zcofk) * zcent
97      !!      with
98      !!        zupsv = upstream flux = wn * (ptb(k) or ptb(k-1) ) [wn>0 or <0]
99      !!        zcenu = centered flux = wn * mk(tn)
100      !!         The surface boundary condition is :
101      !!      variable volume (lk_vvl = T) : zero advective flux
102      !!      lin. free-surf  (lk_vvl = F) : wn(:,:,1) * ptn(:,:,1)
103      !!         Add this trend now to the general trend of tracer (ta,sa):
104      !!             pta = pta + ztra
105      !!         Trend diagnostic ('key_trdtra' defined): the trend is
106      !!      saved for diagnostics. The trends saved is expressed as :
107      !!             save trend =  w.gradz(T) = ztra - ptn divn.
108      !!
109      !! ** Action :  - update pta  with the now advective tracer trends
110      !!              - save trends if needed
111      !!----------------------------------------------------------------------
112      USE oce         , zwx => ua   ! use ua as workspace
113      USE oce         , zwy => va   ! use va as workspace
114      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
115      USE wrk_nemo, ONLY: ztfreez => wrk_2d_1
116      USE wrk_nemo, ONLY: zwz => wrk_3d_1, zind => wrk_3d_2
117      !!
118      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
119      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
120      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
121      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
122      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
123      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
124      !!
125      INTEGER  ::   ji, jj, jk, jn                   ! dummy loop indices
126      INTEGER  ::   ierr                             ! local integer
127      REAL(wp) ::   zbtr, ztra                       ! local scalars
128      REAL(wp) ::   zfp_ui, zfp_vj, zfp_w            !   -      -
129      REAL(wp) ::   zfm_ui, zfm_vj, zfm_w            !   -      -
130      REAL(wp) ::   zcofi , zcofj , zcofk            !   -      -
131      REAL(wp) ::   zupsut, zcenut                   !   -      -
132      REAL(wp) ::   zupsvt, zcenvt                   !   -      -
133      REAL(wp) ::   zupst , zcent                    !   -      -
134      REAL(wp) ::   zice                             !   -      -
135      !!----------------------------------------------------------------------
136
137      IF( wrk_in_use(2, 1) .OR. wrk_in_use(3, 1,2) ) THEN
138         CALL ctl_stop('tra_adv_cen2: ERROR: requested workspace arrays unavailable')   ;   RETURN
139      END IF
140
141      IF( kt == nit000 )  THEN
142         IF(lwp) WRITE(numout,*)
143         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme on ', cdtype
144         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
145         IF(lwp) WRITE(numout,*)
146         !
147         ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
148         IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_cen2: unable to allocate array')
149         !
150         upsmsk(:,:) = 0._wp                             ! not upstream by default
151         !
152         IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits
153         !                                               ! and in closed seas (orca2 and orca4 only)
154         IF( jp_cfg == 2 .AND. .NOT. ln_rstart ) THEN    ! Increase the background in the surface layers
155            avmb(1) = 10.  * avmb(1)      ;      avtb(1) = 10.  * avtb(1)
156            avmb(2) = 10.  * avmb(2)      ;      avtb(2) = 10.  * avtb(2)
157            avmb(3) =  5.  * avmb(3)      ;      avtb(3) =  5.  * avtb(3)
158            avmb(4) =  2.5 * avmb(4)      ;      avtb(4) =  2.5 * avtb(4)
159         ENDIF
160         !
161         l_trd = .FALSE.
162         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) )   l_trd = .TRUE.
163      ENDIF
164      !
165      ! Upstream / centered scheme indicator
166      ! ------------------------------------
167!!gm  not strickly exact : the freezing point should be computed at each ocean levels...
168!!gm  not a big deal since cen2 is no more used in global ice-ocean simulations
169      ztfreez(:,:) = tfreez( tsn(:,:,1,jp_sal) )
170      DO jk = 1, jpk
171         DO jj = 1, jpj
172            DO ji = 1, jpi
173               !                                        ! below ice covered area (if tn < "freezing"+0.1 )
174               IF( tsn(ji,jj,jk,jp_tem) <= ztfreez(ji,jj) + 0.1 ) THEN   ;   zice = 1.e0
175               ELSE                                                      ;   zice = 0.e0
176               ENDIF
177               zind(ji,jj,jk) = MAX (   &
178                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
179                  upsmsk(ji,jj)               ,      &  ! some of some straits
180                  zice                               &  ! below ice covered area (if tn < "freezing"+0.1 )
181                  &                  ) * tmask(ji,jj,jk)
182            END DO
183         END DO
184      END DO
185
186      DO jn = 1, kjpt
187         !
188         ! I. Horizontal advection
189         !    ====================
190         !
191         DO jk = 1, jpkm1
192            !                        ! Second order centered tracer flux at u- and v-points
193            DO jj = 1, jpjm1
194               !
195               DO ji = 1, fs_jpim1   ! vector opt.
196                  ! upstream indicator
197                  zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
198                  zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
199                  !
200                  ! upstream scheme
201                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
202                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
203                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
204                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
205                  zupsut = zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn)
206                  zupsvt = zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn)
207                  ! centered scheme
208                  zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
209                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
210                  ! mixed centered / upstream scheme
211                  zwx(ji,jj,jk) = 0.5 * ( zcofi * zupsut + (1.-zcofi) * zcenut )
212                  zwy(ji,jj,jk) = 0.5 * ( zcofj * zupsvt + (1.-zcofj) * zcenvt )
213               END DO
214            END DO
215         END DO
216
217         ! II. Vertical advection
218         !     ==================
219         !
220         !                                                ! Vertical advective fluxes
221         zwz(:,:,jpk) = 0.e0                                   ! Bottom  value : flux set to zero
222         !                                                     ! Surface value :
223         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                         ! volume variable
224         ELSE                ;   zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn)   ! linear free surface
225         ENDIF
226         !
227         DO jk = 2, jpk              ! Second order centered tracer flux at w-point
228            DO jj = 2, jpjm1
229               DO ji = fs_2, fs_jpim1   ! vector opt.
230                  ! upstream indicator
231                  zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) 
232                  ! mixed centered / upstream scheme
233                  zfp_w = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
234                  zfm_w = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
235                  zupst = zfp_w * ptb(ji,jj,jk,jn) + zfm_w * ptb(ji,jj,jk-1,jn)
236                  ! centered scheme
237                  zcent = pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )
238                  ! mixed centered / upstream scheme
239                  zwz(ji,jj,jk) = 0.5 * ( zcofk * zupst + (1.-zcofk) * zcent )
240               END DO
241            END DO
242         END DO
243
244         ! II. Divergence of advective fluxes
245         ! ----------------------------------
246         DO jk = 1, jpkm1
247            DO jj = 2, jpjm1
248               DO ji = fs_2, fs_jpim1   ! vector opt.
249                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) *  fse3t(ji,jj,jk) )
250                  ! advective trends
251                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
252                  &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
253                  &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  )
254                  ! advective trends added to the general tracer trends
255                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
256               END DO
257            END DO
258         END DO
259
260         !                                 ! trend diagnostics (contribution of upstream fluxes)
261         IF( l_trd ) THEN
262            CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
263            CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
264            CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
265         END IF
266         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
267         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 
268           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
269           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
270         ENDIF
271         !
272      ENDDO
273
274      ! ---------------------------  required in restart file to ensure restartability)
275      ! avmb, avtb will be read in zdfini in restart case as they are used in zdftke, kpp etc...
276      IF( lrst_oce .AND. cdtype == 'TRA' ) THEN
277         CALL iom_rstput( kt, nitrst, numrow, 'avmb', avmb )
278         CALL iom_rstput( kt, nitrst, numrow, 'avtb', avtb )
279      ENDIF
280      !
281      IF( wrk_not_released(2, 1) .OR. wrk_not_released(3, 1,2) )THEN
282         CALL ctl_stop('tra_adv_cen2: ERROR: failed to release workspace arrays')
283      END IF
284      !
285   END SUBROUTINE tra_adv_cen2
286   
287   
288   SUBROUTINE ups_orca_set
289      !!----------------------------------------------------------------------
290      !!                  ***  ROUTINE ups_orca_set  ***
291      !!       
292      !! ** Purpose :   add a portion of upstream scheme in area where the
293      !!                centered scheme generates too strong overshoot
294      !!
295      !! ** Method  :   orca (R4 and R2) confiiguration setting. Set upsmsk
296      !!                array to nozero value in some straith.
297      !!
298      !! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca
299      !!----------------------------------------------------------------------
300      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integers
301      !!----------------------------------------------------------------------
302     
303      ! mixed upstream/centered scheme near river mouths
304      ! ------------------------------------------------
305      SELECT CASE ( jp_cfg )
306      !                                        ! =======================
307      CASE ( 4 )                               !  ORCA_R4 configuration
308         !                                     ! =======================
309         !                                          ! Gibraltar Strait
310         ii0 =  70   ;   ii1 =  71
311         ij0 =  52   ;   ij1 =  53   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
312         !
313         !                                     ! =======================
314      CASE ( 2 )                               !  ORCA_R2 configuration
315         !                                     ! =======================
316         !                                          ! Gibraltar Strait
317         ij0 = 102   ;   ij1 = 102
318         ii0 = 138   ;   ii1 = 138   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20
319         ii0 = 139   ;   ii1 = 139   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
320         ii0 = 140   ;   ii1 = 140   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
321         ij0 = 101   ;   ij1 = 102
322         ii0 = 141   ;   ii1 = 141   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
323         !                                          ! Bab el Mandeb Strait
324         ij0 =  87   ;   ij1 =  88
325         ii0 = 164   ;   ii1 = 164   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10
326         ij0 =  88   ;   ij1 =  88
327         ii0 = 163   ;   ii1 = 163   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
328         ii0 = 162   ;   ii1 = 162   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
329         ii0 = 160   ;   ii1 = 161   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
330         ij0 =  89   ;   ij1 =  89
331         ii0 = 158   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
332         ij0 =  90   ;   ij1 =  90
333         ii0 = 160   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
334         !                                          ! Sound Strait
335         ij0 = 116   ;   ij1 = 116
336         ii0 = 144   ;   ii1 = 144   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
337         ii0 = 145   ;   ii1 = 147   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
338         ii0 = 148   ;   ii1 = 148   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
339         !
340      END SELECT 
341     
342      ! mixed upstream/centered scheme over closed seas
343      ! -----------------------------------------------
344      CALL clo_ups( upsmsk(:,:) )
345      !
346   END SUBROUTINE ups_orca_set
347
348   !!======================================================================
349END MODULE traadv_cen2
Note: See TracBrowser for help on using the repository browser.