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/2014/dev_MERGE_2014/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2014/dev_MERGE_2014/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 4944

Last change on this file since 4944 was 4666, checked in by mathiot, 10 years ago

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