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/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 12555

Last change on this file since 12555 was 12555, checked in by charris, 4 years ago

Changes from GO6 package branch (GMED ticket 450):

svn merge -r 11035:11101 svn+ssh://charris@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_package

File size: 19.1 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         IF(lwp .AND. lflush) CALL flush(numout)
142         !
143         IF( .NOT. ALLOCATED( upsmsk ) )  THEN
144             ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
145             IF( ierr /= 0 )   CALL ctl_stop('STOP', 'tra_adv_cen2: unable to allocate array')
146         ENDIF
147
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      ENDIF
160      !
161      ! Upstream / centered scheme indicator
162      ! ------------------------------------
163!!gm  not strickly exact : the freezing point should be computed at each ocean levels...
164!!gm  not a big deal since cen2 is no more used in global ice-ocean simulations
165!!ch  changes for ice shelf to retain standard behaviour elsewhere, even if not optimal
166      DO jj = 1, jpj 
167         DO ji = 1, jpi 
168            ikt = mikt(ji,jj) 
169            IF (ikt > 1 ) THEN
170               zpres(ji,jj) = grav * rau0 * fsdept(ji,jj,ikt) * 1.e-04 
171            ELSE
172               zpres(ji,jj) = 0.0 
173            ENDIF 
174         END DO
175      END DO
176      CALL eos_fzp( tsn(:,:,1,jp_sal), zfzp(:,:), zpres(:,:) )
177      DO jk = 1, jpk
178         DO jj = 1, jpj
179            DO ji = 1, jpi
180               !                                        ! below ice covered area (if tn < "freezing"+0.1 )
181               IF( tsn(ji,jj,jk,jp_tem) <= zfzp(ji,jj) + 0.1 ) THEN   ;   zice = 1._wp
182               ELSE                                                   ;   zice = 0._wp
183               ENDIF
184               zind(ji,jj,jk) = MAX (   &
185                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
186                  upsmsk(ji,jj)               ,      &  ! some of some straits
187                  zice                               &  ! below ice covered area (if tn < "freezing"+0.1 )
188                  &                  ) * tmask(ji,jj,jk)
189            END DO
190         END DO
191      END DO
192
193      DO jn = 1, kjpt
194         !
195         ! I. Horizontal advection
196         !    ====================
197         !
198         DO jk = 1, jpkm1
199            !                        ! Second order centered tracer flux at u- and v-points
200            DO jj = 1, jpjm1
201               !
202               DO ji = 1, fs_jpim1   ! vector opt.
203                  ! upstream indicator
204                  zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
205                  zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
206                  !
207                  ! upstream scheme
208                  zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
209                  zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
210                  zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
211                  zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
212                  zupsut = zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj  ,jk,jn)
213                  zupsvt = zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji  ,jj+1,jk,jn)
214                  ! centered scheme
215                  zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj  ,jk,jn) )
216                  zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji  ,jj+1,jk,jn) )
217                  ! mixed centered / upstream scheme
218                  zwx(ji,jj,jk) = 0.5 * ( zcofi * zupsut + (1.-zcofi) * zcenut )
219                  zwy(ji,jj,jk) = 0.5 * ( zcofj * zupsvt + (1.-zcofj) * zcenvt )
220               END DO
221            END DO
222         END DO
223
224         ! II. Vertical advection
225         !     ==================
226         !
227         !                                                ! Vertical advective fluxes
228         zwz(:,:,jpk) = 0.e0                                   ! Bottom  value : flux set to zero
229         !                                                     ! Surface value :
230         IF( lk_vvl ) THEN   ;   zwz(:,:, 1 ) = 0.e0                         ! volume variable
231         ELSE
232            DO jj = 1, jpj   ! vector opt.
233               DO ji = 1, jpi   ! vector opt.
234                  ikt = mikt(ji,jj)               
235                  zwz(ji,jj,ikt ) = pwn(ji,jj,ikt) * ptn(ji,jj,ikt,jn)   ! linear free surface
236                  zwz(ji,jj,1:ikt-1) = 0.e0
237               END DO
238            END DO
239         ENDIF
240         !
241         DO jk = 2, jpk              ! Second order centered tracer flux at w-point
242            DO jj = 2, jpjm1
243               DO ji = fs_2, fs_jpim1   ! vector opt.
244                  ! upstream indicator
245                  zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) ) 
246                  ! mixed centered / upstream scheme
247                  zfp_w = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
248                  zfm_w = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
249                  zupst = zfp_w * ptb(ji,jj,jk,jn) + zfm_w * ptb(ji,jj,jk-1,jn)
250                  ! centered scheme
251                  zcent = pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )
252                  ! mixed centered / upstream scheme
253                  zwz(ji,jj,jk) = 0.5 * ( zcofk * zupst + (1.-zcofk) * zcent )
254               END DO
255            END DO
256         END DO
257
258         ! II. Divergence of advective fluxes
259         ! ----------------------------------
260         DO jk = 1, jpkm1
261            DO jj = 2, jpjm1
262               DO ji = fs_2, fs_jpim1   ! vector opt.
263                  zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) *  fse3t(ji,jj,jk) )
264                  ! advective trends
265                  ztra = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk  )   &
266                  &                + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
267                  &                + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1)  )
268                  ! advective trends added to the general tracer trends
269                  pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
270               END DO
271            END DO
272         END DO
273
274         !                                 ! trend diagnostics
275         IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR.    &
276            &( cdtype == 'TRC' .AND. l_trdtrc ) )   THEN
277            CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
278            CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
279            CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
280         END IF
281         !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
282         IF( cdtype == 'TRA' .AND. ln_diaptr ) CALL dia_ptr_ohst_components( jn, 'adv', zwy(:,:,:) )
283         !
284      END DO
285
286      ! ---------------------------  required in restart file to ensure restartability)
287      ! avmb, avtb will be read in zdfini in restart case as they are used in zdftke, kpp etc...
288      IF( lrst_oce .AND. cdtype == 'TRA' ) THEN
289         IF(nn_timing == 2)  CALL timing_start('iom_rstput')
290         CALL iom_rstput( kt, nitrst, numrow, 'avmb', avmb )
291         CALL iom_rstput( kt, nitrst, numrow, 'avtb', avtb )
292         IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
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.