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.
trcadv.F90 in branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/TRP – NEMO

source: branches/NERC/dev_r5518_GO6_conserv_check_up/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90 @ 9260

Last change on this file since 9260 was 8442, checked in by frrh, 7 years ago

Commit changes relating to Met Office GMED ticket 340 for the
tidying of MEDUSA related code and debugging statements in the TOP code.

Only code introduced at revision 8434 of branch
​http://fcm3/projects/NEMO.xm/log/branches/NERC/dev_r5518_GO6_split_trcbiomedusa
is included here, all previous revisions of that branch having been dealt with
under GMED ticket 339.
File size: 11.8 KB
Line 
1MODULE trcadv
2   !!==============================================================================
3   !!                       ***  MODULE  trcadv  ***
4   !! Ocean passive tracers:  advection trend
5   !!==============================================================================
6   !! History :  2.0  !  05-11  (G. Madec)  Original code
7   !!            3.0  !  10-06  (C. Ethe)   Adapted to passive tracers
8   !!----------------------------------------------------------------------
9#if defined key_top
10   !!----------------------------------------------------------------------
11   !!   'key_top'                                                TOP models
12   !!----------------------------------------------------------------------
13   !!   trc_adv      : compute ocean tracer advection trend
14   !!   trc_adv_ctl  : control the different options of advection scheme
15   !!----------------------------------------------------------------------
16   USE oce_trc         ! ocean dynamics and active tracers
17   USE trc             ! ocean passive tracers variables
18   USE trcnam_trp      ! passive tracers transport namelist variables
19   USE traadv_cen2     ! 2nd order centered scheme (tra_adv_cen2   routine)
20   USE traadv_tvd      ! TVD      scheme           (tra_adv_tvd    routine)
21   USE traadv_muscl    ! MUSCL    scheme           (tra_adv_muscl  routine)
22   USE traadv_muscl2   ! MUSCL2   scheme           (tra_adv_muscl2 routine)
23   USE traadv_ubs      ! UBS      scheme           (tra_adv_ubs    routine)
24   USE traadv_qck      ! QUICKEST scheme           (tra_adv_qck    routine)
25   USE traadv_eiv      ! eddy induced velocity     (tra_adv_eiv    routine)
26   USE traadv_mle      ! ML eddy induced velocity  (tra_adv_mle    routine)
27   USE ldftra_oce      ! lateral diffusion coefficient on tracers
28   USE trd_oce
29   USE trdtra
30   USE prtctl_trc      ! Print control
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   trc_adv          ! routine called by step module
36   PUBLIC   trc_adv_alloc    ! routine called by nemogcm module
37
38   INTEGER ::   nadv   ! choice of the type of advection scheme
39   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   r2dt  ! vertical profile time-step, = 2 rdttra
40   !                                                    ! except at nitrrc000 (=rdttra) if neuler=0
41
42   !! * Substitutions
43#  include "domzgr_substitute.h90"
44#  include "vectopt_loop_substitute.h90"
45   !!----------------------------------------------------------------------
46   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
47   !! $Id$
48   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
49   !!----------------------------------------------------------------------
50CONTAINS
51
52   INTEGER FUNCTION trc_adv_alloc()
53      !!----------------------------------------------------------------------
54      !!                  ***  ROUTINE trc_adv_alloc  ***
55      !!----------------------------------------------------------------------
56
57      ALLOCATE( r2dt(jpk), STAT=trc_adv_alloc )
58
59      IF( trc_adv_alloc /= 0 ) CALL ctl_warn('trc_adv_alloc : failed to allocate array.')
60
61   END FUNCTION trc_adv_alloc
62
63
64   SUBROUTINE trc_adv( kt )
65      !!----------------------------------------------------------------------
66      !!                  ***  ROUTINE trc_adv  ***
67      !!
68      !! ** Purpose :   compute the ocean tracer advection trend.
69      !!
70      !! ** Method  : - Update the tracer with the advection term following nadv
71      !!----------------------------------------------------------------------
72      !!
73      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
74      !
75      INTEGER ::   jk, jn 
76      CHARACTER (len=22) ::   charout
77      REAL(wp), POINTER, DIMENSION(:,:,:) :: zun, zvn, zwn  ! effective velocity
78      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrtrd
79      !!----------------------------------------------------------------------
80      !
81      IF( nn_timing == 1 )  CALL timing_start('trc_adv')
82      !
83      CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn )
84      !
85
86      IF( kt == nittrc000 )   CALL trc_adv_ctl          ! initialisation & control of options
87
88      IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN     ! at nittrc000
89         r2dt(:) =  rdttrc(:)           ! = rdttrc (use or restarting with Euler time stepping)
90      ELSEIF( kt <= nittrc000 + nn_dttrc ) THEN          ! at nittrc000 or nittrc000+1
91         r2dt(:) = 2. * rdttrc(:)       ! = 2 rdttrc (leapfrog)
92      ENDIF
93      !                                                   ! effective transport
94      DO jk = 1, jpkm1
95         !                                                ! eulerian transport only
96         zun(:,:,jk) = e2u  (:,:) * fse3u(:,:,jk) * un(:,:,jk)
97         zvn(:,:,jk) = e1v  (:,:) * fse3v(:,:,jk) * vn(:,:,jk)
98         zwn(:,:,jk) = e1e2t(:,:)                 * wn(:,:,jk)
99         !
100      END DO
101      !
102      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
103         zun(:,:,:) = zun(:,:,:) + un_td(:,:,:)
104         zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:)
105      ENDIF
106      !
107      zun(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
108      zvn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
109      zwn(:,:,jpk) = 0._wp                                                     ! no transport trough the bottom
110
111      IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif )   &  ! add the eiv transport (if necessary)
112         &              CALL tra_adv_eiv( kt, nittrc000, zun, zvn, zwn, 'TRC' )
113      !
114      IF( ln_mle    )   CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' )    ! add the mle transport (if necessary)
115      !
116      IF( l_trdtrc )  THEN
117         CALL wrk_alloc( jpi, jpj, jpk, jptra, ztrtrd )
118         ztrtrd(:,:,:,:)  = tra(:,:,:,:)
119      ENDIF
120      !
121      SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==!
122      CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, nittrc000, 'TRC',       zun, zvn, zwn, trb, trn, tra, jptra )   !  2nd order centered
123      CASE ( 2 )   ;    CALL tra_adv_tvd   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  TVD
124      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra, ln_trcadv_msc_ups )   !  MUSCL
125      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  MUSCL2
126      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  UBS
127      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  QUICKEST
128      !
129      CASE (-1 )                                      !==  esopa: test all possibility with control print  ==!
130         CALL tra_adv_cen2  ( kt, nittrc000, 'TRC',       zun, zvn, zwn, trb, trn, tra, jptra )         
131         WRITE(charout, FMT="('adv1')")  ; CALL prt_ctl_trc_info(charout)
132                                           CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
133         CALL tra_adv_tvd   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )         
134         WRITE(charout, FMT="('adv2')")  ; CALL prt_ctl_trc_info(charout)
135                                           CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
136         CALL tra_adv_muscl ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra, ln_trcadv_msc_ups  )         
137         WRITE(charout, FMT="('adv3')")  ; CALL prt_ctl_trc_info(charout)
138                                           CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
139         CALL tra_adv_muscl2( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )         
140         WRITE(charout, FMT="('adv4')")  ; CALL prt_ctl_trc_info(charout)
141                                           CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
142         CALL tra_adv_ubs   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )         
143         WRITE(charout, FMT="('adv5')")  ; CALL prt_ctl_trc_info(charout)
144                                           CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
145         CALL tra_adv_qck   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )         
146         WRITE(charout, FMT="('adv6')")  ; CALL prt_ctl_trc_info(charout)
147                                           CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
148         !
149      END SELECT
150      !
151      IF( l_trdtrc )   THEN                      ! save the advective trends for further diagnostics
152        DO jn = 1, jptra
153           ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn)
154           CALL trd_tra( kt, 'TRC', jn, jptra_totad, ztrtrd(:,:,:,jn) )
155        END DO
156        CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrtrd )
157      ENDIF
158
159      !                                              ! print mean trends (used for debugging)
160      IF( ln_ctl )   THEN
161         WRITE(charout, FMT="('adv ')")  ;  CALL prt_ctl_trc_info(charout)
162                                            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
163      END IF
164      !
165      CALL wrk_dealloc( jpi, jpj, jpk, zun, zvn, zwn )
166      !
167      IF( nn_timing == 1 )  CALL timing_stop('trc_adv')
168      !
169   END SUBROUTINE trc_adv
170
171
172   SUBROUTINE trc_adv_ctl
173      !!---------------------------------------------------------------------
174      !!                  ***  ROUTINE trc_adv_ctl  ***
175      !!               
176      !! ** Purpose : Control the consistency between namelist options for
177      !!              passive tracer advection schemes and set nadv
178      !!----------------------------------------------------------------------
179      INTEGER ::   ioptio
180      !!----------------------------------------------------------------------
181
182      ioptio = 0                      ! Parameter control
183      IF( ln_trcadv_cen2   )   ioptio = ioptio + 1
184      IF( ln_trcadv_tvd    )   ioptio = ioptio + 1
185      IF( ln_trcadv_muscl  )   ioptio = ioptio + 1
186      IF( ln_trcadv_muscl2 )   ioptio = ioptio + 1
187      IF( ln_trcadv_ubs    )   ioptio = ioptio + 1
188      IF( ln_trcadv_qck    )   ioptio = ioptio + 1
189      IF( lk_esopa         )   ioptio =          1
190
191      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namtrc_adv' )
192
193      !                              ! Set nadv
194      IF( ln_trcadv_cen2   )   nadv =  1
195      IF( ln_trcadv_tvd    )   nadv =  2
196      IF( ln_trcadv_muscl  )   nadv =  3
197      IF( ln_trcadv_muscl2 )   nadv =  4
198      IF( ln_trcadv_ubs    )   nadv =  5
199      IF( ln_trcadv_qck    )   nadv =  6
200      IF( lk_esopa         )   nadv = -1
201
202      IF(lwp) THEN                   ! Print the choice
203         WRITE(numout,*)
204         IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is used'
205         IF( nadv ==  2 )   WRITE(numout,*) '         TVD       scheme is used'
206         IF( nadv ==  3 )   WRITE(numout,*) '         MUSCL     scheme is used'
207         IF( nadv ==  4 )   WRITE(numout,*) '         MUSCL2    scheme is used'
208         IF( nadv ==  5 )   WRITE(numout,*) '         UBS       scheme is used'
209         IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  scheme is used'
210         IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection scheme'
211      ENDIF
212      !
213   END SUBROUTINE trc_adv_ctl
214   
215#else
216   !!----------------------------------------------------------------------
217   !!   Default option                                         Empty module
218   !!----------------------------------------------------------------------
219CONTAINS
220   SUBROUTINE trc_adv( kt )
221      INTEGER, INTENT(in) :: kt
222      WRITE(*,*) 'trc_adv: You should not have seen this print! error?', kt
223   END SUBROUTINE trc_adv
224#endif
225
226  !!======================================================================
227END MODULE trcadv
Note: See TracBrowser for help on using the repository browser.