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

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

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