source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

File size: 19.0 KB
Line 
1MODULE traadv_cen2
2   !!======================================================================
3   !!                     ***  MODULE  traadv_cen2  ***
4   !! Ocean  tracers:  horizontal & vertical advective trend
5   !!======================================================================
6   !! History :  OPA  ! 2001-08  (G. Madec, E. Durand) v8.2 trahad+trazad=traadv
7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module
8   !!             -   ! 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 trd_oce         ! trends: ocean variables
24   USE trdtra          ! trends manager: tracers
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 trc_oce         ! share passive tracers/Ocean variables
32   USE lib_mpp         ! MPP library
33   USE wrk_nemo        ! Memory Allocation
34   USE timing          ! Timing
35   USE phycst
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   tra_adv_cen2   ! routine called by traadv.F90
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, kit000, 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 (l_trdtra=T or l_trctra=T): the trend is
86      !!      saved for diagnostics. The trends saved is expressed as
87      !!      Uh.gradh(T), i.e.  save trend = ztra + ptn divn
88      !!
89      !!         Part II : vertical advection
90      !!      For temperature (idem for salinity) the advective trend is com-
91      !!      puted as follows :
92      !!            ztra = 1/e3t dk+1[ zwz ]
93      !!      where the vertical advective flux, zwz, is given by :
94      !!            zwz = zcofk * zupst + (1-zcofk) * zcent
95      !!      with
96      !!        zupsv = upstream flux = wn * (ptb(k) or ptb(k-1) ) [wn>0 or <0]
97      !!        zcenu = centered flux = wn * mk(tn)
98      !!         The surface boundary condition is :
99      !!      variable volume (lk_vvl = T) : zero advective flux
100      !!      lin. free-surf  (lk_vvl = F) : wn(:,:,1) * ptn(:,:,1)
101      !!         Add this trend now to the general trend of tracer (ta,sa):
102      !!             pta = pta + ztra
103      !!         Trend diagnostic (l_trdtra=T or l_trctra=T): the trend is
104      !!      saved for diagnostics. The trends saved is expressed as :
105      !!             save trend =  w.gradz(T) = ztra - ptn divn.
106      !!
107      !! ** Action :  - update pta  with the now advective tracer trends
108      !!              - save trends if needed
109      !!----------------------------------------------------------------------
110      INTEGER                              , INTENT(in   ) ::   kt              ! ocean time-step index
111      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
112      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype          ! =TRA or TRC (tracer indicator)
113      INTEGER                              , INTENT(in   ) ::   kjpt            ! number of tracers
114      REAL(wp), DIMENSION(jpi,jpj,jpk     ), INTENT(in   ) ::   pun, pvn, pwn   ! 3 ocean velocity components
115      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb, ptn        ! before and now tracer fields
116      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta             ! tracer trend
117      !
118      INTEGER  ::   ji, jj, jk, jn, ikt   ! dummy loop indices
119      INTEGER  ::   ierr                 ! local integer
120      REAL(wp) ::   zbtr, ztra                            ! local scalars
121      REAL(wp) ::   zfp_ui, zfp_vj, zfp_w, zcofi          !   -      -
122      REAL(wp) ::   zfm_ui, zfm_vj, zfm_w, zcofj, zcofk   !   -      -
123      REAL(wp) ::   zupsut, zcenut, zupst                 !   -      -
124      REAL(wp) ::   zupsvt, zcenvt, zcent, zice           !   -      -
125      REAL(wp), POINTER, DIMENSION(:,:)   :: zfzp, zpres   ! 2D workspace
126      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy     ! 3D     -
127      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zind    !  -     -
128      !!----------------------------------------------------------------------
129      !
130      IF( nn_timing == 1 )  CALL timing_start('tra_adv_cen2')
131      !
132      CALL wrk_alloc( jpi, jpj, zpres, zfzp )
133      CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy, zwz, zind )
134      !
135
136      IF( kt == kit000 )  THEN
137         IF(lwp) WRITE(numout,*)
138         IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme on ', cdtype
139         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
140         IF(lwp) WRITE(numout,*)
141         !
142         IF( .NOT. ALLOCATED( upsmsk ) )  THEN
143             ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
144             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_cen2: unable to allocate array')
145         ENDIF
146
147         !
148         upsmsk(:,:) = 0._wp                             ! not upstream by default
149         !
150         IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits
151         !                                               ! and in closed seas (orca2 and orca4 only)
152         IF( jp_cfg == 2 .AND. .NOT. ln_rstart ) THEN    ! Increase the background in the surface layers
153            avmb(1) = 10.  * avmb(1)      ;      avtb(1) = 10.  * avtb(1)
154            avmb(2) = 10.  * avmb(2)      ;      avtb(2) = 10.  * avtb(2)
155            avmb(3) =  5.  * avmb(3)      ;      avtb(3) =  5.  * avtb(3)
156            avmb(4) =  2.5 * avmb(4)      ;      avtb(4) =  2.5 * avtb(4)
157         ENDIF
158      ENDIF
159      !
160      ! Upstream / centered scheme indicator
161      ! ------------------------------------
162!!gm  not strickly exact : the freezing point should be computed at each ocean levels...
163!!gm  not a big deal since cen2 is no more used in global ice-ocean simulations
164!!ch  changes for ice shelf to retain standard behaviour elsewhere, even if not optimal
165      DO jj = 1, jpj 
166         DO ji = 1, jpi 
167            ikt = mikt(ji,jj) 
168            IF (ikt > 1 ) THEN
169               zpres(ji,jj) = grav * rau0 * fsdept(ji,jj,ikt) * 1.e-04 
170            ELSE
171               zpres(ji,jj) = 0.0 
172            ENDIF 
173         END DO
174      END DO
175      CALL eos_fzp( tsn(:,:,1,jp_sal), zfzp(:,:), zpres(:,:) )
176      DO jk = 1, jpk
177         DO jj = 1, jpj
178            DO ji = 1, jpi
179               !                                        ! below ice covered area (if tn < "freezing"+0.1 )
180               IF( tsn(ji,jj,jk,jp_tem) <= zfzp(ji,jj) + 0.1 ) THEN   ;   zice = 1._wp
181               ELSE                                                   ;   zice = 0._wp
182               ENDIF
183               zind(ji,jj,jk) = MAX (   &
184                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
185                  upsmsk(ji,jj)               ,      &  ! some of some straits
186                  zice                               &  ! below ice covered area (if tn < "freezing"+0.1 )
187                  &                  ) * tmask(ji,jj,jk)
188            END DO
189         END DO
190      END DO
191
192      DO jn = 1, kjpt
193         !
194         ! I. Horizontal advection
195         !    ====================
196         !
197         DO jk = 1, jpkm1
198            !                        ! Second order centered tracer flux at u- and v-points
199            DO jj = 1, jpjm1
200               !
201               DO ji = 1, fs_jpim1   ! vector opt.
202                  ! upstream indicator
203                  zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
204                  zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
205                  !
206                  ! upstream scheme
207                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
208                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
209                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
210                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
211                  zupsut = zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn)
212                  zupsvt = zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn)
213                  ! centered scheme
214                  zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
215                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
216                  ! mixed centered / upstream scheme
217                  zwx(ji,jj,jk) = 0.5 * ( zcofi * zupsut + (1.-zcofi) * zcenut )
218                  zwy(ji,jj,jk) = 0.5 * ( zcofj * zupsvt + (1.-zcofj) * zcenvt )
219               END DO
220            END DO
221         END DO
222
223         ! II. Vertical advection
224         !     ==================
225         !
226         !                                                ! Vertical advective fluxes
227         zwz(:,:,jpk) = 0.e0                                   ! Bottom  value : flux set to zero
228         !                                                     ! Surface value :
229         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                         ! volume variable
230         ELSE
231            DO jj = 1, jpj   ! vector opt.
232               DO ji = 1, jpi   ! vector opt.
233                  ikt = mikt(ji,jj)               
234                  zwz(ji,jj,ikt ) = pwn(ji,jj,ikt) * ptn(ji,jj,ikt,jn)   ! linear free surface
235                  zwz(ji,jj,1:ikt-1) = 0.e0
236               END DO
237            END DO
238         ENDIF
239         !
240         DO jk = 2, jpk              ! Second order centered tracer flux at w-point
241            DO jj = 2, jpjm1
242               DO ji = fs_2, fs_jpim1   ! vector opt.
243                  ! upstream indicator
244                  zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) 
245                  ! mixed centered / upstream scheme
246                  zfp_w = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
247                  zfm_w = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
248                  zupst = zfp_w * ptb(ji,jj,jk,jn) + zfm_w * ptb(ji,jj,jk-1,jn)
249                  ! centered scheme
250                  zcent = pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )
251                  ! mixed centered / upstream scheme
252                  zwz(ji,jj,jk) = 0.5 * ( zcofk * zupst + (1.-zcofk) * zcent )
253               END DO
254            END DO
255         END DO
256
257         ! II. Divergence of advective fluxes
258         ! ----------------------------------
259         DO jk = 1, jpkm1
260            DO jj = 2, jpjm1
261               DO ji = fs_2, fs_jpim1   ! vector opt.
262                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) *  fse3t(ji,jj,jk) )
263                  ! advective trends
264                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
265                  &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
266                  &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  )
267                  ! advective trends added to the general tracer trends
268                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
269               END DO
270            END DO
271         END DO
272
273         !                                 ! trend diagnostics
274         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.    &
275            &( cdtype == 'TRC' .AND. l_trdtrc ) )   THEN
276            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
277            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
278            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
279         END IF
280         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
281         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 
282           IF( jn == jp_tem )   htr_adv(:) = ptr_sj( zwy(:,:,:) )
283           IF( jn == jp_sal )   str_adv(:) = ptr_sj( zwy(:,:,:) )
284         ENDIF
285         !
286      END DO
287
288      ! ---------------------------  required in restart file to ensure restartability)
289      ! avmb, avtb will be read in zdfini in restart case as they are used in zdftke, kpp etc...
290      IF( lrst_oce .AND. cdtype == 'TRA' ) THEN
291         CALL iom_rstput( kt, nitrst, numrow, 'avmb', avmb )
292         CALL iom_rstput( kt, nitrst, numrow, 'avtb', avtb )
293      ENDIF
294      !
295      CALL wrk_dealloc( jpi, jpj, zpres, zfzp )
296      CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy, zwz, zind )
297      !
298      IF( nn_timing == 1 )  CALL timing_stop('tra_adv_cen2')
299      !
300   END SUBROUTINE tra_adv_cen2
301   
302   
303   SUBROUTINE ups_orca_set
304      !!----------------------------------------------------------------------
305      !!                  ***  ROUTINE ups_orca_set  ***
306      !!       
307      !! ** Purpose :   add a portion of upstream scheme in area where the
308      !!                centered scheme generates too strong overshoot
309      !!
310      !! ** Method  :   orca (R4 and R2) confiiguration setting. Set upsmsk
311      !!                array to nozero value in some straith.
312      !!
313      !! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca
314      !!----------------------------------------------------------------------
315      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integers
316      !!----------------------------------------------------------------------
317      !
318      IF( nn_timing == 1 )  CALL timing_start('ups_orca_set')
319      !
320      ! mixed upstream/centered scheme near river mouths
321      ! ------------------------------------------------
322      SELECT CASE ( jp_cfg )
323      !                                        ! =======================
324      CASE ( 4 )                               !  ORCA_R4 configuration
325         !                                     ! =======================
326         !                                          ! Gibraltar Strait
327         ii0 =  70   ;   ii1 =  71
328         ij0 =  52   ;   ij1 =  53   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
329         !
330         !                                     ! =======================
331      CASE ( 2 )                               !  ORCA_R2 configuration
332         !                                     ! =======================
333         !                                          ! Gibraltar Strait
334         ij0 = 102   ;   ij1 = 102
335         ii0 = 138   ;   ii1 = 138   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20
336         ii0 = 139   ;   ii1 = 139   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
337         ii0 = 140   ;   ii1 = 140   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
338         ij0 = 101   ;   ij1 = 102
339         ii0 = 141   ;   ii1 = 141   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
340         !                                          ! Bab el Mandeb Strait
341         ij0 =  87   ;   ij1 =  88
342         ii0 = 164   ;   ii1 = 164   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10
343         ij0 =  88   ;   ij1 =  88
344         ii0 = 163   ;   ii1 = 163   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
345         ii0 = 162   ;   ii1 = 162   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
346         ii0 = 160   ;   ii1 = 161   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
347         ij0 =  89   ;   ij1 =  89
348         ii0 = 158   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
349         ij0 =  90   ;   ij1 =  90
350         ii0 = 160   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
351         !                                          ! Sound Strait
352         ij0 = 116   ;   ij1 = 116
353         ii0 = 144   ;   ii1 = 144   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
354         ii0 = 145   ;   ii1 = 147   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
355         ii0 = 148   ;   ii1 = 148   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
356         !
357      END SELECT 
358     
359      ! mixed upstream/centered scheme over closed seas
360      ! -----------------------------------------------
361      CALL clo_ups( upsmsk(:,:) )
362      !
363      IF( nn_timing == 1 )  CALL timing_stop('ups_orca_set')
364      !
365   END SUBROUTINE ups_orca_set
366
367   !!======================================================================
368END MODULE traadv_cen2
Note: See TracBrowser for help on using the repository browser.