- Timestamp:
- 2015-10-01T14:48:08+02:00 (9 years ago)
- Location:
- branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO
- Files:
-
- 1 added
- 3 deleted
- 25 edited
- 2 moved
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/dom_oce.F90
r5753 r5770 21 21 22 22 IMPLICIT NONE 23 PUBLIC ! allows the acces to par_oce when dom_oce is used 24 ! ! exception to coding rules... to be suppressed ??? 23 PUBLIC ! allows the acces to par_oce when dom_oce is used (exception to coding rules) 25 24 26 25 PUBLIC dom_oce_alloc ! Called from nemogcm.F90 … … 107 106 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: rdttra !: vertical profile of tracer time step 108 107 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: r2dtra !: = 2*rdttra except at nit000 (=rdttra) if neuler=0 109 110 ! !!* Namelist namcla : cross land advection111 INTEGER, PUBLIC :: nn_cla !: =1 cross land advection for exchanges through some straits (ORCA2)112 108 113 109 !!---------------------------------------------------------------------- -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
r5363 r5770 81 81 ENDIF 82 82 ! 83 CALL dom_nam ! read namelist ( namrun, namdom , namcla)83 CALL dom_nam ! read namelist ( namrun, namdom ) 84 84 CALL dom_clo ! Closed seas and lake 85 85 CALL dom_hgr ! Horizontal mesh … … 131 131 !! ** input : - namrun namelist 132 132 !! - namdom namelist 133 !! - namcla namelist134 133 !! - namnc4 namelist ! "key_netcdf4" only 135 134 !!---------------------------------------------------------------------- … … 146 145 & ppsur, ppa0, ppa1, ppkth, ppacr, ppdzmin, pphmax, ldbletanh, & 147 146 & ppa2, ppkth2, ppacr2 148 NAMELIST/namcla/ nn_cla149 147 #if defined key_netcdf4 150 148 NAMELIST/namnc4/ nn_nchunks_i, nn_nchunks_j, nn_nchunks_k, ln_nc4zip … … 305 303 rdth = rn_rdth 306 304 307 REWIND( numnam_ref ) ! Namelist namcla in reference namelist : Cross land advection308 READ ( numnam_ref, namcla, IOSTAT = ios, ERR = 905)309 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in reference namelist', lwp )310 311 REWIND( numnam_cfg ) ! Namelist namcla in configuration namelist : Cross land advection312 READ ( numnam_cfg, namcla, IOSTAT = ios, ERR = 906 )313 906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namcla in configuration namelist', lwp )314 IF(lwm) WRITE( numond, namcla )315 316 IF(lwp) THEN317 WRITE(numout,*)318 WRITE(numout,*) ' Namelist namcla'319 WRITE(numout,*) ' cross land advection nn_cla = ', nn_cla320 ENDIF321 IF ( nn_cla .EQ. 1 ) THEN322 IF ( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2323 CONTINUE324 ELSE325 CALL ctl_stop( 'STOP', 'Cross land advation iplemented only for ORCA2 configuration: cp_cfg = "orca" and jp_cfg = 2 ' )326 ENDIF327 ENDIF328 329 305 #if defined key_netcdf4 330 306 ! ! NetCDF 4 case ("key_netcdf4" defined) -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/dommsk.F90
r5552 r5770 17 17 !! - ! 2005-11 (V. Garnier) Surface pressure gradient organization 18 18 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 19 !! 3.6 ! 2015-05 (P. Mathiot) ISF: add wmask,wumask and wvmask 19 20 !!---------------------------------------------------------------------- 20 21 … … 129 130 !! tmask_i : interior ocean mask 130 131 !!---------------------------------------------------------------------- 131 ! 132 INTEGER :: ji, jj, jk ! dummy loop indices 132 INTEGER :: ji, jj, jk ! dummy loop indices 133 133 INTEGER :: iif, iil, ii0, ii1, ii ! local integers 134 134 INTEGER :: ijf, ijl, ij0, ij1 ! - - … … 284 284 ! 3. Ocean/land mask at wu-, wv- and w points 285 285 !---------------------------------------------- 286 wmask (:,:,1) = tmask(:,:,1) ! ????????287 wumask(:,:,1) = umask(:,:,1) ! ????????288 wvmask(:,:,1) = vmask(:,:,1) ! ????????289 DO jk =2,jpk290 wmask (:,:,jk) =tmask(:,:,jk) * tmask(:,:,jk-1)291 wumask(:,:,jk) =umask(:,:,jk) * umask(:,:,jk-1)292 wvmask(:,:,jk) =vmask(:,:,jk) * vmask(:,:,jk-1)286 wmask (:,:,1) = tmask(:,:,1) ! surface 287 wumask(:,:,1) = umask(:,:,1) 288 wvmask(:,:,1) = vmask(:,:,1) 289 DO jk = 2, jpk ! interior values 290 wmask (:,:,jk) = tmask(:,:,jk) * tmask(:,:,jk-1) 291 wumask(:,:,jk) = umask(:,:,jk) * umask(:,:,jk-1) 292 wvmask(:,:,jk) = vmask(:,:,jk) * vmask(:,:,jk-1) 293 293 END DO 294 294 … … 377 377 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA_R2 configuration 378 378 ! ! Increased lateral friction near of some straits 379 IF( nn_cla == 0 ) THEN 380 ! ! Gibraltar strait : partial slip (fmask=0.5) 381 ij0 = 101 ; ij1 = 101 382 ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp 383 ij0 = 102 ; ij1 = 102 384 ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp 385 ! 386 ! ! Bab el Mandeb : partial slip (fmask=1) 387 ij0 = 87 ; ij1 = 88 388 ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp 389 ij0 = 88 ; ij1 = 88 390 ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp 391 ! 392 ENDIF 379 ! ! Gibraltar strait : partial slip (fmask=0.5) 380 ij0 = 101 ; ij1 = 101 381 ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp 382 ij0 = 102 ; ij1 = 102 383 ii0 = 139 ; ii1 = 140 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 0.5_wp 384 ! 385 ! ! Bab el Mandeb : partial slip (fmask=1) 386 ij0 = 87 ; ij1 = 88 387 ii0 = 160 ; ii1 = 160 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp 388 ij0 = 88 ; ij1 = 88 389 ii0 = 159 ; ii1 = 159 ; fmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 1:jpk ) = 1._wp 390 ! 393 391 ! ! Danish straits : strong slip (fmask > 2) 394 392 ! We keep this as an example but it is instable in this case -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domngb.F90
r3294 r5770 37 37 !! -> not good if located at too high latitude... 38 38 !!---------------------------------------------------------------------- 39 !40 39 REAL(wp) , INTENT(in ) :: plon, plat ! longitude,latitude of the point 41 40 INTEGER , INTENT( out) :: kii, kjj ! i-,j-index of the closes grid point … … 49 48 IF( nn_timing == 1 ) CALL timing_start('dom_ngb') 50 49 ! 51 CALL wrk_alloc( jpi, jpj,zglam, zgphi, zmask, zdist )50 CALL wrk_alloc( jpi,jpj, zglam, zgphi, zmask, zdist ) 52 51 ! 53 52 zmask(:,:) = 0._wp … … 76 75 ENDIF 77 76 ! 78 CALL wrk_dealloc( jpi, jpj,zglam, zgphi, zmask, zdist )77 CALL wrk_dealloc( jpi,jpj, zglam, zgphi, zmask, zdist ) 79 78 ! 80 79 IF( nn_timing == 1 ) CALL timing_stop('dom_ngb') -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domwri.F90
r5737 r5770 36 36 !!---------------------------------------------------------------------- 37 37 CONTAINS 38 39 40 38 41 39 SUBROUTINE dom_wri_coordinate -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5656 r5770 501 501 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 502 502 ! ! ===================== 503 IF( nn_cla == 0 ) THEN 504 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 505 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 506 DO ji = mi0(ii0), mi1(ii1) 507 DO jj = mj0(ij0), mj1(ij1) 508 mbathy(ji,jj) = 15 509 END DO 503 ! 504 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 505 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 506 DO ji = mi0(ii0), mi1(ii1) 507 DO jj = mj0(ij0), mj1(ij1) 508 mbathy(ji,jj) = 15 510 509 END DO 511 IF(lwp) WRITE(numout,*)512 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0513 !514 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open515 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)516 DO ji = mi0(ii0), mi1(ii1)517 DO jj = mj0(ij0), mj1(ij1)518 mbathy(ji,jj) = 12519 END DO510 END DO 511 IF(lwp) WRITE(numout,*) 512 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 513 ! 514 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 515 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 516 DO ji = mi0(ii0), mi1(ii1) 517 DO jj = mj0(ij0), mj1(ij1) 518 mbathy(ji,jj) = 12 520 519 END DO 521 IF(lwp) WRITE(numout,*)522 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0523 ENDIF520 END DO 521 IF(lwp) WRITE(numout,*) 522 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 524 523 ! 525 524 ENDIF … … 545 544 ! 546 545 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 547 ! 548 IF( nn_cla == 0 ) THEN 549 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 550 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 551 DO ji = mi0(ii0), mi1(ii1) 552 DO jj = mj0(ij0), mj1(ij1) 553 bathy(ji,jj) = 284._wp 554 END DO 546 ! 547 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 548 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 549 DO ji = mi0(ii0), mi1(ii1) 550 DO jj = mj0(ij0), mj1(ij1) 551 bathy(ji,jj) = 284._wp 555 552 END DO 556 IF(lwp) WRITE(numout,*)557 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0558 !559 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open560 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)561 DO ji = mi0(ii0), mi1(ii1)562 DO jj = mj0(ij0), mj1(ij1)563 bathy(ji,jj) = 137._wp564 END DO553 END DO 554 IF(lwp) WRITE(numout,*) 555 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 556 ! 557 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 558 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 559 DO ji = mi0(ii0), mi1(ii1) 560 DO jj = mj0(ij0), mj1(ij1) 561 bathy(ji,jj) = 137._wp 565 562 END DO 566 IF(lwp) WRITE(numout,*)567 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0568 ENDIF563 END DO 564 IF(lwp) WRITE(numout,*) 565 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 569 566 ! 570 567 ENDIF -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/dtatsd.F90
r5215 r5770 174 174 END DO 175 175 END DO 176 IF( nn_cla == 1 ) THEN ! Cross Land advection 177 il0 = 138 ; il1 = 138 ! set T & S profile at Gibraltar Strait 178 ij0 = 101 ; ij1 = 102 179 ii0 = 139 ; ii1 = 139 180 DO jl = mi0(il0), mi1(il1) 181 DO jj = mj0(ij0), mj1(ij1) 182 DO ji = mi0(ii0), mi1(ii1) 183 sf_tsd(jp_tem)%fnow(ji,jj,:) = sf_tsd(jp_tem)%fnow(jl,jj,:) 184 sf_tsd(jp_sal)%fnow(ji,jj,:) = sf_tsd(jp_sal)%fnow(jl,jj,:) 185 END DO 186 END DO 187 END DO 188 il0 = 164 ; il1 = 164 ! set T & S profile at Bab el Mandeb Strait 189 ij0 = 87 ; ij1 = 88 190 ii0 = 161 ; ii1 = 163 191 DO jl = mi0(il0), mi1(il1) 192 DO jj = mj0(ij0), mj1(ij1) 193 DO ji = mi0(ii0), mi1(ii1) 194 sf_tsd(jp_tem)%fnow(ji,jj,:) = sf_tsd(jp_tem)%fnow(jl,jj,:) 195 sf_tsd(jp_sal)%fnow(ji,jj,:) = sf_tsd(jp_sal)%fnow(jl,jj,:) 196 END DO 197 END DO 198 END DO 199 ELSE ! No Cross Land advection 200 ij0 = 87 ; ij1 = 96 ! Reduced temperature in Red Sea 201 ii0 = 148 ; ii1 = 160 202 sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 4:10 ) = 7.0_wp 203 sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp 204 sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp 205 ENDIF 176 ij0 = 87 ; ij1 = 96 ! Reduced temperature in Red Sea 177 ii0 = 148 ; ii1 = 160 178 sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 4:10 ) = 7.0_wp 179 sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 11:13 ) = 6.5_wp 180 sf_tsd(jp_tem)%fnow( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , 14:20 ) = 6.0_wp 206 181 ENDIF 207 182 ! -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/istate.F90
r5766 r5770 34 34 USE dtatsd ! data temperature and salinity (dta_tsd routine) 35 35 USE dtauvd ! data: U & V current (dta_uvd routine) 36 USE zpshde ! partial step: hor. derivative (zps_hde routine)37 USE eosbn2 ! equation of state (eos bn2 routine)38 36 USE domvvl ! varying vertical mesh 39 37 USE dynspg_oce ! pressure gradient schemes … … 119 117 ENDIF 120 118 IF ( ln_uvd_init .AND. lk_c1d ) THEN ! read 3D U and V data at nit000 121 CALL wrk_alloc( jpi, jpj, jpk, 2,zuvd )119 CALL wrk_alloc( jpi,jpj,jpk,2, zuvd ) 122 120 CALL dta_uvd( nit000, zuvd ) 123 121 ub(:,:,:) = zuvd(:,:,:,1) ; un(:,:,:) = ub(:,:,:) 124 122 vb(:,:,:) = zuvd(:,:,:,2) ; vn(:,:,:) = vb(:,:,:) 125 CALL wrk_dealloc( jpi, jpj, jpk, 2,zuvd )123 CALL wrk_dealloc( jpi,jpj,jpk,2, zuvd ) 126 124 ENDIF 127 125 ENDIF 128 ! 126 ! 129 127 ! - ML - sshn could be modified by istate_eel, so that initialization of fse3t_b is done here 130 128 IF( lk_vvl ) THEN 131 129 DO jk = 1, jpk 132 130 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 133 END DO131 END DO 134 132 ENDIF 135 133 ! … … 456 454 !!---------------------------------------------------------------------- 457 455 ! 458 CALL wrk_alloc( jpi, jpj, jpk,zprn)456 CALL wrk_alloc( jpi,jpj,jpk, zprn) 459 457 ! 460 458 IF(lwp) WRITE(numout,*) … … 553 551 rotb (:,:,:) = rotn (:,:,:) ! set the before to the now value 554 552 ! 555 CALL wrk_dealloc( jpi, jpj, jpk,zprn)553 CALL wrk_dealloc( jpi,jpj,jpk, zprn) 556 554 ! 557 555 END SUBROUTINE istate_uvg -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
r5737 r5770 18 18 !! - ! 2010-10 (R. Furner, G. Madec) runoff and cla added directly here 19 19 !! 3.6 ! 2014-11 (P. Mathiot) isf added directly here 20 !! 3.7 ! 2015-10 (G. Madec) remove cross-land advection 20 21 !!---------------------------------------------------------------------- 21 22 … … 29 30 USE sbcrnf ! river runoff 30 31 USE sbcisf ! ice shelf 31 USE cla ! cross land advection (cla_div routine)32 32 USE in_out_manager ! I/O manager 33 33 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 68 68 !! - compute the now divergence given by : 69 69 !! hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 70 !! correct hdiv with runoff inflow (div_rnf), ice shelf melting (div_isf) 71 !! and cross land flow (div_cla) 70 !! correct hdiv with runoff inflow (div_rnf) and ice shelf melting (div_isf) 72 71 !! II. vorticity : 73 72 !! - save the curl computed at the previous time-step … … 226 225 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 227 226 IF( ln_divisf .AND. (nn_isf /= 0) ) CALL sbc_isf_div( hdivn ) ! ice shelf (update hdivn field) 228 IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (Update Hor. divergence)229 227 230 228 ! 4. Lateral boundary conditions on hdivn and rotn … … 256 254 !! - compute the now divergence given by : 257 255 !! hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 258 !! correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla)256 !! correct hdiv with runoff inflow (div_rnf) 259 257 !! - Relavtive Vorticity : 260 258 !! - save the curl computed at the previous time-step (rotb = rotn) … … 325 323 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 326 324 IF( ln_divisf .AND. (nn_isf .GT. 0) ) CALL sbc_isf_div( hdivn ) ! ice shelf (update hdivn field) 327 IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field)328 325 ! 329 326 CALL lbc_lnk( hdivn, 'T', 1. ) ; CALL lbc_lnk( rotn , 'F', 1. ) ! lateral boundary cond. (no sign change) -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r5120 r5770 271 271 CALL solver_init( nit000 ) ! Elliptic solver initialisation 272 272 #endif 273 274 ! ! Control of timestep choice275 IF( lk_dynspg_ts .OR. lk_dynspg_exp ) THEN276 IF( nn_cla == 1 ) CALL ctl_stop( 'Crossland advection not implemented for this free surface formulation' )277 ENDIF278 279 273 ! ! Control of hydrostatic pressure choice 280 IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) THEN 281 CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 282 ENDIF 274 IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 283 275 ! 284 276 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_init') -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90
r4990 r5770 14 14 !! 3.2 ! 2009-03 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 15 15 !! 3.7 ! 2014-04 (F. Roquet, G. Madec) add some trends diag 16 !! - ! 2014-12 (G. Madec) remove cross-land advection (cla) 16 17 !!---------------------------------------------------------------------- 17 18 #if defined key_dynspg_flt || defined key_esopa … … 36 37 USE bdydyn ! ocean open boundary condition on dynamics 37 38 USE bdyvol ! ocean open boundary condition (bdy_vol routine) 38 USE cla ! cross land advection39 39 USE trd_oce ! trends: ocean variables 40 40 USE trddyn ! trend manager: dynamics … … 206 206 CALL Agrif_dyn( kt ) ! Update velocities on each coarse/fine interfaces 207 207 #endif 208 IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 ) CALL cla_dynspg( kt ) ! Cross Land Advection (update (ua,va))209 208 210 209 ! compute the next vertically averaged velocity (effect of the additional force not included) -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r5737 r5770 80 80 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 81 81 ! 82 CALL wrk_alloc( jpi, jpj,zhdiv )82 CALL wrk_alloc( jpi,jpj, zhdiv ) 83 83 ! 84 84 IF( kt == nit000 ) THEN 85 !86 85 IF(lwp) WRITE(numout,*) 87 86 IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 88 87 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 89 !90 88 ENDIF 91 89 ! … … 159 157 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 160 158 !!---------------------------------------------------------------------- 161 ! 162 INTEGER, INTENT(in) :: kt ! time step 159 INTEGER, INTENT(in) :: kt ! time step 160 ! 161 INTEGER :: ji, jj, jk ! dummy loop indices 162 REAL(wp) :: z1_2dt ! local scalars 163 163 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 164 164 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv 165 ! 166 INTEGER :: ji, jj, jk ! dummy loop indices 167 REAL(wp) :: z1_2dt ! local scalars 168 !!---------------------------------------------------------------------- 169 165 !!---------------------------------------------------------------------- 166 ! 170 167 IF( nn_timing == 1 ) CALL timing_start('wzv') 171 168 ! 172 169 IF( kt == nit000 ) THEN 173 !174 170 IF(lwp) WRITE(numout,*) 175 171 IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' … … 177 173 ! 178 174 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 179 !180 175 ENDIF 181 176 ! !------------------------------! … … 216 211 217 212 #if defined key_bdy 218 IF (lk_bdy) THEN213 IF( lk_bdy ) THEN 219 214 DO jk = 1, jpkm1 220 215 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) … … 224 219 ! 225 220 IF( nn_timing == 1 ) CALL timing_stop('wzv') 226 227 221 ! 228 222 END SUBROUTINE wzv 223 229 224 230 225 SUBROUTINE ssh_swp( kt ) … … 265 260 sshb(:,:) = sshn(:,:) ! before <-- now 266 261 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 262 ! 267 263 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 268 264 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv.F90
r5758 r5770 7 7 !! 3.3 ! 2010-09 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport 8 8 !! 3.6 ! 2011-06 (G. Madec) Addition of Mixed Layer Eddy parameterisation 9 !! 4.0! 2014-05 (G. Madec) Add 2nd/4th order cases for CEN and FCT schemes9 !! 3.7 ! 2014-05 (G. Madec) Add 2nd/4th order cases for CEN and FCT schemes 10 10 !! - ! 2014-12 (G. Madec) suppression of cross land advection option 11 11 !!---------------------------------------------------------------------- 12 12 13 13 !!---------------------------------------------------------------------- 14 !! tra_adv : compute ocean tracer advection trend 15 !! tra_adv_ctl : control the different options of advection scheme 16 !!---------------------------------------------------------------------- 17 USE oce ! ocean dynamics and active tracers 18 USE dom_oce ! ocean space and time domain 19 USE domvvl ! variable vertical scale factors 20 USE traadv_cen2 ! 2nd order centered scheme (tra_adv_cen2 routine) 21 USE traadv_tvd ! TVD scheme (tra_adv_tvd routine) 22 USE traadv_muscl ! MUSCL scheme (tra_adv_muscl routine) 23 USE traadv_muscl2 ! MUSCL2 scheme (tra_adv_muscl2 routine) 24 USE traadv_ubs ! UBS scheme (tra_adv_ubs routine) 25 USE traadv_qck ! QUICKEST scheme (tra_adv_qck routine) 26 USE traadv_mle ! ML eddy induced velocity (tra_adv_mle routine) 27 USE cla ! cross land advection (cla_traadv routine) 28 USE ldftra ! lateral diffusion coefficient on tracers 29 USE ldfslp ! Lateral diffusion: slopes of neutral surfaces 14 !! tra_adv : compute ocean tracer advection trend 15 !! tra_adv_ctl : control the different options of advection scheme 16 !!---------------------------------------------------------------------- 17 USE oce ! ocean dynamics and active tracers 18 USE dom_oce ! ocean space and time domain 19 USE domvvl ! variable vertical scale factors 20 USE traadv_cen ! centered scheme (tra_adv_cen routine) 21 USE traadv_fct ! FCT scheme (tra_adv_fct routine) 22 USE traadv_mus ! MUSCL scheme (tra_adv_mus routine) 23 USE traadv_ubs ! UBS scheme (tra_adv_ubs routine) 24 USE traadv_qck ! QUICKEST scheme (tra_adv_qck routine) 25 USE traadv_mle ! ML eddy induced velocity (tra_adv_mle routine) 26 USE ldftra ! lateral diffusion: eddy diffusivity & EIV coeff. 27 USE ldfslp ! Lateral diffusion: slopes of neutral surfaces 30 28 ! 31 USE in_out_manager 32 USE iom 33 USE prtctl 34 USE lib_mpp 35 USE wrk_nemo 36 USE timing 37 USE sbc_oce 29 USE in_out_manager ! I/O manager 30 USE iom ! I/O module 31 USE prtctl ! Print control 32 USE lib_mpp ! MPP library 33 USE wrk_nemo ! Memory Allocation 34 USE timing ! Timing 35 38 36 USE diaptr ! Poleward heat transport 39 40 37 41 38 IMPLICIT NONE … … 45 42 PUBLIC tra_adv_init ! routine called by opa module 46 43 47 ! !!* Namelist namtra_adv * 48 LOGICAL :: ln_traadv_cen2 ! 2nd order centered scheme flag 49 LOGICAL :: ln_traadv_tvd ! TVD scheme flag 50 LOGICAL :: ln_traadv_tvd_zts ! TVD scheme flag with vertical sub time-stepping 51 LOGICAL :: ln_traadv_muscl ! MUSCL scheme flag 52 LOGICAL :: ln_traadv_muscl2 ! MUSCL2 scheme flag 53 LOGICAL :: ln_traadv_ubs ! UBS scheme flag 54 LOGICAL :: ln_traadv_qck ! QUICKEST scheme flag 55 LOGICAL :: ln_traadv_msc_ups ! use upstream scheme within muscl 56 57 58 INTEGER :: nadv ! choice of the type of advection scheme 59 44 ! !!* Namelist namtra_adv * 45 LOGICAL :: ln_traadv_cen ! centered scheme flag 46 INTEGER :: nn_cen_h, nn_cen_v ! =2/4 : horizontal and vertical choices of the order of CEN scheme 47 LOGICAL :: ln_traadv_fct ! FCT scheme flag 48 INTEGER :: nn_fct_h, nn_fct_v ! =2/4 : horizontal and vertical choices of the order of FCT scheme 49 INTEGER :: nn_fct_zts ! >=1 : 2nd order FCT with vertical sub-timestepping 50 LOGICAL :: ln_traadv_mus ! MUSCL scheme flag 51 LOGICAL :: ln_mus_ups ! use upstream scheme in vivcinity of river mouths 52 LOGICAL :: ln_traadv_ubs ! UBS scheme flag 53 INTEGER :: nn_ubs_v ! =2/4 : vertical choice of the order of UBS scheme 54 LOGICAL :: ln_traadv_qck ! QUICKEST scheme flag 55 56 INTEGER :: nadv ! choice of the type of advection scheme 57 ! 58 ! ! associated indices: 59 INTEGER, PARAMETER :: np_NO_adv = 0 ! no T-S advection 60 INTEGER, PARAMETER :: np_CEN = 1 ! 2nd/4th order centered scheme 61 INTEGER, PARAMETER :: np_FCT = 2 ! 2nd/4th order Flux Corrected Transport scheme 62 INTEGER, PARAMETER :: np_FCT_zts = 3 ! 2nd order FCT scheme with vertical sub-timestepping 63 INTEGER, PARAMETER :: np_MUS = 4 ! MUSCL scheme 64 INTEGER, PARAMETER :: np_UBS = 5 ! 3rd order Upstream Biased Scheme 65 INTEGER, PARAMETER :: np_QCK = 6 ! QUICK scheme 66 60 67 !! * Substitutions 61 68 # include "domzgr_substitute.h90" 62 69 # include "vectopt_loop_substitute.h90" 63 70 !!---------------------------------------------------------------------- 64 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)71 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 65 72 !! $Id$ 66 73 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 84 91 IF( nn_timing == 1 ) CALL timing_start('tra_adv') 85 92 ! 86 CALL wrk_alloc( jpi, jpj, jpk,zun, zvn, zwn )93 CALL wrk_alloc( jpi,jpj,jpk, zun, zvn, zwn ) 87 94 ! 88 95 ! ! set time step … … 93 100 ENDIF 94 101 ! 95 IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 ) CALL cla_traadv( kt ) !== Cross Land Advection ==! (hor. advection) 96 ! 97 ! !== effective transport ==! 102 ! !== effective transport ==! 98 103 DO jk = 1, jpkm1 99 104 zun(:,:,jk) = e2u (:,:) * fse3u(:,:,jk) * un(:,:,jk) ! eulerian transport only … … 102 107 END DO 103 108 ! 104 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 109 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! add z-tilde and/or vvl corrections 105 110 zun(:,:,:) = zun(:,:,:) + un_td(:,:,:) 106 111 zvn(:,:,:) = zvn(:,:,:) + vn_td(:,:,:) 107 112 ENDIF 108 113 ! 109 zun(:,:,jpk) = 0._wp ! no transport trough the bottom110 zvn(:,:,jpk) = 0._wp ! no transport trough the bottom111 zwn(:,:,jpk) = 0._wp ! no transport trough the bottom114 zun(:,:,jpk) = 0._wp ! no transport trough the bottom 115 zvn(:,:,jpk) = 0._wp 116 zwn(:,:,jpk) = 0._wp 112 117 ! 113 118 IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad ) & … … 120 125 CALL iom_put( "wocetr_eff", zwn ) 121 126 ! 127 !!gm ??? 122 128 IF( ln_diaptr ) CALL dia_ptr( zvn ) ! diagnose the effective MSF 123 ! 124 125 SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==! 126 CASE ( 1 ) ; CALL tra_adv_cen2 ( kt, nit000, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! 2nd order centered 127 CASE ( 2 ) ; CALL tra_adv_tvd ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! TVD 128 CASE ( 3 ) ; CALL tra_adv_muscl ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsa, jpts, ln_traadv_msc_ups ) ! MUSCL 129 CASE ( 4 ) ; CALL tra_adv_muscl2 ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! MUSCL2 130 CASE ( 5 ) ; CALL tra_adv_ubs ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! UBS 131 CASE ( 6 ) ; CALL tra_adv_qck ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! QUICKEST 132 CASE ( 7 ) ; CALL tra_adv_tvd_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) ! TVD ZTS 133 ! 134 CASE (-1 ) !== esopa: test all possibility with control print ==! 135 CALL tra_adv_cen2 ( kt, nit000, 'TRA', zun, zvn, zwn, tsb, tsn, tsa, jpts ) 136 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv0 - Ta: ', mask1=tmask, & 137 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 138 CALL tra_adv_tvd ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) 139 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv1 - Ta: ', mask1=tmask, & 140 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 141 CALL tra_adv_muscl ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsa, jpts, ln_traadv_msc_ups ) 142 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv3 - Ta: ', mask1=tmask, & 143 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 144 CALL tra_adv_muscl2( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) 145 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv4 - Ta: ', mask1=tmask, & 146 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 147 CALL tra_adv_ubs ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) 148 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv5 - Ta: ', mask1=tmask, & 149 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 150 CALL tra_adv_qck ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) 151 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' adv6 - Ta: ', mask1=tmask, & 152 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 129 !!gm ??? 130 ! 131 SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==! 132 ! 133 CASE ( np_CEN ) ! Centered scheme : 2nd / 4th order 134 CALL tra_adv_cen ( kt, nit000, 'TRA', zun, zvn, zwn , tsn, tsa, jpts, nn_cen_h, nn_cen_v ) 135 CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order 136 CALL tra_adv_fct ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts, nn_fct_h, nn_fct_v ) 137 CASE ( np_FCT_zts ) ! 2nd order FCT with vertical time-splitting 138 CALL tra_adv_fct_zts( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts , nn_fct_zts ) 139 CASE ( np_MUS ) ! MUSCL 140 CALL tra_adv_mus ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsa, jpts , ln_mus_ups ) 141 CASE ( np_UBS ) ! UBS 142 CALL tra_adv_ubs ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts , nn_ubs_v ) 143 CASE ( np_QCK ) ! QUICKEST 144 CALL tra_adv_qck ( kt, nit000, 'TRA', r2dtra, zun, zvn, zwn, tsb, tsn, tsa, jpts ) 145 ! 153 146 END SELECT 154 147 ! … … 159 152 IF( nn_timing == 1 ) CALL timing_stop( 'tra_adv' ) 160 153 ! 161 CALL wrk_dealloc( jpi, jpj, jpk,zun, zvn, zwn )154 CALL wrk_dealloc( jpi,jpj,jpk, zun, zvn, zwn ) 162 155 ! 163 156 END SUBROUTINE tra_adv … … 171 164 !! tracer advection schemes and set nadv 172 165 !!---------------------------------------------------------------------- 173 INTEGER :: ioptio 174 INTEGER :: ios ! Local integer output status for namelist read 175 !! 176 NAMELIST/namtra_adv/ ln_traadv_cen2 , ln_traadv_tvd, & 177 & ln_traadv_muscl, ln_traadv_muscl2, & 178 & ln_traadv_ubs , ln_traadv_qck, & 179 & ln_traadv_msc_ups, ln_traadv_tvd_zts 180 !!---------------------------------------------------------------------- 181 182 REWIND( numnam_ref ) ! Namelist namtra_adv in reference namelist : Tracer advection scheme 166 INTEGER :: ioptio, ios ! Local integers 167 ! 168 NAMELIST/namtra_adv/ ln_traadv_cen, nn_cen_h, nn_cen_v, & ! CEN 169 & ln_traadv_fct, nn_fct_h, nn_fct_v, nn_fct_zts, & ! FCT 170 & ln_traadv_mus, ln_mus_ups, & ! MUSCL 171 & ln_traadv_ubs, nn_ubs_v, & ! UBS 172 & ln_traadv_qck ! QCK 173 !!---------------------------------------------------------------------- 174 ! 175 ! !== Namelist ==! 176 ! 177 REWIND( numnam_ref ) ! Namelist namtra_adv in reference namelist : Tracer advection scheme 183 178 READ ( numnam_ref, namtra_adv, IOSTAT = ios, ERR = 901) 184 179 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in reference namelist', lwp ) 185 186 REWIND( numnam_cfg ) ! Namelist namtra_adv in configuration namelist : Tracer advection scheme180 ! 181 REWIND( numnam_cfg ) ! Namelist namtra_adv in configuration namelist : Tracer advection scheme 187 182 READ ( numnam_cfg, namtra_adv, IOSTAT = ios, ERR = 902 ) 188 183 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtra_adv in configuration namelist', lwp ) 189 184 IF(lwm) WRITE ( numond, namtra_adv ) 190 185 191 IF(lwp) THEN ! Namelist print186 IF(lwp) THEN ! Namelist print 192 187 WRITE(numout,*) 193 188 WRITE(numout,*) 'tra_adv_init : choice/control of the tracer advection scheme' 194 189 WRITE(numout,*) '~~~~~~~~~~~' 195 190 WRITE(numout,*) ' Namelist namtra_adv : chose a advection scheme for tracers' 196 WRITE(numout,*) ' 2nd order advection scheme ln_traadv_cen2 = ', ln_traadv_cen2 197 WRITE(numout,*) ' TVD advection scheme ln_traadv_tvd = ', ln_traadv_tvd 198 WRITE(numout,*) ' MUSCL advection scheme ln_traadv_muscl = ', ln_traadv_muscl 199 WRITE(numout,*) ' MUSCL2 advection scheme ln_traadv_muscl2 = ', ln_traadv_muscl2 200 WRITE(numout,*) ' UBS advection scheme ln_traadv_ubs = ', ln_traadv_ubs 201 WRITE(numout,*) ' QUICKEST advection scheme ln_traadv_qck = ', ln_traadv_qck 202 WRITE(numout,*) ' upstream scheme within muscl ln_traadv_msc_ups = ', ln_traadv_msc_ups 203 WRITE(numout,*) ' TVD advection scheme with zts ln_traadv_tvd_zts = ', ln_traadv_tvd_zts 204 ENDIF 205 206 ioptio = 0 ! Parameter control 207 IF( ln_traadv_cen2 ) ioptio = ioptio + 1 208 IF( ln_traadv_tvd ) ioptio = ioptio + 1 209 IF( ln_traadv_muscl ) ioptio = ioptio + 1 210 IF( ln_traadv_muscl2 ) ioptio = ioptio + 1 211 IF( ln_traadv_ubs ) ioptio = ioptio + 1 212 IF( ln_traadv_qck ) ioptio = ioptio + 1 213 IF( ln_traadv_tvd_zts) ioptio = ioptio + 1 214 IF( lk_esopa ) ioptio = 1 215 216 IF( ( ln_traadv_muscl .OR. ln_traadv_muscl2 .OR. ln_traadv_ubs .OR. ln_traadv_qck .OR. ln_traadv_tvd_zts ) & 217 .AND. ln_isfcav ) CALL ctl_stop( 'Only traadv_cen2 and traadv_tvd is compatible with ice shelf cavity') 218 219 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE advection scheme in namelist namtra_adv' ) 220 221 ! ! Set nadv 222 IF( ln_traadv_cen2 ) nadv = 1 223 IF( ln_traadv_tvd ) nadv = 2 224 IF( ln_traadv_muscl ) nadv = 3 225 IF( ln_traadv_muscl2 ) nadv = 4 226 IF( ln_traadv_ubs ) nadv = 5 227 IF( ln_traadv_qck ) nadv = 6 228 IF( ln_traadv_tvd_zts) nadv = 7 229 IF( lk_esopa ) nadv = -1 230 231 IF(lwp) THEN ! Print the choice 191 WRITE(numout,*) ' centered scheme ln_traadv_cen = ', ln_traadv_cen 192 WRITE(numout,*) ' horizontal 2nd/4th order nn_cen_h = ', nn_fct_h 193 WRITE(numout,*) ' vertical 2nd/4th order nn_cen_v = ', nn_fct_v 194 WRITE(numout,*) ' Flux Corrected Transport scheme ln_traadv_fct = ', ln_traadv_fct 195 WRITE(numout,*) ' horizontal 2nd/4th order nn_fct_h = ', nn_fct_h 196 WRITE(numout,*) ' vertical 2nd/4th order nn_fct_v = ', nn_fct_v 197 WRITE(numout,*) ' 2nd order + vertical sub-timestepping nn_fct_zts = ', nn_fct_zts 198 WRITE(numout,*) ' MUSCL scheme ln_traadv_mus = ', ln_traadv_mus 199 WRITE(numout,*) ' + upstream scheme near river mouths ln_mus_ups = ', ln_mus_ups 200 WRITE(numout,*) ' UBS scheme ln_traadv_ubs = ', ln_traadv_ubs 201 WRITE(numout,*) ' vertical 2nd/4th order nn_ubs_v = ', nn_ubs_v 202 WRITE(numout,*) ' QUICKEST scheme ln_traadv_qck = ', ln_traadv_qck 203 ENDIF 204 205 ioptio = 0 !== Parameter control ==! 206 IF( ln_traadv_cen ) ioptio = ioptio + 1 207 IF( ln_traadv_fct ) ioptio = ioptio + 1 208 IF( ln_traadv_mus ) ioptio = ioptio + 1 209 IF( ln_traadv_ubs ) ioptio = ioptio + 1 210 IF( ln_traadv_qck ) ioptio = ioptio + 1 211 ! 212 IF( ioptio == 0 ) THEN 213 nadv = np_NO_adv 214 CALL ctl_warn( 'tra_adv_init: You are running without tracer advection.' ) 215 ENDIF 216 IF( ioptio /= 1 ) CALL ctl_stop( 'tra_adv_init: Choose ONE advection scheme in namelist namtra_adv' ) 217 ! 218 IF( ln_traadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 ) & ! Centered 219 .AND. ( nn_cen_v /= 2 .AND. nn_cen_v /= 4 ) ) THEN 220 CALL ctl_stop( 'tra_adv_init: CEN scheme, choose 2nd or 4th order' ) 221 ENDIF 222 IF( ln_traadv_fct .AND. ( nn_fct_h /= 2 .AND. nn_fct_h /= 4 ) & ! FCT 223 .AND. ( nn_fct_v /= 2 .AND. nn_fct_v /= 4 ) ) THEN 224 CALL ctl_stop( 'tra_adv_init: FCT scheme, choose 2nd or 4th order' ) 225 ENDIF 226 IF( ln_traadv_fct .AND. nn_fct_zts > 0 ) THEN 227 IF( nn_fct_h == 4 ) THEN 228 nn_fct_h = 2 229 CALL ctl_stop( 'tra_adv_init: force 2nd order FCT scheme, 4th order does not exist with sub-timestepping' ) 230 ENDIF 231 IF( lk_vvl ) THEN 232 CALL ctl_stop( 'tra_adv_init: vertical sub-timestepping not allow in non-linear free surface' ) 233 ENDIF 234 IF( nn_fct_zts == 1 ) CALL ctl_warn( 'tra_adv_init: FCT with ONE sub-timestep = FCT without sub-timestep' ) 235 ENDIF 236 IF( ln_traadv_ubs .AND. ( nn_ubs_v /= 2 .AND. nn_ubs_v /= 4 ) ) THEN ! UBS 237 CALL ctl_stop( 'tra_adv_init: UBS scheme, choose 2nd or 4th order' ) 238 ENDIF 239 IF( ln_traadv_ubs .AND. nn_ubs_v == 4 ) THEN 240 CALL ctl_warn( 'tra_adv_init: UBS scheme, only 2nd FCT scheme available on the vertical. It will be used' ) 241 ENDIF 242 IF( ln_isfcav ) THEN ! ice-shelf cavities 243 IF( ln_traadv_cen .AND. nn_cen_v /= 4 .OR. & ! NO 4th order with ISF 244 & ln_traadv_fct .AND. nn_fct_v /= 4 ) CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' ) 245 ENDIF 246 ! 247 ! !== used advection scheme ==! 248 ! ! set nadv 249 IF( ln_traadv_cen ) nadv = np_CEN 250 IF( ln_traadv_fct ) nadv = np_FCT 251 IF( ln_traadv_fct .AND. nn_fct_zts > 0 ) nadv = np_FCT_zts 252 IF( ln_traadv_mus ) nadv = np_MUS 253 IF( ln_traadv_ubs ) nadv = np_UBS 254 IF( ln_traadv_qck ) nadv = np_QCK 255 256 IF(lwp) THEN ! Print the choice 232 257 WRITE(numout,*) 233 IF( nadv == 1 ) WRITE(numout,*) ' 2nd order scheme is used' 234 IF( nadv == 2 ) WRITE(numout,*) ' TVD scheme is used' 235 IF( nadv == 3 ) WRITE(numout,*) ' MUSCL scheme is used' 236 IF( nadv == 4 ) WRITE(numout,*) ' MUSCL2 scheme is used' 237 IF( nadv == 5 ) WRITE(numout,*) ' UBS scheme is used' 238 IF( nadv == 6 ) WRITE(numout,*) ' QUICKEST scheme is used' 239 IF( nadv == 7 ) WRITE(numout,*) ' TVD ZTS scheme is used' 240 IF( nadv == -1 ) WRITE(numout,*) ' esopa test: use all advection scheme' 241 ENDIF 242 ! 243 CALL tra_adv_mle_init ! initialisation of the Mixed Layer Eddy parametrisation (MLE) 258 IF( nadv == np_NO_adv ) WRITE(numout,*) ' NO T-S advection' 259 IF( nadv == np_CEN ) WRITE(numout,*) ' CEN scheme is used. Horizontal order: ', nn_cen_h, & 260 & ' Vertical order: ', nn_cen_v 261 IF( nadv == np_FCT ) WRITE(numout,*) ' FCT scheme is used. Horizontal order: ', nn_fct_h, & 262 & ' Vertical order: ', nn_fct_v 263 IF( nadv == np_FCT_zts ) WRITE(numout,*) ' use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping' 264 IF( nadv == np_MUS ) WRITE(numout,*) ' MUSCL scheme is used' 265 IF( nadv == np_UBS ) WRITE(numout,*) ' UBS scheme is used' 266 IF( nadv == np_QCK ) WRITE(numout,*) ' QUICKEST scheme is used' 267 ENDIF 268 ! 269 CALL tra_adv_mle_init !== initialisation of the Mixed Layer Eddy parametrisation (MLE) ==! 244 270 ! 245 271 END SUBROUTINE tra_adv_init -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90
r5737 r5770 1 MODULE traadv_ tvd1 MODULE traadv_fct 2 2 !!============================================================================== 3 !! *** MODULE traadv_ tvd***4 !! Ocean tracers: horizontal & vertical advective trend 3 !! *** MODULE traadv_fct *** 4 !! Ocean tracers: horizontal & vertical advective trend (2nd/4th order Flux Corrected Transport method) 5 5 !!============================================================================== 6 !! History : OPA ! 1995-12 (L. Mortier) Original code 7 !! ! 2000-01 (H. Loukos) adapted to ORCA 8 !! ! 2000-10 (MA Foujols E.Kestenare) include file not routine 9 !! ! 2000-12 (E. Kestenare M. Levy) fix bug in trtrd indexes 10 !! ! 2001-07 (E. Durand G. Madec) adaptation to ORCA config 11 !! 8.5 ! 2002-06 (G. Madec) F90: Free form and module 12 !! NEMO 1.0 ! 2004-01 (A. de Miranda, G. Madec, J.M. Molines ): advective bbl 13 !! 2.0 ! 2008-04 (S. Cravatte) add the i-, j- & k- trends computation 14 !! - ! 2009-11 (V. Garnier) Surface pressure gradient organization 15 !! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport 6 !! History : 3.7 ! 2015-09 (L. Debreu, G. Madec) original code (inspired from traadv_tvd.F90) 16 7 !!---------------------------------------------------------------------- 17 8 18 9 !!---------------------------------------------------------------------- 19 !! tra_adv_tvd : update the tracer trend with the 3D advection trends using a TVD scheme 20 !! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm 10 !! tra_adv_fct : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme 11 !! tra_adv_fct_zts: update the tracer trend with a 3D advective trends using a 2nd order FCT scheme 12 !! with sub-time-stepping in the vertical direction 13 !! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm 14 !! interp_4th_cpt : 4th order compact scheme for the vertical component of the advection 21 15 !!---------------------------------------------------------------------- 22 16 USE oce ! ocean dynamics and active tracers … … 28 22 USE diaptr ! poleward transport diagnostics 29 23 ! 24 USE in_out_manager ! I/O manager 30 25 USE lib_mpp ! MPP library 31 26 USE lbclnk ! ocean lateral boundary condition (or mpp link) 32 USE in_out_manager ! I/O manager27 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 33 28 USE wrk_nemo ! Memory Allocation 34 29 USE timing ! Timing 35 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)36 30 37 31 IMPLICIT NONE 38 32 PRIVATE 39 33 40 PUBLIC tra_adv_tvd ! routine called by traadv.F90 41 PUBLIC tra_adv_tvd_zts ! routine called by traadv.F90 42 43 LOGICAL :: l_trd ! flag to compute trends 34 PUBLIC tra_adv_fct ! routine called by traadv.F90 35 PUBLIC tra_adv_fct_zts ! routine called by traadv.F90 36 PUBLIC interp_4th_cpt ! routine called by traadv_cen.F90 37 38 LOGICAL :: l_trd ! flag to compute trends 39 REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6 44 40 45 41 !! * Substitutions … … 47 43 # include "vectopt_loop_substitute.h90" 48 44 !!---------------------------------------------------------------------- 49 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)45 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 50 46 !! $Id$ 51 47 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 53 49 CONTAINS 54 50 55 SUBROUTINE tra_adv_ tvd ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,&56 & ptb, ptn, pta, kjpt)57 !!---------------------------------------------------------------------- 58 !! *** ROUTINE tra_adv_ tvd***51 SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 52 & ptb, ptn, pta, kjpt, kn_fct_h, kn_fct_v ) 53 !!---------------------------------------------------------------------- 54 !! *** ROUTINE tra_adv_fct *** 59 55 !! 60 56 !! ** Purpose : Compute the now trend due to total advection of 61 57 !! tracers and add it to the general trend of tracer equations 62 58 !! 63 !! ** Method : TVD scheme, i.e. 2nd order centered scheme with 64 !! corrected flux (monotonic correction) 65 !! note: - this advection scheme needs a leap-frog time scheme 59 !! ** Method : - 2nd or 4th FCT scheme on the horizontal direction 60 !! (choice through the value of kn_fct) 61 !! - 4th order compact scheme on the vertical 62 !! - corrected flux (monotonic correction) 66 63 !! 67 64 !! ** Action : - update (pta) with the now advective tracer trends 68 !! - save the trends 69 !!---------------------------------------------------------------------- 70 USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace 71 ! 65 !! - send the trends for further diagnostics 66 !!---------------------------------------------------------------------- 72 67 INTEGER , INTENT(in ) :: kt ! ocean time-step index 73 68 INTEGER , INTENT(in ) :: kit000 ! first time step index 74 69 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 75 70 INTEGER , INTENT(in ) :: kjpt ! number of tracers 71 INTEGER , INTENT(in ) :: kn_fct_h ! order of the FCT scheme (=2 or 4) 72 INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4) 76 73 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 77 74 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components … … 79 76 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 80 77 ! 81 INTEGER :: ji, jj, jk, jn ! dummy loop indices 82 INTEGER :: ik 83 REAL(wp) :: z2dtt, zbtr, ztra ! local scalar 84 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - - 85 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - - 86 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz 87 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz 88 !!---------------------------------------------------------------------- 89 ! 90 IF( nn_timing == 1 ) CALL timing_start('tra_adv_tvd') 91 ! 92 CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz ) 78 INTEGER :: ji, jj, jk, jn ! dummy loop indices 79 REAL(wp) :: z2dtt, ztra ! local scalar 80 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u ! - - 81 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v ! - - 82 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 83 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz 84 !!---------------------------------------------------------------------- 85 ! 86 IF( nn_timing == 1 ) CALL timing_start('tra_adv_fct') 87 ! 88 CALL wrk_alloc( jpi,jpj,jpk, zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) 93 89 ! 94 90 IF( kt == kit000 ) THEN 95 91 IF(lwp) WRITE(numout,*) 96 IF(lwp) WRITE(numout,*) 'tra_adv_ tvd : TVDadvection scheme on ', cdtype92 IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype 97 93 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 98 !99 l_trd = .FALSE.100 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.101 94 ENDIF 95 ! 96 l_trd = .FALSE. 97 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 102 98 ! 103 99 IF( l_trd ) THEN 104 100 CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 105 ztrdx(:,:,:) = 0. e0 ; ztrdy(:,:,:) = 0.e0 ; ztrdz(:,:,:) = 0.e0101 ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp 106 102 ENDIF 107 103 ! 108 zwi(:,:,:) = 0.e0 ; 109 ! 104 ! ! surface & bottom value : flux set to zero one for all 105 IF( lk_vvl ) zwz(:,:, 1 ) = 0._wp ! except at the surface in linear free surface case 106 zwx(:,:,jpk) = 0._wp ; zwy(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp 107 ! 108 zwi(:,:,:) = 0._wp 110 109 ! ! =========== 111 110 DO jn = 1, kjpt ! tracer loop 112 111 ! ! =========== 113 ! 1. Bottom and k=1 value : flux set to zero 114 ! ---------------------------------- 115 zwx(:,:,jpk) = 0.e0 ; zwz(:,:,jpk) = 0.e0 116 zwy(:,:,jpk) = 0.e0 ; zwi(:,:,jpk) = 0.e0 117 118 zwz(:,:,1 ) = 0._wp 119 ! 2. upstream advection with initial mass fluxes & intermediate update 120 ! -------------------------------------------------------------------- 121 ! upstream tracer flux in the i and j direction 112 ! 113 ! !== upstream advection with initial mass fluxes & intermediate update ==! 114 ! !* upstream tracer flux in the i and j direction 122 115 DO jk = 1, jpkm1 123 116 DO jj = 1, jpjm1 … … 133 126 END DO 134 127 END DO 135 136 ! upstream tracer flux in the k direction 137 ! Interior value 138 DO jk = 2, jpkm1 128 ! !* upstream tracer flux in the k direction *! 129 DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask) 139 130 DO jj = 1, jpj 140 131 DO ji = 1, jpi … … 145 136 END DO 146 137 END DO 147 ! Surface value 148 IF( lk_vvl ) THEN 149 IF ( ln_isfcav ) THEN 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 zwz(ji,jj, mikt(ji,jj) ) = 0.e0 ! volume variable 153 END DO 154 END DO 155 ELSE 156 zwz(:,:,1) = 0.e0 ! volume variable 157 END IF 158 ELSE 159 IF ( ln_isfcav ) THEN 138 ! 139 IF(.NOT.lk_vvl ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked) 140 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 160 141 DO jj = 1, jpj 161 142 DO ji = 1, jpi … … 163 144 END DO 164 145 END DO 165 ELSE 166 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface167 END 146 ELSE ! no cavities: only at the ocean surface 147 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 148 ENDIF 168 149 ENDIF 169 170 ! total advective trend 171 DO jk = 1, jpkm1 150 ! 151 DO jk = 1, jpkm1 !* trend and after field with monotonic scheme 172 152 z2dtt = p2dt(jk) 173 153 DO jj = 2, jpjm1 174 154 DO ji = fs_2, fs_jpim1 ! vector opt. 175 zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )176 155 ! total intermediate advective trends 177 ztra = - zbtr *( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &178 & 179 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )156 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 157 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 158 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 180 159 ! update and guess with monotonic sheme 160 !!gm why tmask added in the two following lines ??? the mask is done in tranxt ! 181 161 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk) 182 162 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) … … 184 164 END DO 185 165 END DO 186 ! ! Lateral boundary conditions on zwi (unchanged sign) 187 CALL lbc_lnk( zwi, 'T', 1. ) 188 189 ! ! trend diagnostics (contribution of upstream fluxes) 190 IF( l_trd ) THEN 191 ! store intermediate advective trends 166 CALL lbc_lnk( zwi, 'T', 1. ) ! Lateral boundary conditions on zwi (unchanged sign) 167 ! 168 IF( l_trd ) THEN ! trend diagnostics (contribution of upstream fluxes) 192 169 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 193 170 END IF 194 ! 171 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 195 172 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 196 173 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 197 174 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 198 175 ENDIF 199 200 ! 3. antidiffusive flux : high order minus low order 201 ! -------------------------------------------------- 202 ! antidiffusive flux on i and j 203 DO jk = 1, jpkm1 204 DO jj = 1, jpjm1 205 DO ji = 1, fs_jpim1 ! vector opt. 206 zwx(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk) 207 zwy(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk) 208 END DO 209 END DO 210 END DO 211 212 ! antidiffusive flux on k 213 ! Interior value 214 DO jk = 2, jpkm1 215 DO jj = 1, jpj 216 DO ji = 1, jpi 217 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk) 218 END DO 219 END DO 220 END DO 221 ! surface value 222 IF ( ln_isfcav ) THEN 223 DO jj = 1, jpj 224 DO ji = 1, jpi 225 zwz(ji,jj,mikt(ji,jj)) = 0.e0 226 END DO 227 END DO 228 ELSE 229 zwz(:,:,1) = 0.e0 230 END IF 176 ! 177 ! 178 ! !== anti-diffusive flux : high order minus low order ==! 179 ! 180 SELECT CASE( kn_fct_h ) !* horizontal anti-diffusive fluxes 181 ! 182 CASE( 2 ) ! 2nd order centered 183 DO jk = 1, jpkm1 184 DO jj = 1, jpjm1 185 DO ji = 1, fs_jpim1 ! vector opt. 186 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk) 187 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk) 188 END DO 189 END DO 190 END DO 191 ! 192 CASE( 4 ) ! 4th order centered 193 zltu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 194 zltv(:,:,jpk) = 0._wp 195 DO jk = 1, jpkm1 ! Laplacian 196 DO jj = 1, jpjm1 ! First derivative (gradient) 197 DO ji = 1, fs_jpim1 ! vector opt. 198 ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 199 ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 200 END DO 201 END DO 202 DO jj = 2, jpjm1 ! 203 DO ji = fs_2, fs_jpim1 ! vector opt. 204 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6 205 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6 206 END DO 207 END DO 208 END DO 209 ! 210 CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn) 211 ! 212 DO jk = 1, jpkm1 ! Horizontal advective fluxes 213 DO jj = 1, jpjm1 214 DO ji = 1, fs_jpim1 ! vector opt. 215 zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points 216 zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) 217 ! ! C4 minus upstream advective fluxes 218 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk) 219 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk) 220 END DO 221 END DO 222 END DO 223 ! 224 CASE( 41 ) ! 4th order centered ==>> !!gm coding attempt need to be tested 225 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 226 ztv(:,:,jpk) = 0._wp 227 DO jk = 1, jpkm1 ! gradient 228 DO jj = 1, jpjm1 ! First derivative (gradient) 229 DO ji = 1, fs_jpim1 ! vector opt. 230 ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 231 ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 232 END DO 233 END DO 234 END DO 235 CALL lbc_lnk( ztu, 'U', -1. ) ; CALL lbc_lnk( ztv, 'V', -1. ) ! Lateral boundary cond. (unchanged sgn) 236 ! 237 DO jk = 1, jpkm1 ! Horizontal advective fluxes 238 DO jj = 2, jpjm1 239 DO ji = 2, fs_jpim1 ! vector opt. 240 zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points (x2) 241 zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) 242 ! ! C4 interpolation of T at u- & v-points (x2) 243 zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) ) 244 zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) ) 245 ! ! C4 minus upstream advective fluxes 246 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 247 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 248 END DO 249 END DO 250 END DO 251 ! 252 END SELECT 253 ! !* vertical anti-diffusive fluxes 254 SELECT CASE( kn_fct_v ) ! Interior values (w-masked) 255 ! 256 CASE( 2 ) ! 2nd order centered 257 DO jk = 2, jpkm1 258 DO jj = 2, jpjm1 259 DO ji = fs_2, fs_jpim1 260 zwz(ji,jj,jk) = ( 0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) & 261 - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 262 END DO 263 END DO 264 END DO 265 ! 266 CASE( 4 ) ! 4th order COMPACT 267 ! 268 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! COMPACT interpolation of T at w-point 269 ! 270 DO jk = 2, jpkm1 271 DO jj = 2, jpjm1 272 DO ji = fs_2, fs_jpim1 273 zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 274 END DO 275 END DO 276 END DO 277 ! 278 END SELECT 279 ! ! top ocean value: high order = upstream ==>> zwz=0 280 zwz(:,:, 1 ) = 0._wp ! only ocean surface as interior zwz values have been w-masked 281 ! 231 282 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions 232 283 CALL lbc_lnk( zwz, 'W', 1. ) 233 284 234 ! 4. monotonicity algorithm235 ! -------------------------285 ! !== monotonicity algorithm ==! 286 ! 236 287 CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 237 288 238 289 239 ! 5. final trend with corrected fluxes240 ! ------------------------------------290 ! !== final trend with corrected fluxes ==! 291 ! 241 292 DO jk = 1, jpkm1 242 293 DO jj = 2, jpjm1 243 294 DO ji = fs_2, fs_jpim1 ! vector opt. 244 zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 245 ! total advective trends 246 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 247 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 248 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 249 ! add them to the general tracer trends 250 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk) 251 END DO 252 END DO 253 END DO 254 255 ! ! trend diagnostics (contribution of upstream fluxes) 256 IF( l_trd ) THEN 295 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 296 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 297 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) & 298 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 299 END DO 300 END DO 301 END DO 302 ! 303 IF( l_trd ) THEN ! trend diagnostics (contribution of upstream fluxes) 257 304 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed 258 305 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed 259 306 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 260 261 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 262 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 263 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 307 ! 308 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 309 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 310 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 311 ! 312 CALL wrk_dealloc( jpi,jpj,jpk, ztrdx, ztrdy, ztrdz ) 264 313 END IF 265 314 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) … … 271 320 END DO 272 321 ! 273 CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz ) 274 IF( l_trd ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 275 ! 276 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd') 277 ! 278 END SUBROUTINE tra_adv_tvd 279 280 281 SUBROUTINE tra_adv_tvd_zts ( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 282 & ptb, ptn, pta, kjpt ) 283 !!---------------------------------------------------------------------- 284 !! *** ROUTINE tra_adv_tvd_zts *** 322 CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw ) 323 ! 324 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_fct') 325 ! 326 END SUBROUTINE tra_adv_fct 327 328 329 SUBROUTINE tra_adv_fct_zts( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 330 & ptb, ptn, pta, kjpt, kn_fct_zts ) 331 !!---------------------------------------------------------------------- 332 !! *** ROUTINE tra_adv_fct_zts *** 285 333 !! 286 334 !! ** Purpose : Compute the now trend due to total advection of … … 296 344 !! - save the trends 297 345 !!---------------------------------------------------------------------- 298 USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace299 !300 346 INTEGER , INTENT(in ) :: kt ! ocean time-step index 301 347 INTEGER , INTENT(in ) :: kit000 ! first time step index 302 348 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 303 349 INTEGER , INTENT(in ) :: kjpt ! number of tracers 350 INTEGER , INTENT(in ) :: kn_fct_zts ! number of number of vertical sub-timesteps 304 351 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 305 352 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components … … 310 357 REAL(wp), DIMENSION( jpk ) :: zr_p2dt ! reciprocal of tracer timestep 311 358 INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices 312 INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection313 359 INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps 314 360 INTEGER :: jtaken ! toggle for collecting appropriate fluxes from sub timesteps 315 361 REAL(wp) :: z_rzts ! Fractional length of Euler forward sub-timestep for vertical advection 316 REAL(wp) :: z2dtt, z btr, ztra! local scalar362 REAL(wp) :: z2dtt, ztra ! local scalar 317 363 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - - 318 364 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - - 319 REAL(wp), POINTER, DIMENSION(:,: ) ::zwx_sav , zwy_sav320 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz, zhdiv, zwz_sav, zwzts321 REAL(wp), POINTER, DIMENSION(:,:,:) ::ztrdx, ztrdy, ztrdz322 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs323 !!---------------------------------------------------------------------- 324 ! 325 IF( nn_timing == 1 ) CALL timing_start('tra_adv_ tvd_zts')326 ! 327 CALL wrk_alloc( jpi, jpj,zwx_sav, zwy_sav )328 CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz , zhdiv, zwz_sav, zwzts)329 CALL wrk_alloc( jpi, jpj, jpk, 3,ztrs )365 REAL(wp), POINTER, DIMENSION(:,: ) :: zwx_sav , zwy_sav 366 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwx, zwy, zwz, zhdiv, zwzts, zwz_sav 367 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz 368 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs 369 !!---------------------------------------------------------------------- 370 ! 371 IF( nn_timing == 1 ) CALL timing_start('tra_adv_fct_zts') 372 ! 373 CALL wrk_alloc( jpi,jpj, zwx_sav, zwy_sav ) 374 CALL wrk_alloc( jpi,jpj, jpk, zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 375 CALL wrk_alloc( jpi,jpj,jpk,kjpt+1, ztrs ) 330 376 ! 331 377 IF( kt == kit000 ) THEN 332 378 IF(lwp) WRITE(numout,*) 333 IF(lwp) WRITE(numout,*) 'tra_adv_ tvd_zts : TVD ZTS advection schemeon ', cdtype379 IF(lwp) WRITE(numout,*) 'tra_adv_fct_zts : 2nd order FCT scheme with ', kn_fct_zts, ' vertical sub-timestep on ', cdtype 334 380 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 335 381 ENDIF 336 382 ! 337 383 l_trd = .FALSE. 338 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.384 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 339 385 ! 340 386 IF( l_trd ) THEN 341 CALL wrk_alloc( jpi, jpj, jpk,ztrdx, ztrdy, ztrdz )387 CALL wrk_alloc( jpi,jpj,jpk, ztrdx, ztrdy, ztrdz ) 342 388 ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp 343 389 ENDIF 344 390 ! 345 391 zwi(:,:,:) = 0._wp 346 z_rzts = 1._wp / REAL( jnzts, wp )392 z_rzts = 1._wp / REAL( kn_fct_zts, wp ) 347 393 zr_p2dt(:) = 1._wp / p2dt(:) 348 394 ! … … 373 419 374 420 ! upstream tracer flux in the k direction 375 ! Interior value 376 DO jk = 2, jpkm1 421 DO jk = 2, jpkm1 ! Interior value 377 422 DO jj = 1, jpj 378 423 DO ji = 1, jpi 379 424 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 380 425 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 381 zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 382 END DO 383 END DO 384 END DO 385 ! Surface value 386 IF( lk_vvl ) THEN 387 IF ( ln_isfcav ) THEN 426 zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 427 END DO 428 END DO 429 END DO 430 ! ! top value 431 IF( lk_vvl ) THEN ! variable volume: only k=1 as zwz is multiplied by wmask 432 zwz(:,:, 1 ) = 0._wp 433 ELSE ! linear free surface 434 IF( ln_isfcav ) THEN ! ice-shelf cavities 388 435 DO jj = 1, jpj 389 436 DO ji = 1, jpi 390 zwz(ji,jj, mikt(ji,jj) ) = 0.e0 ! volume variable + isf 391 END DO 392 END DO 393 ELSE 394 zwz(:,:,1) = 0.e0 ! volume variable + no isf 395 END IF 396 ELSE 397 IF ( ln_isfcav ) THEN 398 DO jj = 1, jpj 399 DO ji = 1, jpi 400 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface + isf 401 END DO 402 END DO 403 ELSE 404 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface + no isf 405 END IF 437 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) 438 END DO 439 END DO 440 ELSE ! standard case 441 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 442 ENDIF 406 443 ENDIF 407 408 ! total advective trend 409 DO jk = 1, jpkm1 444 ! 445 DO jk = 1, jpkm1 ! total advective trend 410 446 z2dtt = p2dt(jk) 411 447 DO jj = 2, jpjm1 412 448 DO ji = fs_2, fs_jpim1 ! vector opt. 413 zbtr = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )414 449 ! total intermediate advective trends 415 ztra = - zbtr *( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &416 & 417 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )450 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 451 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 452 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 418 453 ! update and guess with monotonic sheme 419 454 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra … … 422 457 END DO 423 458 END DO 424 ! ! Lateral boundary conditions on zwi (unchanged sign) 425 CALL lbc_lnk( zwi, 'T', 1. ) 426 427 ! ! trend diagnostics (contribution of upstream fluxes) 428 IF( l_trd ) THEN 429 ! store intermediate advective trends 459 ! 460 CALL lbc_lnk( zwi, 'T', 1. ) ! Lateral boundary conditions on zwi (unchanged sign) 461 ! 462 IF( l_trd ) THEN ! trend diagnostics (contribution of upstream fluxes) 430 463 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 431 464 END IF 432 ! 465 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 433 466 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 434 467 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) … … 436 469 ENDIF 437 470 438 ! 3. antidiffusive flux : high order minus low order 439 ! -------------------------------------------------- 440 ! antidiffusive flux on i and j 441 442 443 DO jk = 1, jpkm1 444 471 ! 3. anti-diffusive flux : high order minus low order 472 ! --------------------------------------------------- 473 474 DO jk = 1, jpkm1 !* horizontal anti-diffusive fluxes 475 ! 445 476 DO jj = 1, jpjm1 446 477 DO ji = 1, fs_jpim1 ! vector opt. 447 478 zwx_sav(ji,jj) = zwx(ji,jj,jk) 448 479 zwy_sav(ji,jj) = zwy(ji,jj,jk) 449 480 ! 450 481 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) 451 482 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) 452 483 END DO 453 484 END DO 454 455 DO jj = 2, jpjm1 ! partial horizontal divergence485 ! 486 DO jj = 2, jpjm1 ! partial horizontal divergence 456 487 DO ji = fs_2, fs_jpim1 457 488 zhdiv(ji,jj,jk) = ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & … … 459 490 END DO 460 491 END DO 461 492 ! 462 493 DO jj = 1, jpjm1 463 494 DO ji = 1, fs_jpim1 ! vector opt. 464 zwx(ji,jj,jk) = zwx(ji,jj,jk) 465 zwy(ji,jj,jk) = zwy(ji,jj,jk) 495 zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj) 496 zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj) 466 497 END DO 467 498 END DO 468 499 END DO 469 500 470 ! antidiffusive flux on k 471 zwz(:,:,1) = 0._wp ! Surface value 472 zwz_sav(:,:,:) = zwz(:,:,:) 473 ! 474 ztrs(:,:,:,1) = ptb(:,:,:,jn) 475 zwzts(:,:,:) = 0._wp 476 477 DO jl = 1, jnzts ! Start of sub timestepping loop 478 479 IF( jl == 1 ) THEN ! Euler forward to kick things off 480 jtb = 1 ; jtn = 1 ; jta = 2 481 zts(:) = p2dt(:) * z_rzts 482 jtaken = MOD( jnzts + 1 , 2) ! Toggle to collect every second flux 483 ! starting at jl =1 if jnzts is odd; 484 ! starting at jl =2 otherwise 485 ELSEIF( jl == 2 ) THEN ! First leapfrog step 486 jtb = 1 ; jtn = 2 ; jta = 3 487 zts(:) = 2._wp * p2dt(:) * z_rzts 488 ELSE ! Shuffle pointers for subsequent leapfrog steps 489 jtb = MOD(jtb,3) + 1 490 jtn = MOD(jtn,3) + 1 491 jta = MOD(jta,3) + 1 501 ! !* vertical anti-diffusive flux 502 zwz_sav(:,:,:) = zwz(:,:,:) 503 ztrs (:,:,:,1) = ptb(:,:,:,jn) 504 zwzts (:,:,:) = 0._wp 505 IF( lk_vvl ) zwz(:,:, 1 ) = 0._wp ! surface value set to zero in vvl case 506 ! 507 DO jl = 1, kn_fct_zts ! Start of sub timestepping loop 508 ! 509 IF( jl == 1 ) THEN ! Euler forward to kick things off 510 jtb = 1 ; jtn = 1 ; jta = 2 511 zts(:) = p2dt(:) * z_rzts 512 jtaken = MOD( kn_fct_zts + 1 , 2) ! Toggle to collect every second flux 513 ! ! starting at jl =1 if kn_fct_zts is odd; 514 ! ! starting at jl =2 otherwise 515 ELSEIF( jl == 2 ) THEN ! First leapfrog step 516 jtb = 1 ; jtn = 2 ; jta = 3 517 zts(:) = 2._wp * p2dt(:) * z_rzts 518 ELSE ! Shuffle pointers for subsequent leapfrog steps 519 jtb = MOD(jtb,3) + 1 520 jtn = MOD(jtn,3) + 1 521 jta = MOD(jta,3) + 1 492 522 ENDIF 493 DO jk = 2, jpkm1 ! Interior value523 DO jk = 2, jpkm1 ! interior value 494 524 DO jj = 2, jpjm1 495 525 DO ji = fs_2, fs_jpim1 496 zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) 497 IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk)*zts(jk) ! Accumulate time-weighted vertcal flux 498 END DO 499 END DO 500 END DO 501 526 zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) * wmask(ji,jj,jk) 527 IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk) * zts(jk) ! Accumulate time-weighted vertcal flux 528 END DO 529 END DO 530 END DO 531 IF(.NOT.lk_vvl ) THEN ! top value (only in linear free surface case) 532 IF( ln_isfcav ) THEN ! ice-shelf cavities 533 DO jj = 1, jpj 534 DO ji = 1, jpi 535 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface 536 END DO 537 END DO 538 ELSE ! standard case 539 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 540 ENDIF 541 ENDIF 542 ! 502 543 jtaken = MOD( jtaken + 1 , 2 ) 503 504 DO jk = 2, jpkm1 ! Interior value544 ! 545 DO jk = 2, jpkm1 ! total advective trends 505 546 DO jj = 2, jpjm1 506 547 DO ji = fs_2, fs_jpim1 507 zbtr = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 508 ! total advective trends 509 ztra = - zbtr * ( zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 510 ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) + zts(jk) * ztra 511 END DO 512 END DO 513 END DO 514 548 ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) & 549 & - zts(jk) * ( zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) & 550 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 551 END DO 552 END DO 553 END DO 554 ! 515 555 END DO 516 556 … … 518 558 DO jj = 2, jpjm1 519 559 DO ji = fs_2, fs_jpim1 520 zwz(ji,jj,jk) = zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk)560 zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk) 521 561 END DO 522 562 END DO … … 535 575 DO jj = 2, jpjm1 536 576 DO ji = fs_2, fs_jpim1 ! vector opt. 537 zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 538 ! total advective trends 539 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 540 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 541 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 542 ! add them to the general tracer trends 543 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 577 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 578 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) & 579 & / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 544 580 END DO 545 581 END DO … … 551 587 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed 552 588 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 553 589 ! 554 590 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 555 591 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 556 592 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 593 ! 594 CALL wrk_dealloc( jpi,jpj,jpk, ztrdx, ztrdy, ztrdz ) 557 595 END IF 558 596 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) … … 564 602 END DO 565 603 ! 566 CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz, zhdiv, zwz_sav, zwzts ) 567 CALL wrk_dealloc( jpi, jpj, jpk, 3, ztrs ) 568 CALL wrk_dealloc( jpi, jpj, zwx_sav, zwy_sav ) 569 IF( l_trd ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 570 ! 571 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd_zts') 572 ! 573 END SUBROUTINE tra_adv_tvd_zts 604 CALL wrk_alloc( jpi,jpj, zwx_sav, zwy_sav ) 605 CALL wrk_alloc( jpi,jpj, jpk, zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav ) 606 CALL wrk_alloc( jpi,jpj,jpk,kjpt+1, ztrs ) 607 ! 608 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_fct_zts') 609 ! 610 END SUBROUTINE tra_adv_fct_zts 574 611 575 612 … … 683 720 END SUBROUTINE nonosc 684 721 722 723 SUBROUTINE interp_4th_cpt( pt_in, pt_out ) 724 !!---------------------------------------------------------------------- 725 !! *** ROUTINE interp_4th_cpt *** 726 !! 727 !! ** Purpose : Compute the interpolation of tracer at w-point 728 !! 729 !! ** Method : 4th order compact interpolation 730 !!---------------------------------------------------------------------- 731 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! now tracer fields 732 REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! now tracer field interpolated at w-pts 733 ! 734 INTEGER :: ji, jj, jk ! dummy loop integers 735 REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt 736 !!---------------------------------------------------------------------- 737 738 DO jk = 3, jpkm1 !== build the three diagonal matrix ==! 739 DO jj = 1, jpj 740 DO ji = 1, jpi 741 zwd (ji,jj,jk) = 4._wp 742 zwi (ji,jj,jk) = 1._wp 743 zws (ji,jj,jk) = 1._wp 744 zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 745 ! 746 IF( tmask(ji,jj,jk+1) == 0._wp) THEN ! Switch to second order centered at bottom 747 zwd (ji,jj,jk) = 1._wp 748 zwi (ji,jj,jk) = 0._wp 749 zws (ji,jj,jk) = 0._wp 750 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 751 ENDIF 752 END DO 753 END DO 754 END DO 755 ! 756 jk=2 ! Switch to second order centered at top 757 DO jj=1,jpj 758 DO ji=1,jpi 759 zwd (ji,jj,jk) = 1._wp 760 zwi (ji,jj,jk) = 0._wp 761 zws (ji,jj,jk) = 0._wp 762 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 763 END DO 764 END DO 765 ! 766 ! !== tridiagonal solve ==! 767 DO jj = 1, jpj ! first recurrence 768 DO ji = 1, jpi 769 zwt(ji,jj,2) = zwd(ji,jj,2) 770 END DO 771 END DO 772 DO jk = 3, jpkm1 773 DO jj = 1, jpj 774 DO ji = 1, jpi 775 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 776 END DO 777 END DO 778 END DO 779 ! 780 DO jj = 1, jpj ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 781 DO ji = 1, jpi 782 pt_out(ji,jj,2) = zwrm(ji,jj,2) 783 END DO 784 END DO 785 DO jk = 3, jpkm1 786 DO jj = 1, jpj 787 DO ji = 1, jpi 788 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 789 END DO 790 END DO 791 END DO 792 793 DO jj = 1, jpj ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 794 DO ji = 1, jpi 795 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 796 END DO 797 END DO 798 DO jk = jpk-2, 2, -1 799 DO jj = 1, jpj 800 DO ji = 1, jpi 801 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 802 END DO 803 END DO 804 END DO 805 ! 806 END SUBROUTINE interp_4th_cpt 807 685 808 !!====================================================================== 686 END MODULE traadv_ tvd809 END MODULE traadv_fct -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90
r5737 r5770 1 MODULE traadv_mus cl1 MODULE traadv_mus 2 2 !!====================================================================== 3 !! *** MODULE traadv_mus cl***3 !! *** MODULE traadv_mus *** 4 4 !! Ocean tracers: horizontal & vertical advective trend 5 5 !!====================================================================== … … 9 9 !! 3.2 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport 10 10 !! 3.4 ! 2012-06 (P. Oddo, M. Vichi) include the upstream where needed 11 !!---------------------------------------------------------------------- 12 13 !!---------------------------------------------------------------------- 14 !! tra_adv_muscl : update the tracer trend with the horizontal 11 !! 3.7 ! 2015-09 (G. Madec) add the ice-shelf cavities boundary condition 12 !!---------------------------------------------------------------------- 13 14 !!---------------------------------------------------------------------- 15 !! tra_adv_mus : update the tracer trend with the horizontal 15 16 !! and vertical advection trends using MUSCL scheme 16 17 !!---------------------------------------------------------------------- … … 34 35 PRIVATE 35 36 36 PUBLIC tra_adv_mus cl! routine called by traadv.F9037 PUBLIC tra_adv_mus ! routine called by traadv.F90 37 38 38 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits … … 50 51 CONTAINS 51 52 52 SUBROUTINE tra_adv_mus cl( kt, kit000, cdtype, p2dt, pun, pvn, pwn,&53 & 53 SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 54 & ptb, pta, kjpt, ld_msc_ups ) 54 55 !!---------------------------------------------------------------------- 55 !! *** ROUTINE tra_adv_mus cl***56 !! 57 !! ** Purpose : Compute the now trend due to total advection of T and58 !! Susing a MUSCL scheme (Monotone Upstream-centered Scheme for59 !! Conservation Laws) and add it to the general tracer trend.56 !! *** ROUTINE tra_adv_mus *** 57 !! 58 !! ** Purpose : Compute the now trend due to total advection of tracers 59 !! using a MUSCL scheme (Monotone Upstream-centered Scheme for 60 !! Conservation Laws) and add it to the general tracer trend. 60 61 !! 61 62 !! ** Method : MUSCL scheme plus centered scheme at ocean boundaries 63 !! ld_msc_ups=T : 62 64 !! 63 65 !! ** Action : - update (ta,sa) with the now advective tracer trends 64 !! - s ave trends66 !! - send trends to trdtra module for further diagnostcs 65 67 !! 66 68 !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation … … 77 79 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 78 80 ! 79 INTEGER :: ji, jj, jk, jn 80 INTEGER :: ierr 81 REAL(wp) :: zu, z0u, zzwx, zw 82 REAL(wp) :: zv, z0v, zzwy, z0w 83 REAL(wp) :: z tra, zbtr, zdt, zalpha! - -81 INTEGER :: ji, jj, jk, jn ! dummy loop indices 82 INTEGER :: ierr ! local integer 83 REAL(wp) :: zu, z0u, zzwx, zw ! local scalars 84 REAL(wp) :: zv, z0v, zzwy, z0w ! - - 85 REAL(wp) :: zdt, zalpha ! - - 84 86 REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy ! 3D workspace 85 87 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx , zwy ! - - 86 88 !!---------------------------------------------------------------------- 87 89 ! 88 IF( nn_timing == 1 ) CALL timing_start('tra_adv_mus cl')89 ! 90 CALL wrk_alloc( jpi, jpj, jpk,zslpx, zslpy, zwx, zwy )90 IF( nn_timing == 1 ) CALL timing_start('tra_adv_mus') 91 ! 92 CALL wrk_alloc( jpi,jpj,jpk, zslpx, zslpy, zwx, zwy ) 91 93 ! 92 94 IF( kt == kit000 ) THEN … … 97 99 IF(lwp) WRITE(numout,*) 98 100 ! 99 ! 100 IF( ld_msc_ups ) THEN 101 IF( .NOT. ALLOCATED( upsmsk ) ) THEN 102 ALLOCATE( upsmsk(jpi,jpj), STAT=ierr ) 103 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate upsmsk array') 104 ENDIF 101 ! Upstream / MUSCL scheme indicator 102 ! 103 ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr ) 104 xind(:,:,:) = 1._wp ! set equal to 1 where up-stream is not needed 105 ! 106 IF( ld_msc_ups ) THEN ! define the upstream indicator (if asked) 107 ALLOCATE( upsmsk(jpi,jpj), STAT=ierr ) 105 108 upsmsk(:,:) = 0._wp ! not upstream by default 106 ENDIF 107 108 IF( .NOT. ALLOCATED( xind ) ) THEN 109 ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr ) 110 IF( ierr /= 0 ) CALL ctl_stop('STOP', 'tra_adv_muscl: unable to allocate zind array') 111 ENDIF 112 ! 113 114 ! 115 ! Upstream / MUSCL scheme indicator 116 ! ------------------------------------ 117 !!gm useless 118 xind(:,:,:) = 1._wp ! set equal to 1 where up-stream is not needed 119 !!gm 120 ! 121 IF( ld_msc_ups ) THEN 109 ! 122 110 DO jk = 1, jpkm1 123 111 xind(:,:,jk) = 1._wp & ! =>1 where up-stream is not needed 124 112 & - MAX ( rnfmsk(:,:) * rnfmsk_z(jk), & ! =>0 near runoff mouths (& closed sea outflows) 125 & upsmsk(:,:) ) * tmask(:,:,jk) ! =>0 near some straits113 & upsmsk(:,:) ) * tmask(:,:,jk) ! =>0 in some user defined area 126 114 END DO 127 115 ENDIF … … 172 160 END DO 173 161 END DO 174 END DO ! interior values175 162 END DO 163 ! 176 164 ! !-- MUSCL horizontal advective fluxes 177 165 DO jk = 1, jpkm1 ! interior values … … 196 184 END DO 197 185 END DO 198 ! ! lateral boundary conditions on zwx, zwy (changed sign) 199 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) 200 ! 201 ! Tracer flux divergence at t-point added to the general trend 202 DO jk = 1, jpkm1 186 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! lateral boundary conditions (changed sign) 187 ! 188 DO jk = 1, jpkm1 ! Tracer flux divergence at t-point added to the general trend 203 189 DO jj = 2, jpjm1 204 190 DO ji = fs_2, fs_jpim1 ! vector opt. 205 zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 206 ! horizontal advective trends 207 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 208 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 209 ! add it to the general tracer trends 210 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 191 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 192 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) & 193 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 211 194 END DO 212 195 END DO … … 226 209 ! II. Vertical advective fluxes 227 210 ! ----------------------------- 228 ! !-- first guess of the slopes 229 zwx (:,:, 1 ) = 0.e0 ; zwx (:,:,jpk) = 0.e0 ! surface & bottom boundary conditions 211 ! !-- first guess of the slopes 212 zwx(:,:, 1 ) = 0._wp ! surface & bottom boundary conditions 213 zwx(:,:,jpk) = 0._wp ! surface & bottom boundary conditions 230 214 DO jk = 2, jpkm1 ! interior values 231 215 zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) ) 232 216 END DO 233 217 234 ! 235 zslpx(:,:,1) = 0. e0! surface values236 DO jk = 2, jpkm1 218 ! !-- Slopes of tracer 219 zslpx(:,:,1) = 0._wp ! surface values 220 DO jk = 2, jpkm1 ! interior value 237 221 DO jj = 1, jpj 238 222 DO ji = 1, jpi … … 242 226 END DO 243 227 END DO 244 ! 245 DO jk = 2, jpkm1 228 ! !-- Slopes limitation 229 DO jk = 2, jpkm1 ! interior values 246 230 DO jj = 1, jpj 247 231 DO ji = 1, jpi … … 252 236 END DO 253 237 END DO 254 ! !-- vertical advective flux 255 ! ! surface values (bottom already set to zero) 256 IF( lk_vvl ) THEN ; zwx(:,:, 1 ) = 0.e0 ! variable volume 257 ELSE ; zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface 258 ENDIF 259 ! 260 DO jk = 1, jpkm1 ! interior values 238 ! !-- vertical advective flux 239 DO jk = 1, jpkm1 ! interior values 261 240 zdt = p2dt(jk) 262 241 DO jj = 2, jpjm1 263 242 DO ji = fs_2, fs_jpim1 ! vector opt. 264 zbtr = 1. / ( e1e2t(ji,jj) * fse3w(ji,jj,jk+1) )265 243 z0w = SIGN( 0.5, pwn(ji,jj,jk+1) ) 266 244 zalpha = 0.5 + z0w 267 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr245 zw = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt / ( e1e2t(ji,jj) * fse3w(ji,jj,jk+1) ) 268 246 zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1) 269 247 zzwy = ptb(ji,jj,jk ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk ) 270 zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 248 zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk) 271 249 END DO 272 250 END DO 273 251 END DO 274 252 ! ! top values (bottom already set to zero) 253 IF( lk_vvl ) THEN ! variable volume 254 zwx(:,:, 1 ) = 0._wp ! k=1 only as zwx has been multiplied by wmask 255 ELSE ! linear free surface 256 IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean) 257 DO jj = 1, jpj 258 DO ji = 1, jpi 259 zwx(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) 260 END DO 261 END DO 262 ELSE ! no cavities: only at the ocean surface 263 zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 264 ENDIF 265 ENDIF 266 ! 275 267 DO jk = 1, jpkm1 ! Compute & add the vertical advective trend 276 268 DO jj = 2, jpjm1 277 269 DO ji = fs_2, fs_jpim1 ! vector opt. 278 zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 279 ! vertical advective trends 280 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) 281 ! add it to the general tracer trends 282 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 270 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 283 271 END DO 284 272 END DO … … 291 279 END DO 292 280 ! 293 CALL wrk_dealloc( jpi, jpj, jpk,zslpx, zslpy, zwx, zwy )294 ! 295 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_mus cl')296 ! 297 END SUBROUTINE tra_adv_mus cl281 CALL wrk_dealloc( jpi,jpj,jpk, zslpx, zslpy, zwx, zwy ) 282 ! 283 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_mus') 284 ! 285 END SUBROUTINE tra_adv_mus 298 286 299 287 !!====================================================================== 300 END MODULE traadv_mus cl288 END MODULE traadv_mus -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90
r5737 r5770 102 102 IF(lwp) WRITE(numout,*) 103 103 ENDIF 104 ! 104 105 l_trd = .FALSE. 105 106 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. … … 130 131 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 131 132 !! 132 INTEGER :: ji, jj, jk, jn ! dummy loop indices133 REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk ! local scalars134 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zfu, zfc, zfd133 INTEGER :: ji, jj, jk, jn ! dummy loop indices 134 REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk ! local scalars 135 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zfu, zfc, zfd 135 136 !---------------------------------------------------------------------- 136 137 ! … … 139 140 DO jn = 1, kjpt ! tracer loop 140 141 ! ! =========== 141 zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0 142 zfd(:,:,:) = 0.0 ; zwx(:,:,:) = 0.0 143 ! 144 DO jk = 1, jpkm1 145 ! 146 !--- Computation of the ustream and downstream value of the tracer and the mask 142 zfu(:,:,:) = 0._wp ; zfc(:,:,:) = 0._wp 143 zfd(:,:,:) = 0._wp ; zwx(:,:,:) = 0._wp 144 ! 145 !!gm why not using a SHIFT instruction... 146 DO jk = 1, jpkm1 !--- Computation of the ustream and downstream value of the tracer and the mask 147 147 DO jj = 2, jpjm1 148 148 DO ji = fs_2, fs_jpim1 ! vector opt. 149 ! Upstream in the x-direction for the tracer 150 zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn) 151 ! Downstream in the x-direction for the tracer 152 zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn) 149 zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn) ! Upstream in the x-direction for the tracer 150 zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn) ! Downstream in the x-direction for the tracer 153 151 END DO 154 152 END DO … … 159 157 ! Horizontal advective fluxes 160 158 ! --------------------------- 161 !162 159 DO jk = 1, jpkm1 163 160 DO jj = 2, jpjm1 … … 380 377 ! 381 378 INTEGER :: ji, jj, jk, jn ! dummy loop indices 382 REAL(wp) :: zbtr , ztra ! local scalars383 379 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz 384 380 !!---------------------------------------------------------------------- 385 381 ! 386 CALL wrk_alloc( jpi, jpj, jpk, zwz ) 382 CALL wrk_alloc( jpi,jpj,jpk, zwz ) 383 ! 384 ! ! surface & bottom values 385 IF( lk_vvl ) zwz(:,:, 1 ) = 0._wp ! set to zero one for all 386 zwz(:,:,jpk) = 0._wp ! except at the surface in linear free surface 387 ! 387 388 ! ! =========== 388 389 DO jn = 1, kjpt ! tracer loop 389 390 ! ! =========== 390 ! 1. Bottom value : flux set to zero 391 zwz(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 392 ! 393 ! ! Surface value 394 IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! Variable volume : flux set to zero 395 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) ! Constant volume : advective flux through the surface 391 ! 392 DO jk = 2, jpkm1 !* Interior point (w-masked 2nd order centered flux) 393 DO jj = 2, jpjm1 394 DO ji = fs_2, fs_jpim1 ! vector opt. 395 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk) 396 END DO 397 END DO 398 END DO 399 IF(.NOT.lk_vvl ) THEN !* top value (only in linear free surf. as zwz is multiplied by wmask) 400 IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean) 401 DO jj = 1, jpj 402 DO ji = 1, jpi 403 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn) ! linear free surface 404 END DO 405 END DO 406 ELSE ! no ice-shelf cavities (only ocean surface) 407 zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn) 408 ENDIF 396 409 ENDIF 397 410 ! 398 DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point411 DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==! 399 412 DO jj = 2, jpjm1 400 413 DO ji = fs_2, fs_jpim1 ! vector opt. 401 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) 402 END DO 403 END DO 404 END DO 405 ! 406 DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==! 407 DO jj = 2, jpjm1 408 DO ji = fs_2, fs_jpim1 ! vector opt. 409 zbtr = 1. / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 410 ! k- vertical advective trends 411 ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) 412 ! added to the general tracer trends 413 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 414 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) & 415 & / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) ) 414 416 END DO 415 417 END DO -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_ubs.F90
r5737 r5770 16 16 USE trc_oce ! share passive tracers/Ocean variables 17 17 USE trd_oce ! trends: ocean variables 18 USE traadv_fct ! acces to routine interp_4th_cpt 18 19 USE trdtra ! trends manager: tracers 19 20 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient … … 38 39 # include "vectopt_loop_substitute.h90" 39 40 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)41 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 41 42 !! $Id$ 42 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 44 45 CONTAINS 45 46 46 SUBROUTINE tra_adv_ubs ( kt, kit000, cdtype, p2dt, pun, pvn, pwn,&47 & ptb, ptn, pta, kjpt)47 SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 48 & ptb, ptn, pta, kjpt, kn_ubs_v ) 48 49 !!---------------------------------------------------------------------- 49 50 !! *** ROUTINE tra_adv_ubs *** … … 52 53 !! and add it to the general trend of passive tracer equations. 53 54 !! 54 !! ** Method : The upstream biased scheme (UBS) is based on a 3rd order55 !! ** Method : The 3rd order Upstream Biased Scheme (UBS) is based on an 55 56 !! upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005) 56 57 !! It is only used in the horizontal direction. … … 61 62 !! where zltu is the second derivative of the before temperature field: 62 63 !! zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ] 63 !! This results in a dissipatively dominant (i.e. hyper-diffusive)64 !! This results in a dissipatively dominant (i.e. hyper-diffusive) 64 65 !! truncation error. The overall performance of the advection scheme 65 66 !! is similar to that reported in (Farrow and Stevens, 1995). 66 !! For stability reasons, the first term of the fluxes which corresponds67 !! For stability reasons, the first term of the fluxes which corresponds 67 68 !! to a second order centered scheme is evaluated using the now velocity 68 69 !! (centered in time) while the second term which is the diffusive part 69 70 !! of the scheme, is evaluated using the before velocity (forward in time). 70 71 !! Note that UBS is not positive. Do not use it on passive tracers. 71 !! On the vertical, the advection is evaluated using a TVD scheme, 72 !! as the UBS have been found to be too diffusive. 72 !! On the vertical, the advection is evaluated using a FCT scheme, 73 !! as the UBS have been found to be too diffusive. 74 !!gm !! kn_ubs_v argument (not coded for the moment) 75 !! controles whether the FCT is based on a 2nd order centrered scheme (kn_ubs_v=2) 76 !! or on a 4th order compact scheme (kn_ubs_v=4). 73 77 !! 74 78 !! ** Action : - update (pta) with the now advective tracer trends … … 81 85 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 82 86 INTEGER , INTENT(in ) :: kjpt ! number of tracers 87 INTEGER , INTENT(in ) :: kn_ubs_v ! number of tracers 83 88 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 84 89 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean transport components … … 95 100 IF( nn_timing == 1 ) CALL timing_start('tra_adv_ubs') 96 101 ! 97 CALL wrk_alloc( jpi, jpj, jpk,ztu, ztv, zltu, zltv, zti, ztw )102 CALL wrk_alloc( jpi,jpj,jpk, ztu, ztv, zltu, zltv, zti, ztw ) 98 103 ! 99 104 IF( kt == kit000 ) THEN … … 106 111 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 107 112 ! 113 zltu(:,:,jpk) = 0._wp ; zltv(:,:,jpk) = 0._wp ! Bottom value : set to zero one for all 114 ztw (:,:,jpk) = 0._wp ; zti (:,:,jpk) = 0._wp 115 IF( lk_vvl ) ztw(:,:, 1 ) = 0._wp ! surface value: set to zero only in vvl case 116 ! 108 117 ! ! =========== 109 118 DO jn = 1, kjpt ! tracer loop 110 119 ! ! =========== 111 ! 1. Bottom value : flux set to zero112 ! ----------------------------------113 zltu(:,:,jpk) = 0.e0 ; zltv(:,:,jpk) = 0.e0114 120 ! 115 DO jk = 1, jpkm1 ! Horizontal slab 116 ! 117 ! Laplacian 118 DO jj = 1, jpjm1 ! First derivative (gradient) 121 DO jk = 1, jpkm1 !== horizontal laplacian of before tracer ==! 122 DO jj = 1, jpjm1 ! First derivative (masked gradient) 119 123 DO ji = 1, fs_jpim1 ! vector opt. 120 124 zeeu = e2_e1u(ji,jj) * fse3u(ji,jj,jk) * umask(ji,jj,jk) … … 124 128 END DO 125 129 END DO 126 DO jj = 2, jpjm1 ! Second derivative (divergence)130 DO jj = 2, jpjm1 ! Second derivative (divergence) 127 131 DO ji = fs_2, fs_jpim1 ! vector opt. 128 zcoef = 1. / ( 6.* fse3t(ji,jj,jk) )132 zcoef = 1._wp / ( 6._wp * fse3t(ji,jj,jk) ) 129 133 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef 130 134 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef … … 132 136 END DO 133 137 ! 134 END DO ! End of slab138 END DO 135 139 CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn) 136 137 140 ! 138 ! Horizontal advective fluxes 139 DO jk = 1, jpkm1 ! Horizontal slab 141 DO jk = 1, jpkm1 !== Horizontal advective fluxes ==! (UBS) 140 142 DO jj = 1, jpjm1 141 143 DO ji = 1, fs_jpim1 ! vector opt. 142 ! upstream transport (x2) 143 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 144 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) ! upstream transport (x2) 144 145 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 145 146 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 146 147 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 147 ! 2nd order centered advective fluxes (x2)148 ! ! 2nd order centered advective fluxes (x2) 148 149 zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ) 149 150 zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) ) 150 ! UBS advective fluxes151 ! ! UBS advective fluxes 151 152 ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) ) 152 153 ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) ) 153 154 END DO 154 155 END DO 155 END DO ! End of slab156 157 zltu(:,:,:) = pta(:,:,:,jn) ! store pta trends158 159 DO jk = 1, jpkm1 ! Horizontal advective trends156 END DO 157 ! 158 zltu(:,:,:) = pta(:,:,:,jn) ! store the initial trends before its update 159 ! 160 DO jk = 1, jpkm1 !== add the horizontal advective trend ==! 160 161 DO jj = 2, jpjm1 161 162 DO ji = fs_2, fs_jpim1 ! vector opt. … … 166 167 END DO 167 168 ! 168 END DO ! End of slab 169 170 ! Horizontal trend used in tra_adv_ztvd subroutine 171 zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) 172 169 END DO 170 ! 171 zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:) ! Horizontal advective trend used in vertical 2nd order FCT case 172 ! ! and/or in trend diagnostic (l_trd=T) 173 173 ! 174 174 IF( l_trd ) THEN ! trend diagnostics … … 181 181 IF( jn == jp_sal ) str_adv(:) = ptr_sj( ztv(:,:,:) ) 182 182 ENDIF 183 184 ! TVD scheme for the vertical direction 185 ! ---------------------- 186 IF( l_trd ) zltv(:,:,:) = pta(:,:,:,jn) ! store pta if trend diag. 187 188 ! Bottom value : flux set to zero 189 ztw(:,:,jpk) = 0.e0 ; zti(:,:,jpk) = 0.e0 190 191 ! Surface value 192 IF( lk_vvl ) THEN ; ztw(:,:,1) = 0.e0 ! variable volume : flux set to zero 193 ELSE ; ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! free constant surface 194 ENDIF 195 ! upstream advection with initial mass fluxes & intermediate update 196 ! ------------------------------------------------------------------- 197 ! Interior value 198 DO jk = 2, jpkm1 199 DO jj = 1, jpj 200 DO ji = 1, jpi 201 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 202 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 203 ztw(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 204 END DO 205 END DO 206 END DO 207 ! update and guess with monotonic sheme 208 DO jk = 1, jpkm1 209 z2dtt = p2dt(jk) 210 DO jj = 2, jpjm1 211 DO ji = fs_2, fs_jpim1 ! vector opt. 212 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 213 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr 214 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztak 215 zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 216 END DO 217 END DO 218 END DO 219 ! 220 CALL lbc_lnk( zti, 'T', 1. ) ! Lateral boundary conditions on zti, zsi (unchanged sign) 221 222 ! antidiffusive flux : high order minus low order 223 ztw(:,:,1) = 0.e0 ! Surface value 224 DO jk = 2, jpkm1 ! Interior value 225 DO jj = 1, jpj 226 DO ji = 1, jpi 227 ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - ztw(ji,jj,jk) 228 END DO 229 END DO 230 END DO 231 ! 232 CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt ) ! monotonicity algorithm 233 234 ! final trend with corrected fluxes 235 DO jk = 1, jpkm1 183 ! 184 ! !== vertical advective trend ==! 185 ! 186 SELECT CASE( kn_ubs_v ) ! select the vertical advection scheme 187 ! 188 CASE( 2 ) ! 2nd order FCT 189 ! 190 IF( l_trd ) zltv(:,:,:) = pta(:,:,:,jn) ! store pta if trend diag. 191 ! 192 ! !* upstream advection with initial mass fluxes & intermediate update ==! 193 DO jk = 2, jpkm1 ! Interior value (w-masked) 194 DO jj = 1, jpj 195 DO ji = 1, jpi 196 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 197 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 198 ztw(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 199 END DO 200 END DO 201 END DO 202 IF(.NOT.lk_vvl ) THEN ! top ocean value (only in linear free surface as ztw has been w-masked) 203 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 204 DO jj = 1, jpj 205 DO ji = 1, jpi 206 ztw(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface 207 END DO 208 END DO 209 ELSE ! no cavities: only at the ocean surface 210 ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 211 ENDIF 212 ENDIF 213 ! 214 DO jk = 1, jpkm1 !* trend and after field with monotonic scheme 215 z2dtt = p2dt(jk) 216 DO jj = 2, jpjm1 217 DO ji = fs_2, fs_jpim1 ! vector opt. 218 ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 219 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztak 220 zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk) 221 END DO 222 END DO 223 END DO 224 CALL lbc_lnk( zti, 'T', 1. ) ! Lateral boundary conditions on zti, zsi (unchanged sign) 225 ! 226 ! !* anti-diffusive flux : high order minus low order 227 DO jk = 2, jpkm1 ! Interior value (w-masked) 228 DO jj = 1, jpj 229 DO ji = 1, jpi 230 ztw(ji,jj,jk) = ( 0.5_wp * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) & 231 & - ztw(ji,jj,jk) ) * wmask(ji,jj,jk) 232 END DO 233 END DO 234 END DO 235 ! ! top ocean value: high order == upstream ==>> zwz=0 236 IF(.NOT.lk_vvl ) ztw(:,:, 1 ) = 0._wp ! only ocean surface as interior zwz values have been w-masked 237 ! 238 CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt ) ! monotonicity algorithm 239 ! 240 CASE( 4 ) ! 4th order COMPACT 241 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! 4th order compact interpolation of T at w-point 242 DO jk = 2, jpkm1 243 DO jj = 2, jpjm1 244 DO ji = fs_2, fs_jpim1 245 ztw(ji,jj,jk) = pwn(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk) 246 END DO 247 END DO 248 END DO 249 IF(.NOT.lk_vvl ) ztw(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) !!gm ISF & 4th COMPACT doesn't work 250 ! 251 END SELECT 252 ! 253 DO jk = 1, jpkm1 ! final trend with corrected fluxes 236 254 DO jj = 2, jpjm1 237 255 DO ji = fs_2, fs_jpim1 ! vector opt. 238 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 239 ! k- vertical advective trends 240 ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) 241 ! added to the general tracer trends 242 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 243 END DO 244 END DO 245 END DO 246 247 ! Save the final vertical advective trends 248 IF( l_trd ) THEN ! vertical advective trend diagnostics 256 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 257 END DO 258 END DO 259 END DO 260 ! 261 IF( l_trd ) THEN ! vertical advective trend diagnostics 249 262 DO jk = 1, jpkm1 ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w]) 250 263 DO jj = 2, jpjm1 251 264 DO ji = fs_2, fs_jpim1 ! vector opt. 252 z btr = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )253 z_hdivn = ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) * zbtr254 zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn265 zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) & 266 & + ptn(ji,jj,jk,jn) * ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) & 267 & / ( e1e2t(ji,jj) * fse3t_n(ji,jj,jk) ) 255 268 END DO 256 269 END DO … … 261 274 END DO 262 275 ! 263 CALL wrk_dealloc( jpi, jpj, jpk,ztu, ztv, zltu, zltv, zti, ztw )276 CALL wrk_dealloc( jpi,jpj,jpk, ztu, ztv, zltu, zltv, zti, ztw ) 264 277 ! 265 278 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_ubs') … … 294 307 IF( nn_timing == 1 ) CALL timing_start('nonosc_z') 295 308 ! 296 CALL wrk_alloc( jpi, jpj, jpk,zbetup, zbetdo )309 CALL wrk_alloc( jpi,jpj,jpk, zbetup, zbetdo ) 297 310 ! 298 311 zbig = 1.e+40_wp 299 312 zrtrn = 1.e-15_wp 300 313 zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp 301 314 ! 302 315 ! Search local extrema 303 316 ! -------------------- 304 ! large negative value (-zbig) inside land317 ! ! large negative value (-zbig) inside land 305 318 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 306 319 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) ) 307 ! search maximum in neighbourhood308 DO jk = 1, jpkm1 320 ! 321 DO jk = 1, jpkm1 ! search maximum in neighbourhood 309 322 ikm1 = MAX(jk-1,1) 310 323 DO jj = 2, jpjm1 … … 316 329 END DO 317 330 END DO 318 ! large positive value (+zbig) inside land331 ! ! large positive value (+zbig) inside land 319 332 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 320 333 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) ) 321 ! search minimum in neighbourhood322 DO jk = 1, jpkm1 334 ! 335 DO jk = 1, jpkm1 ! search minimum in neighbourhood 323 336 ikm1 = MAX(jk-1,1) 324 337 DO jj = 2, jpjm1 … … 330 343 END DO 331 344 END DO 332 333 ! restore masked values to zero 345 ! ! restore masked values to zero 334 346 pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) 335 347 paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) 336 337 338 ! 2. Positive and negative part of fluxes and beta terms 339 ! ------------------------------------------------------ 340 348 ! 349 ! Positive and negative part of fluxes and beta terms 350 ! --------------------------------------------------- 341 351 DO jk = 1, jpkm1 342 352 z2dtt = p2dt(jk) … … 347 357 zneg = MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) ) 348 358 ! up & down beta terms 349 zbt = e1 t(ji,jj) *e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt359 zbt = e1e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt 350 360 zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt 351 361 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt … … 353 363 END DO 354 364 END DO 365 ! 355 366 ! monotonic flux in the k direction, i.e. pcc 356 367 ! ------------------------------------------- … … 366 377 END DO 367 378 ! 368 CALL wrk_dealloc( jpi, jpj, jpk,zbetup, zbetdo )379 CALL wrk_dealloc( jpi,jpj,jpk, zbetup, zbetdo ) 369 380 ! 370 381 IF( nn_timing == 1 ) CALL timing_stop('nonosc_z') -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf.F90
r5766 r5770 73 73 SELECT CASE ( nldf ) !* compute lateral mixing trend and add it to the general trend 74 74 ! 75 CASE ( n _lap )! laplacian: iso-level operator75 CASE ( np_lap ) ! laplacian: iso-level operator 76 76 CALL tra_ldf_lap ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsa, jpts, 1 ) 77 77 ! 78 CASE ( n _lap_i )! laplacian: standard iso-neutral operator (Madec)78 CASE ( np_lap_i ) ! laplacian: standard iso-neutral operator (Madec) 79 79 CALL tra_ldf_iso ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts, 1 ) 80 80 ! 81 CASE ( n _lap_it )! laplacian: triad iso-neutral operator (griffies)81 CASE ( np_lap_it ) ! laplacian: triad iso-neutral operator (griffies) 82 82 CALL tra_ldf_triad( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb, tsb, tsa, jpts, 1 ) 83 83 ! 84 CASE ( n _blp , n_blp_i , n_blp_it )! bilaplacian: iso-level & iso-neutral operators84 CASE ( np_blp , np_blp_i , np_blp_it ) ! bilaplacian: iso-level & iso-neutral operators 85 85 CALL tra_ldf_blp ( kt, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, tsb , tsa, jpts, nldf ) 86 86 END SELECT … … 126 126 IF( ln_traldf_blp ) ioptio = ioptio + 1 127 127 IF( ioptio > 1 ) CALL ctl_stop( 'tra_ldf_init: use ONE or NONE of the 2 lap/bilap operator type on tracer' ) 128 IF( ioptio == 0 ) nldf = n _no_ldf ! No lateral diffusion128 IF( ioptio == 0 ) nldf = np_no_ldf ! No lateral diffusion 129 129 ioptio = 0 130 130 IF( ln_traldf_lev ) ioptio = ioptio + 1 … … 137 137 IF( ln_traldf_lap ) THEN ! laplacian operator 138 138 IF ( ln_zco ) THEN ! z-coordinate 139 IF ( ln_traldf_lev ) nldf = n _lap ! iso-level = horizontal (no rotation)140 IF ( ln_traldf_hor ) nldf = n _lap ! iso-level = horizontal (no rotation)141 IF ( ln_traldf_iso ) nldf = n _lap_i ! iso-neutral: standard ( rotation)142 IF ( ln_traldf_triad ) nldf = n _lap_it ! iso-neutral: triad ( rotation)139 IF ( ln_traldf_lev ) nldf = np_lap ! iso-level = horizontal (no rotation) 140 IF ( ln_traldf_hor ) nldf = np_lap ! iso-level = horizontal (no rotation) 141 IF ( ln_traldf_iso ) nldf = np_lap_i ! iso-neutral: standard ( rotation) 142 IF ( ln_traldf_triad ) nldf = np_lap_it ! iso-neutral: triad ( rotation) 143 143 ENDIF 144 144 IF ( ln_zps ) THEN ! z-coordinate with partial step 145 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed146 IF ( ln_traldf_hor ) nldf = n _lap ! horizontal(no rotation)147 IF ( ln_traldf_iso ) nldf = n _lap_i ! iso-neutral: standard(rotation)148 IF ( ln_traldf_triad ) nldf = n _lap_it ! iso-neutral: triad(rotation)145 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 146 IF ( ln_traldf_hor ) nldf = np_lap ! horizontal (no rotation) 147 IF ( ln_traldf_iso ) nldf = np_lap_i ! iso-neutral: standard (rotation) 148 IF ( ln_traldf_triad ) nldf = np_lap_it ! iso-neutral: triad (rotation) 149 149 ENDIF 150 150 IF ( ln_sco ) THEN ! s-coordinate 151 IF ( ln_traldf_lev ) nldf = n _lap ! iso-level(no rotation)152 IF ( ln_traldf_hor ) nldf = n _lap_i ! horizontal ( rotation)153 IF ( ln_traldf_iso ) nldf = n _lap_i ! iso-neutral: standard (rotation)154 IF ( ln_traldf_triad ) nldf = n _lap_it ! iso-neutral: triad (rotation)151 IF ( ln_traldf_lev ) nldf = np_lap ! iso-level (no rotation) 152 IF ( ln_traldf_hor ) nldf = np_lap_i ! horizontal ( rotation) 153 IF ( ln_traldf_iso ) nldf = np_lap_i ! iso-neutral: standard ( rotation) 154 IF ( ln_traldf_triad ) nldf = np_lap_it ! iso-neutral: triad ( rotation) 155 155 ENDIF 156 156 ENDIF … … 158 158 IF( ln_traldf_blp ) THEN ! bilaplacian operator 159 159 IF ( ln_zco ) THEN ! z-coordinate 160 IF ( ln_traldf_lev ) nldf = n _blp ! iso-level = horizontal (no rotation)161 IF ( ln_traldf_hor ) nldf = n _blp ! iso-level = horizontal (no rotation)162 IF ( ln_traldf_iso ) nldf = n _blp_i ! iso-neutral: standard (rotation)163 IF ( ln_traldf_triad ) nldf = n _blp_it ! iso-neutral: triad (rotation)160 IF ( ln_traldf_lev ) nldf = np_blp ! iso-level = horizontal (no rotation) 161 IF ( ln_traldf_hor ) nldf = np_blp ! iso-level = horizontal (no rotation) 162 IF ( ln_traldf_iso ) nldf = np_blp_i ! iso-neutral: standard ( rotation) 163 IF ( ln_traldf_triad ) nldf = np_blp_it ! iso-neutral: triad ( rotation) 164 164 ENDIF 165 165 IF ( ln_zps ) THEN ! z-coordinate with partial step 166 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed167 IF ( ln_traldf_hor ) nldf = n _blp ! horizontal(no rotation)168 IF ( ln_traldf_iso ) nldf = n _blp_i ! iso-neutral: standard (rotation)169 IF ( ln_traldf_triad ) nldf = n _blp_it ! iso-neutral: triad (rotation)166 IF ( ln_traldf_lev ) ierr = 1 ! iso-level not allowed 167 IF ( ln_traldf_hor ) nldf = np_blp ! horizontal (no rotation) 168 IF ( ln_traldf_iso ) nldf = np_blp_i ! iso-neutral: standard ( rotation) 169 IF ( ln_traldf_triad ) nldf = np_blp_it ! iso-neutral: triad ( rotation) 170 170 ENDIF 171 171 IF ( ln_sco ) THEN ! s-coordinate 172 IF ( ln_traldf_lev ) nldf = n _blp ! iso-level(no rotation)173 IF ( ln_traldf_hor ) nldf = n _blp_it ! horizontal ( rotation) !!gm a checker....174 IF ( ln_traldf_iso ) nldf = n _blp_i ! iso-neutral: standard (rotation)175 IF ( ln_traldf_triad ) nldf = n _blp_it ! iso-neutral: triad (rotation)172 IF ( ln_traldf_lev ) nldf = np_blp ! iso-level (no rotation) 173 IF ( ln_traldf_hor ) nldf = np_blp_it ! horizontal ( rotation) 174 IF ( ln_traldf_iso ) nldf = np_blp_i ! iso-neutral: standard ( rotation) 175 IF ( ln_traldf_triad ) nldf = np_blp_it ! iso-neutral: triad ( rotation) 176 176 ENDIF 177 177 ENDIF … … 181 181 & CALL ctl_stop( ' eddy induced velocity on tracers requires isopycnal', & 182 182 & ' laplacian diffusion' ) 183 IF( nldf == n _lap_i .OR. nldf == n_lap_it .OR. &184 & nldf == n _blp_i .OR. nldf == n_blp_it ) l_ldfslp = .TRUE. ! slope of neutral surfaces required183 IF( nldf == np_lap_i .OR. nldf == np_lap_it .OR. & 184 & nldf == np_blp_i .OR. nldf == np_blp_it ) l_ldfslp = .TRUE. ! slope of neutral surfaces required 185 185 ! 186 186 IF(lwp) THEN 187 187 WRITE(numout,*) 188 IF( nldf == n _no_ldf ) WRITE(numout,*) ' NO lateral diffusion'189 IF( nldf == n _lap ) WRITE(numout,*) ' laplacian iso-level operator'190 IF( nldf == n _lap_i ) WRITE(numout,*) ' Rotated laplacian operator (standard)'191 IF( nldf == n _lap_it ) WRITE(numout,*) ' Rotated laplacian operator (triad)'192 IF( nldf == n _blp ) WRITE(numout,*) ' bilaplacian iso-level operator'193 IF( nldf == n _blp_i ) WRITE(numout,*) ' Rotated bilaplacian operator (standard)'194 IF( nldf == n _blp_it ) WRITE(numout,*) ' Rotated bilaplacian operator (triad)'188 IF( nldf == np_no_ldf ) WRITE(numout,*) ' NO lateral diffusion' 189 IF( nldf == np_lap ) WRITE(numout,*) ' laplacian iso-level operator' 190 IF( nldf == np_lap_i ) WRITE(numout,*) ' Rotated laplacian operator (standard)' 191 IF( nldf == np_lap_it ) WRITE(numout,*) ' Rotated laplacian operator (triad)' 192 IF( nldf == np_blp ) WRITE(numout,*) ' bilaplacian iso-level operator' 193 IF( nldf == np_blp_i ) WRITE(numout,*) ' Rotated bilaplacian operator (standard)' 194 IF( nldf == np_blp_it ) WRITE(numout,*) ' Rotated bilaplacian operator (triad)' 195 195 ENDIF 196 196 ! -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_blp.F90
r5758 r5770 8 8 9 9 !!---------------------------------------------------------------------- 10 !! traldf_blp : update the tracer trend with the bilaplacian lateral mixing trend10 !! traldf_blp : update the tracer trend with the bilaplacian lateral mixing trend 11 11 !!---------------------------------------------------------------------- 12 USE oce ! ocean dynamics and active tracers13 USE dom_oce ! ocean space and time domain14 USE phycst ! physical constants15 USE trc_oce ! share passive tracers/Ocean variables16 USE zdf_oce ! ocean vertical physics17 USE ldftra ! lateral physics: eddy diffusivity18 USE ldfslp ! lateral physics: iso-neutral slopes19 USE traldf_lap ! lateral diffusion (Standard operator) (tra_ldf_lap routine)20 USE traldf_iso ! lateral diffusion (Standard operator) (tra_ldf_iso routine)21 USE traldf_triad ! lateral diffusion (Standard operator) (tra_ldf_triad routine)22 USE diaptr ! poleward transport diagnostics23 USE zpshde ! partial step: hor. derivative (zps_hde routine)12 USE oce ! ocean dynamics and active tracers 13 USE dom_oce ! ocean space and time domain 14 USE phycst ! physical constants 15 USE trc_oce ! share passive tracers/Ocean variables 16 USE zdf_oce ! ocean vertical physics 17 USE ldftra ! lateral physics: eddy diffusivity 18 USE ldfslp ! lateral physics: iso-neutral slopes 19 USE traldf_lap ! lateral diffusion (Standard operator) (tra_ldf_lap routine) 20 USE traldf_iso ! lateral diffusion (Standard operator) (tra_ldf_iso routine) 21 USE traldf_triad ! lateral diffusion (Standard operator) (tra_ldf_triad routine) 22 USE diaptr ! poleward transport diagnostics 23 USE zpshde ! partial step: hor. derivative (zps_hde routine) 24 24 ! 25 25 USE in_out_manager ! I/O manager … … 33 33 PRIVATE 34 34 35 PUBLIC tra_ldf_blp 35 PUBLIC tra_ldf_blp ! routine called by traldf.F90 36 36 37 ! 38 INTEGER, PARAMETER, PUBLIC :: n _no_ldf = 00! without operator (i.e. no lateral diffusive trend)39 ! !! laplacian ! bilaplacian!40 INTEGER, PARAMETER, PUBLIC :: n _lap = 10 , n_blp = 20 ! iso-level operator41 INTEGER, PARAMETER, PUBLIC :: n _lap_i = 11 , n_blp_i = 21 ! standard iso-neutral or geopotential operator42 INTEGER, PARAMETER, PUBLIC :: n _lap_it = 12 , n_blp_it = 22 ! triad iso-neutral or geopotential operator37 ! ! Flag to control the type of lateral diffusive operator 38 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 ! without operator (i.e. no lateral diffusive trend) 39 ! !! laplacian ! bilaplacian ! 40 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 ! iso-level operator 41 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 , np_blp_i = 21 ! standard iso-neutral or geopotential operator 42 INTEGER, PARAMETER, PUBLIC :: np_lap_it = 12 , np_blp_it = 22 ! triad iso-neutral or geopotential operator 43 43 44 44 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zdkt3d !: vertical tracer gradient at 2 levels … … 94 94 WRITE(numout,*) 95 95 SELECT CASE ( kldf ) 96 CASE ( n _blp ) ; WRITE(numout,*) 'tra_ldf_blp : iso-level bilaplacian operator on ', cdtype97 CASE ( n _blp_i ) ; WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)'98 CASE ( n _blp_it ) ; WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)'96 CASE ( np_blp ) ; WRITE(numout,*) 'tra_ldf_blp : iso-level bilaplacian operator on ', cdtype 97 CASE ( np_blp_i ) ; WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (Standard)' 98 CASE ( np_blp_it ) ; WRITE(numout,*) 'tra_ldf_blp : iso-neutral bilaplacian operator on ', cdtype, ' (triad)' 99 99 END SELECT 100 100 WRITE(numout,*) '~~~~~~~~~~~' … … 105 105 SELECT CASE ( kldf ) !== 1st laplacian applied to ptb (output in zlap) ==! 106 106 ! 107 CASE ( n _blp )! iso-level bilaplacian107 CASE ( np_blp ) ! iso-level bilaplacian 108 108 CALL tra_ldf_lap ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, zlap, kjpt, 1 ) 109 109 ! 110 CASE ( n _blp_i )! rotated bilaplacian : standard operator (Madec)110 CASE ( np_blp_i ) ! rotated bilaplacian : standard operator (Madec) 111 111 CALL tra_ldf_iso ( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 ) 112 112 ! 113 CASE ( n _blp_it )! rotated bilaplacian : triad operator (griffies)113 CASE ( np_blp_it ) ! rotated bilaplacian : triad operator (griffies) 114 114 CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, pgu, pgv, pgui, pgvi, ptb, ptb, zlap, kjpt, 1 ) 115 115 ! … … 126 126 SELECT CASE ( kldf ) !== 2nd laplacian applied to zlap (output in pta) ==! 127 127 ! 128 CASE ( n _blp )! iso-level bilaplacian128 CASE ( np_blp ) ! iso-level bilaplacian 129 129 CALL tra_ldf_lap ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, pta, kjpt, 2 ) 130 130 ! 131 CASE ( n _blp_i )! rotated bilaplacian : standard operator (Madec)131 CASE ( np_blp_i ) ! rotated bilaplacian : standard operator (Madec) 132 132 CALL tra_ldf_iso ( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 ) 133 133 ! 134 CASE ( n _blp_it )! rotated bilaplacian : triad operator (griffies)134 CASE ( np_blp_it ) ! rotated bilaplacian : triad operator (griffies) 135 135 CALL tra_ldf_triad( kt, kit000, cdtype, pahu, pahv, zglu, zglv, zgui, zgvi, zlap, ptb, pta, kjpt, 2 ) 136 136 ! -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf.F90
r5758 r5770 83 83 CASE ( 0 ) ; CALL tra_zdf_exp( kt, nit000, 'TRA', r2dtra, nn_zdfexp, tsb, tsa, jpts ) ! explicit scheme 84 84 CASE ( 1 ) ; CALL tra_zdf_imp( kt, nit000, 'TRA', r2dtra, tsb, tsa, jpts ) ! implicit scheme 85 CASE ( -1 ) ! esopa: test all possibility with control print86 CALL tra_zdf_exp( kt, nit000, 'TRA', r2dtra, nn_zdfexp, tsb, tsa, jpts )87 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf0 - Ta: ', mask1=tmask, &88 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )89 CALL tra_zdf_imp( kt, nit000, 'TRA', r2dtra, tsb, tsa, jpts )90 CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' zdf1 - Ta: ', mask1=tmask, &91 & tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )92 85 END SELECT 93 86 !!gm WHY here ! and I don't like that ! … … 149 142 & ' GLS or TKE scheme, the implicit scheme is required, set ln_zdfexp = .false.' ) 150 143 151 ! Test: esopa152 IF( lk_esopa ) nzdf = -1 ! All schemes used153 154 144 IF(lwp) THEN 155 145 WRITE(numout,*) 156 146 WRITE(numout,*) 'tra_zdf_init : vertical tracer physics scheme' 157 147 WRITE(numout,*) '~~~~~~~~~~~' 158 IF( nzdf == -1 ) WRITE(numout,*) ' ESOPA test All scheme used'159 148 IF( nzdf == 0 ) WRITE(numout,*) ' Explicit time-splitting scheme' 160 149 IF( nzdf == 1 ) WRITE(numout,*) ' Implicit (euler backward) scheme' -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/TRD/trdmxl.F90
r5758 r5770 800 800 END IF 801 801 802 IF( nn_cla == 1 ) CALL ctl_warn( ' You set n_cla = 1. Note that the Mixed-Layer diagnostics ', &803 & ' are not exact along the corresponding straits. ')804 805 802 ! ! allocate trdmxl arrays 806 803 IF( trd_mxl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'trd_mxl_init : unable to allocate trdmxl arrays' ) … … 809 806 810 807 811 812 nkstp = nit000 - 1 ! current time step indicator initialization 808 nkstp = nit000 - 1 ! current time step indicator initialization 813 809 814 810 -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r5758 r5770 33 33 !! - ! 2013-06 (I. Epicoco, S. Mocavero, CMCC) nemo_northcomms: setup avoiding MPI communication 34 34 !! 3.7 ! 2014-12 (G. Madec) suppression of cross land advection option 35 !! - ! 2014-12 (G. Madec) remove KPP scheme 35 !! - ! 2014-12 (G. Madec) remove KPP scheme and cross-land advection (cla) 36 36 !!---------------------------------------------------------------------- 37 37 … … 46 46 !!---------------------------------------------------------------------- 47 47 USE step_oce ! module used in the ocean time stepping module 48 USE cla ! cross land advection (tra_cla routine)49 48 USE domcfg ! domain configuration (dom_cfg routine) 50 49 USE mppini ! shared/distributed memory setting (mpp_init routine) … … 458 457 459 458 ! ! Misc. options 460 IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 ) CALL cla_init ! Cross Land Advection461 459 CALL sto_par_init ! Stochastic parametrization 462 460 IF( ln_sto_eos ) CALL sto_pts_init ! RRandom T/S fluctuations -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/step.F90
r5758 r5770 109 109 IF( lk_bdy ) CALL bdy_dta ( kstp, time_offset=+1 ) ! update dynamic & tracer data at open boundaries 110 110 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) 111 CALL FLUSH ( numout ) 111 112 112 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 113 113 ! Update stochastic parameters and random T/S fluctuations 114 114 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 115 115 CALL sto_par( kstp ) ! Stochastic parameters 116 CALL FLUSH ( numout )117 116 118 117 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 124 123 CALL bn2 ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency 125 124 CALL bn2 ( tsn, rab_n, rn2 ) ! now Brunt-Vaisala frequency 126 CALL FLUSH ( numout ) 125 127 126 ! 128 127 ! VERTICAL PHYSICS … … 137 136 avmv(:,:,:) = rn_avm0 * wvmask(:,:,:) 138 137 ENDIF 139 CALL FLUSH ( numout ) 138 140 139 IF( ln_rnf_mouth ) THEN ! increase diffusivity at rivers mouths 141 140 DO jk = 2, nkrnf ; avt(:,:,jk) = avt(:,:,jk) + 2._wp * rn_avt_rnf * rnfmsk(:,:) * tmask(:,:,jk) ; END DO … … 152 151 IF( lrst_oce .AND. lk_zdftke ) CALL tke_rst( kstp, 'WRITE' ) 153 152 IF( lrst_oce .AND. lk_zdfgls ) CALL gls_rst( kstp, 'WRITE' ) 154 CALL FLUSH ( numout )155 153 ! 156 154 ! LATERAL PHYSICS … … 179 177 ! ! eddy diffusivity coeff. and/or eiv coeff. 180 178 IF( l_ldftra_time .OR. l_ldfeiv_time ) CALL ldf_tra( kstp ) 181 write(*,*) 'after ldf_slp'182 CALL FLUSH ( numout )183 179 184 180 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 188 184 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 189 185 CALL wzv ( kstp ) ! now cross-level velocity 190 write(*,*) 'after wzv'191 CALL FLUSH ( numout )192 186 193 187 IF( lk_dynspg_ts ) THEN … … 231 225 CALL wzv ( kstp ) ! now cross-level velocity 232 226 ENDIF 233 write(*,*) 'after wzv 2'234 CALL FLUSH ( numout )235 227 236 228 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 253 245 CALL trc_stp( kstp ) ! time-stepping 254 246 #endif 255 write(*,*) 'end dyn '256 CALL FLUSH ( numout )257 247 258 248 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 270 260 IF( lk_bdy ) CALL bdy_tra_dmp( kstp ) ! bdy damping trends 271 261 CALL tra_adv ( kstp ) ! horizontal & vertical advection 272 write(*,*) 'before tra_ldf'273 CALL FLUSH ( numout )274 262 CALL tra_ldf ( kstp ) ! lateral mixing 275 write(*,*) 'after tra_ldf'276 CALL FLUSH ( numout )277 263 278 264 !!gm : why CALL to dia_ptr has been moved here??? (use trends info?) … … 316 302 CALL tra_nxt( kstp ) ! tracer fields at next time step 317 303 ENDIF 318 write(*,*) 'after tra_nxt'319 CALL FLUSH ( numout )320 304 321 305 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 360 344 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 361 345 ! 362 write(*,*) 'after dom_vvl'363 CALL FLUSH ( numout )364 365 346 366 347 !!gm : This does not only concern the dynamics ==>>> add a new title -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv.F90
r5766 r5770 4 4 !! Ocean passive tracers: advection trend 5 5 !!============================================================================== 6 !! History : 2.0 ! 05-11 (G. Madec) Original code 7 !! 3.0 ! 10-06 (C. Ethe) Adapted to passive tracers 6 !! History : 2.0 ! 2005-11 (G. Madec) Original code 7 !! 3.0 ! 2010-06 (C. Ethe) Adapted to passive tracers 8 !! 3.7 ! 2014-05 (G. Madec, C. Ethe) Add 2nd/4th order cases for CEN and FCT schemes 8 9 !!---------------------------------------------------------------------- 9 10 #if defined key_top … … 11 12 !! 'key_top' TOP models 12 13 !!---------------------------------------------------------------------- 13 !! trc_adv : compute ocean tracer advection trend 14 !! trc_adv_ini : 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 traadv_cen2 ! 2nd order centered scheme (tra_adv_cen2 routine) 19 USE traadv_tvd ! TVD scheme (tra_adv_tvd routine) 20 USE traadv_muscl ! MUSCL scheme (tra_adv_muscl routine) 21 USE traadv_muscl2 ! MUSCL2 scheme (tra_adv_muscl2 routine) 22 USE traadv_ubs ! UBS scheme (tra_adv_ubs routine) 23 USE traadv_qck ! QUICKEST scheme (tra_adv_qck routine) 24 USE traadv_mle ! ML eddy induced velocity (tra_adv_mle routine) 25 USE ldftra ! lateral diffusion coefficient on tracers 26 USE prtctl_trc ! Print control 14 !! trc_adv : compute ocean tracer advection trend 15 !! trc_adv_ini : control the different options of advection scheme 16 !!---------------------------------------------------------------------- 17 USE oce_trc ! ocean dynamics and active tracers 18 USE trc ! ocean passive tracers variables 19 USE traadv_cen ! centered scheme (tra_adv_cen routine) 20 USE traadv_fct ! FCT scheme (tra_adv_fct routine) 21 USE traadv_mus ! MUSCL scheme (tra_adv_mus routine) 22 USE traadv_ubs ! UBS scheme (tra_adv_ubs routine) 23 USE traadv_qck ! QUICKEST scheme (tra_adv_qck routine) 24 USE traadv_mle ! ML eddy induced velocity (tra_adv_mle routine) 25 USE ldftra ! lateral diffusion coefficient on tracers 26 USE ldfslp ! Lateral diffusion: slopes of neutral surfaces 27 ! 28 USE prtctl_trc ! Print control 27 29 28 30 IMPLICIT NONE … … 33 35 PUBLIC trc_adv_ini 34 36 35 ! !!: ** Advection (namtrc_adv) ** 36 LOGICAL , PUBLIC :: ln_trcadv_cen2 ! 2nd order centered scheme flag 37 LOGICAL , PUBLIC :: ln_trcadv_tvd ! TVD scheme flag 38 LOGICAL , PUBLIC :: ln_trcadv_muscl ! MUSCL scheme flag 39 LOGICAL , PUBLIC :: ln_trcadv_muscl2 ! MUSCL2 scheme flag 40 LOGICAL , PUBLIC :: ln_trcadv_ubs ! UBS scheme flag 41 LOGICAL , PUBLIC :: ln_trcadv_qck ! QUICKEST scheme flag 42 LOGICAL , PUBLIC :: ln_trcadv_msc_ups ! use upstream scheme within muscl 43 44 INTEGER :: nadv ! choice of the type of advection scheme 37 ! !!* Namelist namtrc_adv * 38 LOGICAL :: ln_trcadv_cen ! centered scheme flag 39 INTEGER :: nn_cen_h, nn_cen_v ! =2/4 : horizontal and vertical choices of the order of CEN scheme 40 LOGICAL :: ln_trcadv_fct ! FCT scheme flag 41 INTEGER :: nn_fct_h, nn_fct_v ! =2/4 : horizontal and vertical choices of the order of FCT scheme 42 INTEGER :: nn_fct_zts ! >=1 : 2nd order FCT with vertical sub-timestepping 43 LOGICAL :: ln_trcadv_mus ! MUSCL scheme flag 44 LOGICAL :: ln_mus_ups ! use upstream scheme in vivcinity of river mouths 45 LOGICAL :: ln_trcadv_ubs ! UBS scheme flag 46 INTEGER :: nn_ubs_v ! =2/4 : vertical choice of the order of UBS scheme 47 LOGICAL :: ln_trcadv_qck ! QUICKEST scheme flag 48 49 ! ! choices of advection scheme: 50 INTEGER, PARAMETER :: np_NO_adv = 0 ! no T-S advection 51 INTEGER, PARAMETER :: np_CEN = 1 ! 2nd/4th order centered scheme 52 INTEGER, PARAMETER :: np_FCT = 2 ! 2nd/4th order Flux Corrected Transport scheme 53 INTEGER, PARAMETER :: np_FCT_zts = 3 ! 2nd order FCT scheme with vertical sub-timestepping 54 INTEGER, PARAMETER :: np_MUS = 4 ! MUSCL scheme 55 INTEGER, PARAMETER :: np_UBS = 5 ! 3rd order Upstream Biased Scheme 56 INTEGER, PARAMETER :: np_QCK = 6 ! QUICK scheme 57 58 INTEGER :: nadv ! chosen advection scheme 59 ! 45 60 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: r2dt ! vertical profile time-step, = 2 rdttra 46 61 ! ! except at nitrrc000 (=rdttra) if neuler=0 … … 50 65 # include "vectopt_loop_substitute.h90" 51 66 !!---------------------------------------------------------------------- 52 !! NEMO/TOP 3. 3 , NEMO Consortium (2010)67 !! NEMO/TOP 3.7 , NEMO Consortium (2015) 53 68 !! $Id$ 54 69 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 60 75 !! *** ROUTINE trc_adv_alloc *** 61 76 !!---------------------------------------------------------------------- 62 77 ! 63 78 ALLOCATE( r2dt(jpk), STAT=trc_adv_alloc ) 64 79 ! 65 80 IF( trc_adv_alloc /= 0 ) CALL ctl_warn('trc_adv_alloc : failed to allocate array.') 66 81 ! 67 82 END FUNCTION trc_adv_alloc 68 83 … … 78 93 INTEGER, INTENT(in) :: kt ! ocean time-step index 79 94 ! 80 INTEGER :: jk 95 INTEGER :: jk ! dummy loop index 81 96 CHARACTER (len=22) :: charout 82 97 REAL(wp), POINTER, DIMENSION(:,:,:) :: zun, zvn, zwn ! effective velocity … … 93 108 ENDIF 94 109 ! !== effective transport ==! 95 zun(:,:,jpk) = 0._wp ! no transport trough the bottom96 zvn(:,:,jpk) = 0._wp97 zwn(:,:,jpk) = 0._wp98 110 DO jk = 1, jpkm1 99 111 zun(:,:,jk) = e2u (:,:) * fse3u(:,:,jk) * un(:,:,jk) ! eulerian transport … … 112 124 IF( ln_mle ) CALL tra_adv_mle( kt, nittrc000, zun, zvn, zwn, 'TRC' ) ! add the mle transport 113 125 ! 114 ! 115 SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==! 116 ! 117 CASE ( 1 ) ; CALL tra_adv_cen2 ( kt, nittrc000, 'TRC', zun, zvn, zwn, trb, trn, tra, jptra ) ! 2nd order centered 118 CASE ( 2 ) ; CALL tra_adv_tvd ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra ) ! TVD 119 CASE ( 3 ) ; CALL tra_adv_muscl ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, tra, jptra, ln_trcadv_msc_ups ) ! MUSCL 120 CASE ( 4 ) ; CALL tra_adv_muscl2( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra ) ! MUSCL2 121 CASE ( 5 ) ; CALL tra_adv_ubs ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra ) ! UBS 122 CASE ( 6 ) ; CALL tra_adv_qck ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra ) ! QUICKEST 126 zun(:,:,jpk) = 0._wp ! no transport trough the bottom 127 zvn(:,:,jpk) = 0._wp 128 zwn(:,:,jpk) = 0._wp 129 ! 130 ! 131 SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==! 132 ! 133 CASE ( np_CEN ) ! Centered : 2nd / 4th order 134 CALL tra_adv_cen ( kt, nittrc000,'TRC', zun, zvn, zwn , trn, tra, jptra, nn_cen_h, nn_cen_v ) 135 CASE ( np_FCT ) ! FCT : 2nd / 4th order 136 CALL tra_adv_fct ( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra, nn_fct_h, nn_fct_v ) 137 CASE ( np_FCT_zts ) ! 2nd order FCT with vertical time-splitting 138 CALL tra_adv_fct_zts( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra , nn_fct_zts ) 139 CASE ( np_MUS ) ! MUSCL 140 CALL tra_adv_mus ( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, tra, jptra , ln_mus_ups ) 141 CASE ( np_UBS ) ! UBS 142 CALL tra_adv_ubs ( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra , nn_ubs_v ) 143 CASE ( np_QCK ) ! QUICKEST 144 CALL tra_adv_qck ( kt, nittrc000,'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra ) 123 145 ! 124 146 END SELECT … … 146 168 INTEGER :: ios ! Local integer output status for namelist read 147 169 !! 148 NAMELIST/namtrc_adv/ ln_trcadv_cen2 , ln_trcadv_tvd , & 149 & ln_trcadv_muscl, ln_trcadv_muscl2, & 150 & ln_trcadv_ubs , ln_trcadv_qck, ln_trcadv_msc_ups 170 NAMELIST/namtrc_adv/ ln_trcadv_cen, nn_cen_h, nn_cen_v, & ! CEN 171 & ln_trcadv_fct, nn_fct_h, nn_fct_v, nn_fct_zts, & ! FCT 172 & ln_trcadv_mus, ln_mus_ups, & ! MUSCL 173 & ln_trcadv_ubs, nn_ubs_v, & ! UBS 174 & ln_trcadv_qck ! QCK 151 175 !!---------------------------------------------------------------------- 152 176 ! … … 165 189 WRITE(numout,*) '~~~~~~~~~~~' 166 190 WRITE(numout,*) ' Namelist namtrc_adv : chose a advection scheme for tracers' 167 WRITE(numout,*) ' 2nd order advection scheme ln_trcadv_cen2 = ', ln_trcadv_cen2 168 WRITE(numout,*) ' TVD advection scheme ln_trcadv_tvd = ', ln_trcadv_tvd 169 WRITE(numout,*) ' MUSCL advection scheme ln_trcadv_muscl = ', ln_trcadv_muscl 170 WRITE(numout,*) ' MUSCL2 advection scheme ln_trcadv_muscl2 = ', ln_trcadv_muscl2 171 WRITE(numout,*) ' UBS advection scheme ln_trcadv_ubs = ', ln_trcadv_ubs 172 WRITE(numout,*) ' QUICKEST advection scheme ln_trcadv_qck = ', ln_trcadv_qck 191 WRITE(numout,*) ' centered scheme ln_trcadv_cen = ', ln_trcadv_cen 192 WRITE(numout,*) ' horizontal 2nd/4th order nn_cen_h = ', nn_fct_h 193 WRITE(numout,*) ' vertical 2nd/4th order nn_cen_v = ', nn_fct_v 194 WRITE(numout,*) ' Flux Corrected Transport scheme ln_trcadv_fct = ', ln_trcadv_fct 195 WRITE(numout,*) ' horizontal 2nd/4th order nn_fct_h = ', nn_fct_h 196 WRITE(numout,*) ' vertical 2nd/4th order nn_fct_v = ', nn_fct_v 197 WRITE(numout,*) ' 2nd order + vertical sub-timestepping nn_fct_zts = ', nn_fct_zts 198 WRITE(numout,*) ' MUSCL scheme ln_trcadv_mus = ', ln_trcadv_mus 199 WRITE(numout,*) ' + upstream scheme near river mouths ln_mus_ups = ', ln_mus_ups 200 WRITE(numout,*) ' UBS scheme ln_trcadv_ubs = ', ln_trcadv_ubs 201 WRITE(numout,*) ' vertical 2nd/4th order nn_ubs_v = ', nn_ubs_v 202 WRITE(numout,*) ' QUICKEST scheme ln_trcadv_qck = ', ln_trcadv_qck 173 203 ENDIF 174 204 ! 175 205 176 206 ioptio = 0 ! Parameter control 177 IF( ln_trcadv_cen2 ) ioptio = ioptio + 1 178 IF( ln_trcadv_tvd ) ioptio = ioptio + 1 179 IF( ln_trcadv_muscl ) ioptio = ioptio + 1 180 IF( ln_trcadv_muscl2 ) ioptio = ioptio + 1 181 IF( ln_trcadv_ubs ) ioptio = ioptio + 1 182 IF( ln_trcadv_qck ) ioptio = ioptio + 1 183 ! 207 IF( ln_trcadv_cen ) nadv = np_CEN 208 IF( ln_trcadv_fct ) nadv = np_FCT 209 IF( ln_trcadv_fct .AND. nn_fct_zts > 0 ) nadv = np_FCT_zts 210 IF( ln_trcadv_mus ) nadv = np_MUS 211 IF( ln_trcadv_ubs ) nadv = np_UBS 212 IF( ln_trcadv_qck ) nadv = np_QCK 213 ! 214 IF( ioptio == 0 ) THEN 215 nadv = np_NO_adv 216 CALL ctl_warn( 'trc_adv_init: You are running without tracer advection.' ) 217 ENDIF 184 218 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE advection scheme in namelist namtrc_adv' ) 185 219 ! 220 IF( ln_trcadv_cen .AND. ( nn_cen_h /= 2 .AND. nn_cen_h /= 4 ) & 221 .AND. ( nn_cen_v /= 2 .AND. nn_cen_v /= 4 ) ) THEN 222 CALL ctl_stop( 'trc_adv_init: CEN scheme, choose 2nd or 4th order' ) 223 ENDIF 224 IF( ln_trcadv_fct .AND. ( nn_fct_h /= 2 .AND. nn_fct_h /= 4 ) & 225 .AND. ( nn_fct_v /= 2 .AND. nn_fct_v /= 4 ) ) THEN 226 CALL ctl_stop( 'trc_adv_init: FCT scheme, choose 2nd or 4th order' ) 227 ENDIF 228 IF( ln_trcadv_fct .AND. nn_fct_zts > 0 ) THEN 229 IF( nn_fct_h == 4 ) THEN 230 nn_fct_h = 2 231 CALL ctl_stop( 'trc_adv_init: force 2nd order FCT scheme, 4th order does not exist with sub-timestepping' ) 232 ENDIF 233 IF( lk_vvl ) THEN 234 CALL ctl_stop( 'trc_adv_init: vertical sub-timestepping not allow in non-linear free surface' ) 235 ENDIF 236 IF( nn_fct_zts == 1 ) CALL ctl_warn( 'trc_adv_init: FCT with ONE sub-timestep = FCT without sub-timestep' ) 237 ENDIF 238 IF( ln_trcadv_ubs .AND. ( nn_ubs_v /= 2 .AND. nn_ubs_v /= 4 ) ) THEN 239 CALL ctl_stop( 'trc_adv_init: UBS scheme, choose 2nd or 4th order' ) 240 ENDIF 241 IF( ln_trcadv_ubs .AND. nn_ubs_v == 4 ) THEN 242 CALL ctl_warn( 'trc_adv_init: UBS scheme, only 2nd FCT scheme available on the vertical. It will be used' ) 243 ENDIF 244 IF( ln_isfcav ) THEN ! ice-shelf cavities 245 IF( ln_trcadv_cen .AND. nn_cen_v /= 4 .OR. & ! NO 4th order with ISF 246 & ln_trcadv_fct .AND. nn_fct_v /= 4 ) CALL ctl_stop( 'tra_adv_init: 4th order COMPACT scheme not allowed with ISF' ) 247 ENDIF 248 ! 186 249 ! ! Set nadv 187 IF( ln_trcadv_cen 2 ) nadv = 1188 IF( ln_trcadv_ tvd ) nadv = 2189 IF( ln_trcadv_muscl ) nadv = 3190 IF( ln_trcadv_mus cl2 ) nadv = 4191 IF( ln_trcadv_ubs ) nadv = 5192 IF( ln_trcadv_qck ) nadv = 6250 IF( ln_trcadv_cen ) nadv = np_CEN 251 IF( ln_trcadv_fct ) nadv = np_FCT 252 IF( nn_fct_zts > 0 ) nadv = np_FCT_zts 253 IF( ln_trcadv_mus ) nadv = np_MUS 254 IF( ln_trcadv_ubs ) nadv = np_UBS 255 IF( ln_trcadv_qck ) nadv = np_QCK 193 256 ! 194 257 IF(lwp) THEN ! Print the choice 195 258 WRITE(numout,*) 196 IF( nadv == 1 ) WRITE(numout,*) ' 2nd order scheme is used' 197 IF( nadv == 2 ) WRITE(numout,*) ' TVD scheme is used' 198 IF( nadv == 3 ) WRITE(numout,*) ' MUSCL scheme is used' 199 IF( nadv == 4 ) WRITE(numout,*) ' MUSCL2 scheme is used' 200 IF( nadv == 5 ) WRITE(numout,*) ' UBS scheme is used' 201 IF( nadv == 6 ) WRITE(numout,*) ' QUICKEST scheme is used' 259 IF( nadv == np_NO_adv ) WRITE(numout,*) ' NO passive tracer advection' 260 IF( nadv == np_CEN ) WRITE(numout,*) ' CEN scheme is used. Horizontal order: ', nn_cen_h, & 261 & ' Vertical order: ', nn_cen_v 262 IF( nadv == np_FCT ) WRITE(numout,*) ' FCT scheme is used. Horizontal order: ', nn_fct_h, & 263 & ' Vertical order: ', nn_fct_v 264 IF( nadv == np_FCT_zts ) WRITE(numout,*) ' use 2nd order FCT with ', nn_fct_zts,'vertical sub-timestepping' 265 IF( nadv == np_MUS ) WRITE(numout,*) ' MUSCL scheme is used' 266 IF( nadv == np_UBS ) WRITE(numout,*) ' UBS scheme is used' 267 IF( nadv == np_QCK ) WRITE(numout,*) ' QUICKEST scheme is used' 202 268 ENDIF 203 269 ! -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf.F90
r5766 r5770 73 73 IF( nn_timing == 1 ) CALL timing_start('trc_ldf') 74 74 ! 75 76 75 IF( l_trdtrc ) THEN 77 76 CALL wrk_alloc( jpi,jpj,jpk,jptra, ztrtrd ) 78 77 ztrtrd(:,:,:,:) = tra(:,:,:,:) 79 78 ENDIF 80 81 ! ! set the lateral diffusivity coef. for passive tracer79 ! 80 ! !* set the lateral diffusivity coef. for passive tracer 82 81 CALL wrk_alloc( jpi,jpj,jpk, zahu, zahv ) 83 82 zahu(:,:,:) = rldf * ahtu(:,:,:) … … 86 85 SELECT CASE ( nldf ) !* compute lateral mixing trend and add it to the general trend 87 86 ! 88 CASE ( n _lap )! iso-level laplacian87 CASE ( np_lap ) ! iso-level laplacian 89 88 CALL tra_ldf_lap ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, tra, jptra, 1 ) 90 89 ! 91 CASE ( n _lap_i )! laplacian : standard iso-neutral operator (Madec)90 CASE ( np_lap_i ) ! laplacian : standard iso-neutral operator (Madec) 92 91 CALL tra_ldf_iso ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, trb, tra, jptra, 1 ) 93 92 ! 94 CASE ( n _lap_it )! laplacian : triad iso-neutral operator (griffies)93 CASE ( np_lap_it ) ! laplacian : triad iso-neutral operator (griffies) 95 94 CALL tra_ldf_triad( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb, trb, tra, jptra, 1 ) 96 95 ! 97 CASE ( n _blp , n_blp_i , n_blp_it )! bilaplacian: all operator (iso-level, -neutral)96 CASE ( np_blp , np_blp_i , np_blp_it ) ! bilaplacian: all operator (iso-level, -neutral) 98 97 CALL tra_ldf_blp ( kt, nittrc000,'TRC', zahu, zahv, gtru, gtrv, gtrui, gtrvi, trb , tra, jptra, nldf ) 99 98 ! … … 119 118 END SUBROUTINE trc_ldf 120 119 121 !!gm ldf_ctl should be called in trcini so that l_ldfslp=T cause the slope init and calculation122 120 123 121 SUBROUTINE trc_ldf_ini … … 128 126 !! 129 127 !! ** Method : set nldf from the namtra_ldf logicals 130 !! nldf == -2 No lateral diffusion131 128 !! nldf == 0 laplacian operator 132 129 !! nldf == 1 Rotated laplacian operator … … 134 131 !! nldf == 3 Rotated bilaplacian 135 132 !!---------------------------------------------------------------------- 136 INTEGER :: ioptio, ierr 137 INTEGER :: ios! Local integer output status for namelist read138 ! 139 NAMELIST/namtrc_ldf/ ln_trcldf_lap, ln_trcldf_blp, &133 INTEGER :: ioptio, ierr ! temporary integers 134 INTEGER :: ios ! Local integer output status for namelist read 135 ! 136 NAMELIST/namtrc_ldf/ ln_trcldf_lap, ln_trcldf_blp, & 140 137 & ln_trcldf_lev, ln_trcldf_hor, ln_trcldf_iso, ln_trcldf_triad, & 141 138 & rn_ahtrc_0 , rn_bhtrc_0 … … 172 169 IF( ln_trcldf_lap ) ioptio = ioptio + 1 173 170 IF( ln_trcldf_blp ) ioptio = ioptio + 1 174 IF( ioptio > 1 ) CALL ctl_stop( 'trc_ldf_ctl: use ONE or NONE of the 2 lap/bilap operator type on tracer' )175 IF( ioptio == 0 ) nldf = n_no_ldf ! No lateral diffusion171 IF( ioptio > 1 ) CALL ctl_stop( 'trc_ldf_ctl: use ONE or NONE of the 2 lap/bilap operator type on tracer' ) 172 IF( ioptio == 0 ) nldf = np_no_ldf ! No lateral diffusion 176 173 177 174 IF( ln_trcldf_lap .AND. ln_trcldf_blp ) CALL ctl_stop( 'trc_ldf_ctl: bilaplacian should be used on both TRC and TRA' ) … … 189 186 IF( ln_trcldf_lap ) THEN !== laplacian operator ==! 190 187 IF ( ln_zco ) THEN ! z-coordinate 191 IF ( ln_trcldf_lev ) nldf = n _lap ! iso-level = horizontal (no rotation)192 IF ( ln_trcldf_hor ) nldf = n _lap ! iso-level = horizontal (no rotation)193 IF ( ln_trcldf_iso ) nldf = n _lap_i ! iso-neutral: standard ( rotation)194 IF ( ln_trcldf_triad ) nldf = n _lap_it ! iso-neutral: triad ( rotation)188 IF ( ln_trcldf_lev ) nldf = np_lap ! iso-level = horizontal (no rotation) 189 IF ( ln_trcldf_hor ) nldf = np_lap ! iso-level = horizontal (no rotation) 190 IF ( ln_trcldf_iso ) nldf = np_lap_i ! iso-neutral: standard ( rotation) 191 IF ( ln_trcldf_triad ) nldf = np_lap_it ! iso-neutral: triad ( rotation) 195 192 ENDIF 196 193 IF ( ln_zps ) THEN ! z-coordinate with partial step 197 194 IF ( ln_trcldf_lev ) ierr = 1 ! iso-level not allowed 198 IF ( ln_trcldf_hor ) nldf = n _lap ! horizontal (no rotation)199 IF ( ln_trcldf_iso ) nldf = n _lap_i ! iso-neutral: standard (rotation)200 IF ( ln_trcldf_triad ) nldf = n _lap_it ! iso-neutral: triad (rotation)195 IF ( ln_trcldf_hor ) nldf = np_lap ! horizontal (no rotation) 196 IF ( ln_trcldf_iso ) nldf = np_lap_i ! iso-neutral: standard (rotation) 197 IF ( ln_trcldf_triad ) nldf = np_lap_it ! iso-neutral: triad (rotation) 201 198 ENDIF 202 199 IF ( ln_sco ) THEN ! s-coordinate 203 IF ( ln_trcldf_lev ) nldf = n _lap ! iso-level (no rotation)204 IF ( ln_trcldf_hor ) nldf = n _lap_it ! horizontal ( rotation) !!gm a checker....205 IF ( ln_trcldf_iso ) nldf = n _lap_i ! iso-neutral: standard (rotation)206 IF ( ln_trcldf_triad ) nldf = n _lap_it ! iso-neutral: triad (rotation)200 IF ( ln_trcldf_lev ) nldf = np_lap ! iso-level (no rotation) 201 IF ( ln_trcldf_hor ) nldf = np_lap_it ! horizontal ( rotation) !!gm a checker.... 202 IF ( ln_trcldf_iso ) nldf = np_lap_i ! iso-neutral: standard (rotation) 203 IF ( ln_trcldf_triad ) nldf = np_lap_it ! iso-neutral: triad (rotation) 207 204 ENDIF 208 205 ! ! diffusivity ratio: passive / active tracers … … 217 214 ENDIF 218 215 ENDIF 219 216 ! 220 217 IF( ln_trcldf_blp ) THEN !== bilaplacian operator ==! 221 218 IF ( ln_zco ) THEN ! z-coordinate 222 IF ( ln_trcldf_lev ) nldf = n _blp ! iso-level = horizontal (no rotation)223 IF ( ln_trcldf_hor ) nldf = n _blp ! iso-level = horizontal (no rotation)224 IF ( ln_trcldf_iso ) nldf = n _blp_i ! iso-neutral: standard (rotation)225 IF ( ln_trcldf_triad ) nldf = n _blp_it ! iso-neutral: triad (rotation)219 IF ( ln_trcldf_lev ) nldf = np_blp ! iso-level = horizontal (no rotation) 220 IF ( ln_trcldf_hor ) nldf = np_blp ! iso-level = horizontal (no rotation) 221 IF ( ln_trcldf_iso ) nldf = np_blp_i ! iso-neutral: standard (rotation) 222 IF ( ln_trcldf_triad ) nldf = np_blp_it ! iso-neutral: triad (rotation) 226 223 ENDIF 227 224 IF ( ln_zps ) THEN ! z-coordinate with partial step 228 225 IF ( ln_trcldf_lev ) ierr = 1 ! iso-level not allowed 229 IF ( ln_trcldf_hor ) nldf = n _blp ! horizontal (no rotation)230 IF ( ln_trcldf_iso ) nldf = n _blp_i ! iso-neutral: standard (rotation)231 IF ( ln_trcldf_triad ) nldf = n _blp_it ! iso-neutral: triad (rotation)226 IF ( ln_trcldf_hor ) nldf = np_blp ! horizontal (no rotation) 227 IF ( ln_trcldf_iso ) nldf = np_blp_i ! iso-neutral: standard (rotation) 228 IF ( ln_trcldf_triad ) nldf = np_blp_it ! iso-neutral: triad (rotation) 232 229 ENDIF 233 230 IF ( ln_sco ) THEN ! s-coordinate 234 IF ( ln_trcldf_lev ) nldf = n _blp ! iso-level (no rotation)235 IF ( ln_trcldf_hor ) nldf = n _blp_it ! horizontal ( rotation) !!gm a checker....236 IF ( ln_trcldf_iso ) nldf = n _blp_i ! iso-neutral: standard (rotation)237 IF ( ln_trcldf_triad ) nldf = n _blp_it ! iso-neutral: triad (rotation)231 IF ( ln_trcldf_lev ) nldf = np_blp ! iso-level (no rotation) 232 IF ( ln_trcldf_hor ) nldf = np_blp_it ! horizontal ( rotation) !!gm a checker.... 233 IF ( ln_trcldf_iso ) nldf = np_blp_i ! iso-neutral: standard (rotation) 234 IF ( ln_trcldf_triad ) nldf = np_blp_it ! iso-neutral: triad (rotation) 238 235 ENDIF 239 236 ! ! diffusivity ratio: passive / active tracers … … 248 245 ENDIF 249 246 ENDIF 250 247 ! 251 248 IF( ierr == 1 ) CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' ) 252 249 IF( ln_ldfeiv .AND. .NOT.ln_trcldf_iso ) & … … 256 253 IF( .NOT.l_ldfslp ) CALL ctl_stop( ' the rotation of the diffusive tensor require l_ldfslp' ) 257 254 ENDIF 258 255 ! 259 256 IF(lwp) THEN 260 257 WRITE(numout,*) 261 IF( nldf == n _no_ldf ) WRITE(numout,*) ' NO lateral diffusion'262 IF( nldf == n _lap ) WRITE(numout,*) ' laplacian iso-level operator'263 IF( nldf == n _lap_i ) WRITE(numout,*) ' Rotated laplacian operator (standard)'264 IF( nldf == n _lap_it ) WRITE(numout,*) ' Rotated laplacian operator (triad)'265 IF( nldf == n _blp ) WRITE(numout,*) ' bilaplacian iso-level operator'266 IF( nldf == n _blp_i ) WRITE(numout,*) ' Rotated bilaplacian operator (standard)'267 IF( nldf == n _blp_it ) WRITE(numout,*) ' Rotated bilaplacian operator (triad)'258 IF( nldf == np_no_ldf ) WRITE(numout,*) ' NO lateral diffusion' 259 IF( nldf == np_lap ) WRITE(numout,*) ' laplacian iso-level operator' 260 IF( nldf == np_lap_i ) WRITE(numout,*) ' Rotated laplacian operator (standard)' 261 IF( nldf == np_lap_it ) WRITE(numout,*) ' Rotated laplacian operator (triad)' 262 IF( nldf == np_blp ) WRITE(numout,*) ' bilaplacian iso-level operator' 263 IF( nldf == np_blp_i ) WRITE(numout,*) ' Rotated bilaplacian operator (standard)' 264 IF( nldf == np_blp_it ) WRITE(numout,*) ' Rotated bilaplacian operator (triad)' 268 265 ENDIF 269 266 ! -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf.F90
r5766 r5770 11 11 !! 'key_top' TOP models 12 12 !!---------------------------------------------------------------------- 13 !! trc_zdf : update the tracer trend with the lateral diffusion14 !! trc_zdf_ini : initialization, namelist read, and parameters control13 !! trc_zdf : update the tracer trend with the lateral diffusion 14 !! trc_zdf_ini : initialization, namelist read, and parameters control 15 15 !!---------------------------------------------------------------------- 16 USE oce_trc ! ocean dynamics and active tracers17 USE trc ! ocean passive tracers variables18 USE tr azdf_exp ! vertical diffusion: explicit (tra_zdf_exp routine)19 USE trazdf_ imp ! vertical diffusion: implicit (tra_zdf_imp routine)20 USE tr cldf21 USE tr d_oce22 USE trdtra 23 USE prtctl_trc 16 USE trc ! ocean passive tracers variables 17 USE oce_trc ! ocean dynamics and active tracers 18 USE trd_oce ! trends: ocean variables 19 USE trazdf_exp ! vertical diffusion: explicit (tra_zdf_exp routine) 20 USE trazdf_imp ! vertical diffusion: implicit (tra_zdf_imp routine) 21 USE trcldf ! passive tracers: lateral diffusion 22 USE trdtra ! trends manager: tracers 23 USE prtctl_trc ! Print control 24 24 25 25 IMPLICIT NONE 26 26 PRIVATE 27 27 28 PUBLIC trc_zdf 29 PUBLIC trc_zdf_ alloc! called by nemogcm.F9030 PUBLIC trc_zdf_ ini! called by nemogcm.F9031 ! !!: ** Vertical diffusion32 ! (nam_trczdf) **28 PUBLIC trc_zdf ! called by step.F90 29 PUBLIC trc_zdf_ini ! called by nemogcm.F90 30 PUBLIC trc_zdf_alloc ! called by nemogcm.F90 31 32 ! !!** Vertical diffusion (nam_trczdf) ** 33 33 LOGICAL , PUBLIC :: ln_trczdf_exp !: explicit vertical diffusion scheme flag 34 34 INTEGER , PUBLIC :: nn_trczdf_exp !: number of sub-time step (explicit time stepping) … … 44 44 # include "vectopt_loop_substitute.h90" 45 45 !!---------------------------------------------------------------------- 46 !! NEMO/TOP 3. 3 , NEMO Consortium (2010)46 !! NEMO/TOP 3.7 , NEMO Consortium (2015) 47 47 !! $Id$ 48 48 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 59 59 ! 60 60 END FUNCTION trc_zdf_alloc 61 61 62 62 63 SUBROUTINE trc_zdf( kt ) … … 87 88 88 89 SELECT CASE ( nzdf ) ! compute lateral mixing trend and add it to the general trend 89 CASE ( -1 ) ! esopa: test all possibility with control print90 CALL tra_zdf_exp( kt, nittrc000, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra )91 WRITE(charout, FMT="('zdf1 ')") ; CALL prt_ctl_trc_info(charout)92 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )93 CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dt, trb, tra, jptra )94 WRITE(charout, FMT="('zdf2 ')") ; CALL prt_ctl_trc_info(charout)95 CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )96 90 CASE ( 0 ) ; CALL tra_zdf_exp( kt, nittrc000, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra ) ! explicit scheme 97 91 CASE ( 1 ) ; CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dt, trb, tra, jptra ) ! implicit scheme 98 99 92 END SELECT 100 93 … … 118 111 END SUBROUTINE trc_zdf 119 112 113 120 114 SUBROUTINE trc_zdf_ini 121 115 !!---------------------------------------------------------------------- … … 132 126 !!---------------------------------------------------------------------- 133 127 INTEGER :: ios ! Local integer output status for namelist read 134 ! 128 !! 135 129 NAMELIST/namtrc_zdf/ ln_trczdf_exp , nn_trczdf_exp 136 130 !!---------------------------------------------------------------------- 137 138 ! ! Vertical mixing 131 ! 139 132 REWIND( numnat_ref ) ! namtrc_zdf in reference namelist 140 133 READ ( numnat_ref, namtrc_zdf, IOSTAT = ios, ERR = 905) 141 134 905 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_zdf in reference namelist', lwp ) 142 135 ! 143 136 REWIND( numnat_cfg ) ! namtrc_zdf in configuration namelist 144 137 READ ( numnat_cfg, namtrc_zdf, IOSTAT = ios, ERR = 906 ) 145 138 906 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namtrc_zdf in configuration namelist', lwp ) 146 139 IF(lwm) WRITE ( numont, namtrc_zdf ) 147 148 IF(lwp) THEN ! !Control print140 ! 141 IF(lwp) THEN ! Control print 149 142 WRITE(numout,*) 150 143 WRITE(numout,*) ' Namelist namtrc_zdf : set vertical diffusion parameters' … … 153 146 ENDIF 154 147 155 ! Define the vertical tracer physics scheme156 IF( ln_trczdf_exp ) THEN ;nzdf = 0 ! explicit scheme157 ELSE ;nzdf = 1 ! implicit scheme148 ! ! Define the vertical tracer physics scheme 149 IF( ln_trczdf_exp ) THEN ; nzdf = 0 ! explicit scheme 150 ELSE ; nzdf = 1 ! implicit scheme 158 151 ENDIF 159 152 160 ! Force implicit schemes161 IF( ln_trcldf_iso 162 IF( ln_trcldf_hor .AND. ln_sco 153 ! ! Force implicit schemes 154 IF( ln_trcldf_iso ) nzdf = 1 ! iso-neutral lateral physics 155 IF( ln_trcldf_hor .AND. ln_sco ) nzdf = 1 ! horizontal lateral physics in s-coordinate 163 156 #if defined key_zdftke || defined key_zdfgls 164 157 nzdf = 1 ! TKE or GLS physics 165 158 #endif 166 159 IF( ln_trczdf_exp .AND. nzdf == 1 ) & … … 175 168 IF( nzdf == 1 ) WRITE(numout,*) ' Implicit (euler backward) scheme' 176 169 ENDIF 177 170 ! 178 171 END SUBROUTINE trc_zdf_ini 172 179 173 #else 180 174 !!---------------------------------------------------------------------- -
branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/TOP_SRC/trcini.F90
r5766 r5770 53 53 !! or read data or analytical formulation 54 54 !!--------------------------------------------------------------------- 55 !!---------------------------------------------------------------------56 55 ! 57 56 IF( nn_timing == 1 ) CALL timing_start('trc_init') … … 92 91 ! 93 92 END SUBROUTINE trc_init 93 94 94 95 95 SUBROUTINE trc_ini_ctl … … 103 103 l_trcdm2dc = ln_dm2dc .OR. ( ln_cpl .AND. ncpl_qsr_freq /= 1 ) 104 104 l_trcdm2dc = l_trcdm2dc .AND. .NOT. lk_offline 105 IF( l_trcdm2dc .AND. lwp ) & 106 & CALL ctl_warn(' Coupling with passive tracers and used of diurnal cycle. & 107 & Computation of a daily mean shortwave for some biogeochemical models) ') 108 109 ! 110 IF( nn_cla == 1 ) & 111 & CALL ctl_stop( ' Cross Land Advection not yet implemented with passive tracer ; nn_cla must be 0' ) 105 IF( l_trcdm2dc .AND. lwp ) CALL ctl_warn( 'Coupling with passive tracers and used of diurnal cycle.', & 106 & 'Computation of a daily mean shortwave for some biogeochemical models ' ) 112 107 ! 113 108 END SUBROUTINE trc_ini_ctl 109 114 110 115 111 SUBROUTINE trc_ini_inv … … 157 153 END SUBROUTINE trc_ini_inv 158 154 155 159 156 SUBROUTINE trc_ini_sms 160 157 !!---------------------------------------------------------------------- … … 196 193 ! 197 194 END SUBROUTINE trc_ini_trp 195 198 196 199 197 SUBROUTINE trc_ini_state … … 248 246 END SUBROUTINE trc_ini_state 249 247 248 250 249 SUBROUTINE top_alloc 251 250 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.