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

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

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