source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 7 years ago

Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.

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