- Timestamp:
- 2017-12-13T15:58:53+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90
r7753 r9019 7 7 !! 3.0 ! 2010-06 (C. Ethe) Adapted to passive tracers 8 8 !! 3.7 ! 2014-05 (G. Madec, C. Ethe) Add 2nd/4th order cases for CEN and FCT schemes 9 !! 4.0 ! 2017-09 (G. Madec) remove vertical time-splitting option 9 10 !!---------------------------------------------------------------------- 10 11 #if defined key_top … … 17 18 USE oce_trc ! ocean dynamics and active tracers 18 19 USE trc ! ocean passive tracers variables 20 USE sbcwave ! wave module 21 USE sbc_oce ! surface boundary condition: ocean 19 22 USE traadv_cen ! centered scheme (tra_adv_cen routine) 20 23 USE traadv_fct ! FCT scheme (tra_adv_fct routine) … … 23 26 USE traadv_qck ! QUICKEST scheme (tra_adv_qck routine) 24 27 USE traadv_mle ! ML eddy induced velocity (tra_adv_mle routine) 25 USE ldftra ! lateral diffusion coefficient on tracers28 USE ldftra ! lateral diffusion: eddy diffusivity & EIV coeff. 26 29 USE ldfslp ! Lateral diffusion: slopes of neutral surfaces 27 30 ! 28 USE prtctl_trc ! Print control 31 USE prtctl_trc ! control print 32 USE timing ! Timing 29 33 30 34 IMPLICIT NONE 31 35 PRIVATE 32 36 33 PUBLIC trc_adv 34 PUBLIC trc_adv_ini 37 PUBLIC trc_adv ! called by trctrp.F90 38 PUBLIC trc_adv_ini ! called by trcini.F90 35 39 36 40 ! !!* Namelist namtrc_adv * 41 LOGICAL :: ln_trcadv_NONE ! no advection on passive tracers 37 42 LOGICAL :: ln_trcadv_cen ! centered scheme flag 38 43 INTEGER :: nn_cen_h, nn_cen_v ! =2/4 : horizontal and vertical choices of the order of CEN scheme 39 44 LOGICAL :: ln_trcadv_fct ! FCT scheme flag 40 45 INTEGER :: nn_fct_h, nn_fct_v ! =2/4 : horizontal and vertical choices of the order of FCT scheme 41 INTEGER :: nn_fct_zts ! >=1 : 2nd order FCT with vertical sub-timestepping42 46 LOGICAL :: ln_trcadv_mus ! MUSCL scheme flag 43 47 LOGICAL :: ln_mus_ups ! use upstream scheme in vivcinity of river mouths … … 46 50 LOGICAL :: ln_trcadv_qck ! QUICKEST scheme flag 47 51 48 ! ! choices of advection scheme: 52 INTEGER :: nadv ! choice of the type of advection scheme 53 ! ! associated indices: 49 54 INTEGER, PARAMETER :: np_NO_adv = 0 ! no T-S advection 50 55 INTEGER, PARAMETER :: np_CEN = 1 ! 2nd/4th order centered scheme 51 56 INTEGER, PARAMETER :: np_FCT = 2 ! 2nd/4th order Flux Corrected Transport scheme 52 INTEGER, PARAMETER :: np_FCT_zts = 3 ! 2nd order FCT scheme with vertical sub-timestepping 53 INTEGER, PARAMETER :: np_MUS = 4 ! MUSCL scheme 54 INTEGER, PARAMETER :: np_UBS = 5 ! 3rd order Upstream Biased Scheme 55 INTEGER, PARAMETER :: np_QCK = 6 ! QUICK scheme 56 57 INTEGER :: nadv ! chosen advection scheme 58 ! 57 INTEGER, PARAMETER :: np_MUS = 3 ! MUSCL scheme 58 INTEGER, PARAMETER :: np_UBS = 4 ! 3rd order Upstream Biased Scheme 59 INTEGER, PARAMETER :: np_QCK = 5 ! QUICK scheme 60 59 61 !! * Substitutions 60 62 # include "vectopt_loop_substitute.h90" 61 63 !!---------------------------------------------------------------------- 62 !! NEMO/TOP 3.7 , NEMO Consortium (2015)64 !! NEMO/TOP 4.0 , NEMO Consortium (2017) 63 65 !! $Id$ 64 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)66 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 65 67 !!---------------------------------------------------------------------- 66 68 CONTAINS … … 72 74 !! ** Purpose : compute the ocean tracer advection trend. 73 75 !! 74 !! ** Method : - Update the tracerwith the advection term following nadv76 !! ** Method : - Update after tracers (tra) with the advection term following nadv 75 77 !!---------------------------------------------------------------------- 76 78 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 78 80 INTEGER :: jk ! dummy loop index 79 81 CHARACTER (len=22) :: charout 80 REAL(wp), POINTER, DIMENSION(:,:,:) :: zun, zvn, zwn ! effective velocity 81 !!---------------------------------------------------------------------- 82 ! 83 IF( nn_timing == 1 ) CALL timing_start('trc_adv') 84 ! 85 CALL wrk_alloc( jpi,jpj,jpk, zun, zvn, zwn ) 86 ! !== effective transport ==! 82 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zun, zvn, zwn ! effective velocity 83 !!---------------------------------------------------------------------- 84 ! 85 IF( ln_timing ) CALL timing_start('trc_adv') 86 ! 87 ! !== effective transport ==! 87 88 IF( l_offline ) THEN 88 zun(:,:,:) = un(:,:,:) ! effective transport already in un/vn/wn89 zun(:,:,:) = un(:,:,:) ! already in (un,vn,wn) 89 90 zvn(:,:,:) = vn(:,:,:) 90 91 zwn(:,:,:) = wn(:,:,:) 91 ELSE 92 ! 93 DO jk = 1, jpkm1 94 zun(:,:,jk) = e2u (:,:) * e3u_n(:,:,jk) * un(:,:,jk) ! eulerian transport 95 zvn(:,:,jk) = e1v (:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 96 zwn(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) 97 END DO 92 ELSE ! build the effective transport 93 zun(:,:,jpk) = 0._wp 94 zvn(:,:,jpk) = 0._wp 95 zwn(:,:,jpk) = 0._wp 96 IF( ln_wave .AND. ln_sdw ) THEN 97 DO jk = 1, jpkm1 ! eulerian transport + Stokes Drift 98 zun(:,:,jk) = e2u (:,:) * e3u_n(:,:,jk) * ( un(:,:,jk) + usd(:,:,jk) ) 99 zvn(:,:,jk) = e1v (:,:) * e3v_n(:,:,jk) * ( vn(:,:,jk) + vsd(:,:,jk) ) 100 zwn(:,:,jk) = e1e2t(:,:) * ( wn(:,:,jk) + wsd(:,:,jk) ) 101 END DO 102 ELSE 103 DO jk = 1, jpkm1 104 zun(:,:,jk) = e2u (:,:) * e3u_n(:,:,jk) * un(:,:,jk) ! eulerian transport 105 zvn(:,:,jk) = e1v (:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 106 zwn(:,:,jk) = e1e2t(:,:) * wn(:,:,jk) 107 END DO 108 ENDIF 98 109 ! 99 110 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! add z-tilde and/or vvl corrections … … 107 118 IF( ln_mle ) CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' ) ! add the mle transport 108 119 ! 109 zun(:,:,jpk) = 0._wp ! no transport trough the bottom110 zvn(:,:,jpk) = 0._wp111 zwn(:,:,jpk) = 0._wp112 !113 120 ENDIF 114 121 ! 115 122 SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==! 116 123 ! 117 CASE ( np_CEN ) ! Centered : 2nd / 4th order 118 CALL tra_adv_cen ( kt, nittrc000,'TRC', zun, zvn, zwn , trn, tra, jptra, nn_cen_h, nn_cen_v ) 119 CASE ( np_FCT ) ! FCT : 2nd / 4th order 120 CALL tra_adv_fct ( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra, nn_fct_h, nn_fct_v ) 121 CASE ( np_FCT_zts ) ! 2nd order FCT with vertical time-splitting 122 CALL tra_adv_fct_zts( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra , nn_fct_zts ) 123 CASE ( np_MUS ) ! MUSCL 124 CALL tra_adv_mus ( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb, tra, jptra , ln_mus_ups ) 125 CASE ( np_UBS ) ! UBS 126 CALL tra_adv_ubs ( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra , nn_ubs_v ) 127 CASE ( np_QCK ) ! QUICKEST 128 CALL tra_adv_qck ( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra ) 124 CASE ( np_CEN ) ! Centered : 2nd / 4th order 125 CALL tra_adv_cen( kt, nittrc000,'TRC', zun, zvn, zwn , trn, tra, jptra, nn_cen_h, nn_cen_v ) 126 CASE ( np_FCT ) ! FCT : 2nd / 4th order 127 CALL tra_adv_fct( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra, nn_fct_h, nn_fct_v ) 128 CASE ( np_MUS ) ! MUSCL 129 CALL tra_adv_mus( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb, tra, jptra , ln_mus_ups ) 130 CASE ( np_UBS ) ! UBS 131 CALL tra_adv_ubs( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra , nn_ubs_v ) 132 CASE ( np_QCK ) ! QUICKEST 133 CALL tra_adv_qck( kt, nittrc000,'TRC', r2dttrc, zun, zvn, zwn, trb, trn, tra, jptra ) 129 134 ! 130 135 END SELECT 131 136 ! 132 IF( ln_ctl ) THEN !== print mean trends (used for debugging) 133 WRITE(charout, FMT="('adv ')") ; CALL prt_ctl_trc_info(charout) 134 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) 137 IF( ln_ctl ) THEN !== print mean trends (used for debugging) 138 WRITE(charout, FMT="('adv ')") 139 CALL prt_ctl_trc_info(charout) 140 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' ) 135 141 END IF 136 142 ! 137 CALL wrk_dealloc( jpi,jpj,jpk, zun, zvn, zwn ) 138 ! 139 IF( nn_timing == 1 ) CALL timing_stop('trc_adv') 143 IF( ln_timing ) CALL timing_stop('trc_adv') 140 144 ! 141 145 END SUBROUTINE trc_adv … … 146 150 !! *** ROUTINE trc_adv_ini *** 147 151 !! 148 !! ** Purpose : Control the consistency between namelist options for152 !! ** Purpose : Control the consistency between namelist options for 149 153 !! passive tracer advection schemes and set nadv 150 154 !!---------------------------------------------------------------------- … … 152 156 INTEGER :: ios ! Local integer output status for namelist read 153 157 !! 154 NAMELIST/namtrc_adv/ ln_trcadv_cen, nn_cen_h, nn_cen_v, & ! CEN 155 & ln_trcadv_fct, nn_fct_h, nn_fct_v, nn_fct_zts, & ! FCT 156 & ln_trcadv_mus, ln_mus_ups, & ! MUSCL 157 & ln_trcadv_ubs, nn_ubs_v, & ! UBS 158 & ln_trcadv_qck ! QCK 159 !!---------------------------------------------------------------------- 160 ! 161 REWIND( numnat_ref ) ! namtrc_adv in reference namelist 158 NAMELIST/namtrc_adv/ ln_trcadv_NONE, & ! No advection 159 & ln_trcadv_cen, nn_cen_h, nn_cen_v, & ! CEN 160 & ln_trcadv_fct, nn_fct_h, nn_fct_v, & ! FCT 161 & ln_trcadv_mus, ln_mus_ups, & ! MUSCL 162 & ln_trcadv_ubs, nn_ubs_v, & ! UBS 163 & ln_trcadv_qck ! QCK 164 !!---------------------------------------------------------------------- 165 ! 166 ! !== Namelist ==! 167 REWIND( numnat_ref ) ! namtrc_adv in reference namelist 162 168 READ ( numnat_ref, namtrc_adv, IOSTAT = ios, ERR = 901) 163 169 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_adv in reference namelist', lwp ) 164 165 REWIND( numnat_cfg ) ! namtrc_adv in configuration namelist170 ! 171 REWIND( numnat_cfg ) ! namtrc_adv in configuration namelist 166 172 READ ( numnat_cfg, namtrc_adv, IOSTAT = ios, ERR = 902 ) 167 173 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_adv in configuration namelist', lwp ) 168 174 IF(lwm) WRITE ( numont, namtrc_adv ) 169 170 IF(lwp) THEN ! Namelist print175 ! 176 IF(lwp) THEN ! Namelist print 171 177 WRITE(numout,*) 172 178 WRITE(numout,*) 'trc_adv_ini : choice/control of the tracer advection scheme' 173 179 WRITE(numout,*) '~~~~~~~~~~~' 174 180 WRITE(numout,*) ' Namelist namtrc_adv : chose a advection scheme for tracers' 181 WRITE(numout,*) ' No advection on passive tracers ln_trcadv_NONE= ', ln_trcadv_NONE 175 182 WRITE(numout,*) ' centered scheme ln_trcadv_cen = ', ln_trcadv_cen 176 183 WRITE(numout,*) ' horizontal 2nd/4th order nn_cen_h = ', nn_fct_h … … 179 186 WRITE(numout,*) ' horizontal 2nd/4th order nn_fct_h = ', nn_fct_h 180 187 WRITE(numout,*) ' vertical 2nd/4th order nn_fct_v = ', nn_fct_v 181 WRITE(numout,*) ' 2nd order + vertical sub-timestepping nn_fct_zts = ', nn_fct_zts182 188 WRITE(numout,*) ' MUSCL scheme ln_trcadv_mus = ', ln_trcadv_mus 183 189 WRITE(numout,*) ' + upstream scheme near river mouths ln_mus_ups = ', ln_mus_ups … … 187 193 ENDIF 188 194 ! 189 190 ioptio = 0 !== Parameter control ==! 191 IF( ln_trcadv_cen ) ioptio = ioptio + 1 192 IF( ln_trcadv_fct ) ioptio = ioptio + 1 193 IF( ln_trcadv_mus ) ioptio = ioptio + 1 194 IF( ln_trcadv_ubs ) ioptio = ioptio + 1 195 IF( ln_trcadv_qck ) ioptio = ioptio + 1 196 197 ! 198 IF( ioptio == 0 ) THEN 199 nadv = np_NO_adv 200 CALL ctl_warn( 'trc_adv_init: You are running without tracer advection.' ) 201 ENDIF 202 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE advection scheme in namelist namtrc_adv' ) 195 ! !== Parameter control & set nadv ==! 196 ioptio = 0 197 IF( ln_trcadv_NONE ) THEN ; ioptio = ioptio + 1 ; nadv = np_NO_adv ; ENDIF 198 IF( ln_trcadv_cen ) THEN ; ioptio = ioptio + 1 ; nadv = np_CEN ; ENDIF 199 IF( ln_trcadv_fct ) THEN ; ioptio = ioptio + 1 ; nadv = np_FCT ; ENDIF 200 IF( ln_trcadv_mus ) THEN ; ioptio = ioptio + 1 ; nadv = np_MUS ; ENDIF 201 IF( ln_trcadv_ubs ) THEN ; ioptio = ioptio + 1 ; nadv = np_UBS ; ENDIF 202 IF( ln_trcadv_qck ) THEN ; ioptio = ioptio + 1 ; nadv = np_QCK ; ENDIF 203 ! 204 IF( ioptio /= 1 ) CALL ctl_stop( 'trc_adv_ini: Choose ONE advection option in namelist namtrc_adv' ) 203 205 ! 204 206 IF( ln_trcadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 ) & 205 207 .AND. ( nn_cen_v /= 2 .AND. nn_cen_v /= 4 ) ) THEN 206 CALL ctl_stop( 'trc_adv_ini t: CEN scheme, choose 2nd or 4th order' )208 CALL ctl_stop( 'trc_adv_ini: CEN scheme, choose 2nd or 4th order' ) 207 209 ENDIF 208 210 IF( ln_trcadv_fct .AND. ( nn_fct_h /= 2 .AND. nn_fct_h /= 4 ) & 209 211 .AND. ( nn_fct_v /= 2 .AND. nn_fct_v /= 4 ) ) THEN 210 CALL ctl_stop( 'trc_adv_init: FCT scheme, choose 2nd or 4th order' ) 211 ENDIF 212 IF( ln_trcadv_fct .AND. nn_fct_zts > 0 ) THEN 213 IF( nn_fct_h == 4 ) THEN 214 nn_fct_h = 2 215 CALL ctl_stop( 'trc_adv_init: force 2nd order FCT scheme, 4th order does not exist with sub-timestepping' ) 216 ENDIF 217 IF( .NOT.ln_linssh ) THEN 218 CALL ctl_stop( 'trc_adv_init: vertical sub-timestepping not allow in non-linear free surface' ) 219 ENDIF 220 IF( nn_fct_zts == 1 ) CALL ctl_warn( 'trc_adv_init: FCT with ONE sub-timestep = FCT without sub-timestep' ) 212 CALL ctl_stop( 'trc_adv_ini: FCT scheme, choose 2nd or 4th order' ) 221 213 ENDIF 222 214 IF( ln_trcadv_ubs .AND. ( nn_ubs_v /= 2 .AND. nn_ubs_v /= 4 ) ) THEN 223 CALL ctl_stop( 'trc_adv_ini t: UBS scheme, choose 2nd or 4th order' )215 CALL ctl_stop( 'trc_adv_ini: UBS scheme, choose 2nd or 4th order' ) 224 216 ENDIF 225 217 IF( ln_trcadv_ubs .AND. nn_ubs_v == 4 ) THEN 226 CALL ctl_warn( 'trc_adv_ini t: UBS scheme, only 2nd FCT scheme available on the vertical. It will be used' )218 CALL ctl_warn( 'trc_adv_ini: UBS scheme, only 2nd FCT scheme available on the vertical. It will be used' ) 227 219 ENDIF 228 220 IF( ln_isfcav ) THEN ! ice-shelf cavities 229 IF( ln_trcadv_cen .AND. nn_cen_v /= 4 .OR. & ! NO 4th order with ISF 230 & ln_trcadv_fct .AND. nn_fct_v /= 4 ) CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' ) 231 ENDIF 232 ! 233 ! !== used advection scheme ==! 234 ! ! set nadv 235 IF( ln_trcadv_cen ) nadv = np_CEN 236 IF( ln_trcadv_fct ) nadv = np_FCT 237 IF( ln_trcadv_fct .AND. nn_fct_zts > 0 ) nadv = np_FCT_zts 238 IF( ln_trcadv_mus ) nadv = np_MUS 239 IF( ln_trcadv_ubs ) nadv = np_UBS 240 IF( ln_trcadv_qck ) nadv = np_QCK 241 ! 242 IF(lwp) THEN ! Print the choice 221 IF( ln_trcadv_cen .AND. nn_cen_v == 4 .OR. & ! NO 4th order with ISF 222 & ln_trcadv_fct .AND. nn_fct_v == 4 ) CALL ctl_stop( 'tra_adv_ini: 4th order COMPACT scheme not allowed with ISF' ) 223 ENDIF 224 ! 225 ! !== Print the choice ==! 226 IF(lwp) THEN 243 227 WRITE(numout,*) 244 IF( nadv == np_NO_adv ) WRITE(numout,*) ' NO passive tracer advection' 245 IF( nadv == np_CEN ) WRITE(numout,*) ' CEN scheme is used. Horizontal order: ', nn_cen_h, & 246 & ' Vertical order: ', nn_cen_v 247 IF( nadv == np_FCT ) WRITE(numout,*) ' FCT scheme is used. Horizontal order: ', nn_fct_h, & 248 & ' Vertical order: ', nn_fct_v 249 IF( nadv == np_FCT_zts ) WRITE(numout,*) ' use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping' 250 IF( nadv == np_MUS ) WRITE(numout,*) ' MUSCL scheme is used' 251 IF( nadv == np_UBS ) WRITE(numout,*) ' UBS scheme is used' 252 IF( nadv == np_QCK ) WRITE(numout,*) ' QUICKEST scheme is used' 228 SELECT CASE ( nadv ) 229 CASE( np_NO_adv ) ; WRITE(numout,*) ' ===>> NO passive tracer advection' 230 CASE( np_CEN ) ; WRITE(numout,*) ' ===>> CEN scheme is used. Horizontal order: ', nn_cen_h, & 231 & ' Vertical order: ', nn_cen_v 232 CASE( np_FCT ) ; WRITE(numout,*) ' ===>> FCT scheme is used. Horizontal order: ', nn_fct_h, & 233 & ' Vertical order: ', nn_fct_v 234 CASE( np_MUS ) ; WRITE(numout,*) ' ===>> MUSCL scheme is used' 235 CASE( np_UBS ) ; WRITE(numout,*) ' ===>> UBS scheme is used' 236 CASE( np_QCK ) ; WRITE(numout,*) ' ===>> QUICKEST scheme is used' 237 END SELECT 253 238 ENDIF 254 239 ! 255 240 END SUBROUTINE trc_adv_ini 256 241 257 #else258 !!----------------------------------------------------------------------259 !! Default option Empty module260 !!----------------------------------------------------------------------261 CONTAINS262 SUBROUTINE trc_adv( kt )263 INTEGER, INTENT(in) :: kt264 WRITE(*,*) 'trc_adv: You should not have seen this print! error?', kt265 END SUBROUTINE trc_adv266 242 #endif 267 243
Note: See TracChangeset
for help on using the changeset viewer.
