Changeset 14834 for NEMO/trunk/src/OCE/TRA
- Timestamp:
- 2021-05-11T11:24:44+02:00 (3 years ago)
- Location:
- NEMO/trunk/src/OCE/TRA
- Files:
-
- 2 deleted
- 20 edited
- 3 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/TRA/eosbn2.F90
r14131 r14834 577 577 578 578 SUBROUTINE eos_insitu_pot_2d( pts, prhop ) 579 !! 580 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 581 ! ! 2 : salinity [psu] 582 REAL(wp), DIMENSION(:,:) , INTENT( out) :: prhop ! potential density (surface referenced) 583 !! 584 CALL eos_insitu_pot_2d_t( pts, is_tile(pts), prhop, is_tile(prhop) ) 585 END SUBROUTINE eos_insitu_pot_2d 586 587 588 SUBROUTINE eos_insitu_pot_2d_t( pts, ktts, prhop, ktrhop ) 579 589 !!---------------------------------------------------------------------- 580 590 !! *** ROUTINE eos_insitu_pot *** … … 589 599 !! 590 600 !!---------------------------------------------------------------------- 591 REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 601 INTEGER , INTENT(in ) :: ktts, ktrhop 602 REAL(wp), DIMENSION(A2D_T(ktts),JPTS), INTENT(in ) :: pts ! 1 : potential temperature [Celsius] 592 603 ! ! 2 : salinity [psu] 593 REAL(wp), DIMENSION( jpi,jpj), INTENT( out) :: prhop ! potential density (surface referenced)604 REAL(wp), DIMENSION(A2D_T(ktrhop) ), INTENT( out) :: prhop ! potential density (surface referenced) 594 605 ! 595 606 INTEGER :: ji, jj, jk, jsmp ! dummy loop indices … … 606 617 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 607 618 ! 608 DO_2D( 1, 1, 1, 1)609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 619 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 620 ! 621 zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature 622 zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity 623 ztm = tmask(ji,jj,1) ! tmask 624 ! 625 zn0 = (((((EOS060*zt & 626 & + EOS150*zs+EOS050)*zt & 627 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 628 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 629 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 630 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 631 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 632 ! 633 ! 634 prhop(ji,jj) = zn0 * ztm ! potential density referenced at the surface 635 ! 636 END_2D 626 637 627 638 CASE( np_seos ) !== simplified EOS ==! 628 639 ! 629 DO_2D( 1, 1, 1, 1)640 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 630 641 zt = pts (ji,jj,jp_tem) - 10._wp 631 642 zs = pts (ji,jj,jp_sal) - 35._wp … … 646 657 IF( ln_timing ) CALL timing_stop('eos-pot') 647 658 ! 648 END SUBROUTINE eos_insitu_pot_2d 659 END SUBROUTINE eos_insitu_pot_2d_t 649 660 650 661 -
NEMO/trunk/src/OCE/TRA/traadv.F90
r14433 r14834 18 18 USE oce ! ocean dynamics and active tracers 19 19 USE dom_oce ! ocean space and time domain 20 ! TEMP: [tiling] This change not necessary after extended haloes development20 ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 21 21 USE domtile 22 22 USE domvvl ! variable vertical scale factors … … 25 25 USE traadv_cen ! centered scheme (tra_adv_cen routine) 26 26 USE traadv_fct ! FCT scheme (tra_adv_fct routine) 27 USE traadv_fct_lf ! FCT scheme (tra_adv_fct routine - loop fusion version)28 27 USE traadv_mus ! MUSCL scheme (tra_adv_mus routine) 29 USE traadv_mus_lf ! MUSCL scheme (tra_adv_mus routine - loop fusion version)30 28 USE traadv_ubs ! UBS scheme (tra_adv_ubs routine) 31 29 USE traadv_qck ! QUICKEST scheme (tra_adv_qck routine) … … 61 59 LOGICAL :: ln_traadv_qck ! QUICKEST scheme flag 62 60 63 INTEGER :: nadv ! choice of the type of advection scheme61 INTEGER, PUBLIC :: nadv ! choice of the type of advection scheme 64 62 ! ! associated indices: 65 INTEGER, PARAMETER :: np_NO_adv = 0 ! no T-S advection66 INTEGER, PARAMETER :: np_CEN = 1 ! 2nd/4th order centered scheme67 INTEGER, PARAMETER :: np_FCT = 2 ! 2nd/4th order Flux Corrected Transport scheme68 INTEGER, PARAMETER :: np_MUS = 3 ! MUSCL scheme69 INTEGER, PARAMETER :: np_UBS = 4 ! 3rd order Upstream Biased Scheme70 INTEGER, PARAMETER :: np_QCK = 5 ! QUICK scheme63 INTEGER, PARAMETER, PUBLIC :: np_NO_adv = 0 ! no T-S advection 64 INTEGER, PARAMETER, PUBLIC :: np_CEN = 1 ! 2nd/4th order centered scheme 65 INTEGER, PARAMETER, PUBLIC :: np_FCT = 2 ! 2nd/4th order Flux Corrected Transport scheme 66 INTEGER, PARAMETER, PUBLIC :: np_MUS = 3 ! MUSCL scheme 67 INTEGER, PARAMETER, PUBLIC :: np_UBS = 4 ! 3rd order Upstream Biased Scheme 68 INTEGER, PARAMETER, PUBLIC :: np_QCK = 5 ! QUICK scheme 71 69 72 70 !! * Substitutions … … 93 91 ! 94 92 INTEGER :: ji, jj, jk ! dummy loop index 95 ! TEMP: [tiling] This change not necessary and can be A2D(nn_hls) if using XIOS (subdomain support)93 ! TEMP: [tiling] This change not necessary and can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 96 94 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zuu, zvv, zww ! 3D workspace 97 95 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds 98 ! TEMP: [tiling] This change not necessary after extra haloes development96 ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 99 97 LOGICAL :: lskip 100 98 !!---------------------------------------------------------------------- … … 104 102 lskip = .FALSE. 105 103 106 ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support)107 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile104 ! TEMP: [tiling] These changes not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 105 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 108 106 ALLOCATE( zuu(jpi,jpj,jpk), zvv(jpi,jpj,jpk), zww(jpi,jpj,jpk) ) 109 107 ENDIF 110 108 111 ! TEMP: [tiling] These changes not necessary after extra haloes development (lbc_lnk removed from tra_adv_*) and if XIOS has subdomain support (ldf_eiv_dia) 112 IF( nadv /= np_CEN .OR. (nadv == np_CEN .AND. nn_cen_h == 4) .OR. ln_ldfeiv_dia ) THEN 113 IF( ln_tile ) THEN 114 IF( ntile == 1 ) THEN 115 CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) 116 ELSE 117 lskip = .TRUE. 118 ENDIF 109 ! TEMP: [tiling] These changes not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 110 IF( ln_tile .AND. nadv == np_FCT ) THEN 111 IF( ntile == 1 ) THEN 112 CALL dom_tile_stop( ldhold=.TRUE. ) 113 ELSE 114 lskip = .TRUE. 119 115 ENDIF 120 116 ENDIF … … 122 118 ! !== effective transport ==! 123 119 IF( ln_wave .AND. ln_sdw ) THEN 124 DO_3D ( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )120 DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 125 121 zuu(ji,jj,jk) = e2u (ji,jj) * e3u(ji,jj,jk,Kmm) * ( uu(ji,jj,jk,Kmm) + usd(ji,jj,jk) ) 126 122 zvv(ji,jj,jk) = e1v (ji,jj) * e3v(ji,jj,jk,Kmm) * ( vv(ji,jj,jk,Kmm) + vsd(ji,jj,jk) ) … … 128 124 END_3D 129 125 ELSE 130 DO_3D ( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )126 DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 131 127 zuu(ji,jj,jk) = e2u (ji,jj) * e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) ! eulerian transport only 132 128 zvv(ji,jj,jk) = e1v (ji,jj) * e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) … … 136 132 ! 137 133 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! add z-tilde and/or vvl corrections 138 DO_3D ( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )134 DO_3D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 139 135 zuu(ji,jj,jk) = zuu(ji,jj,jk) + un_td(ji,jj,jk) 140 136 zvv(ji,jj,jk) = zvv(ji,jj,jk) + vn_td(ji,jj,jk) … … 142 138 ENDIF 143 139 ! 144 DO_2D ( nn_hls, nn_hls, nn_hls, nn_hls)140 DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 145 141 zuu(ji,jj,jpk) = 0._wp ! no transport trough the bottom 146 142 zvv(ji,jj,jpk) = 0._wp … … 148 144 END_2D 149 145 ! 150 ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support)151 146 IF( ln_ldfeiv .AND. .NOT. ln_traldf_triad ) & 152 & CALL ldf_eiv_trp( kt, nit000, zuu(A2D(nn_hls),:), zvv(A2D(nn_hls),:), zww(A2D(nn_hls),:), & 153 & 'TRA', Kmm, Krhs ) ! add the eiv transport (if necessary) 154 ! 155 IF( ln_mle ) CALL tra_mle_trp( kt, nit000, zuu(A2D(nn_hls),:), zvv(A2D(nn_hls),:), zww(A2D(nn_hls),:), & 156 & 'TRA', Kmm ) ! add the mle transport (if necessary) 157 ! 158 ! TEMP: [tiling] This change not necessary if using XIOS (subdomain support) 159 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 147 & CALL ldf_eiv_trp( kt, nit000, zuu, zvv, zww, 'TRA', Kmm, Krhs ) ! add the eiv transport (if necessary) 148 ! 149 IF( ln_mle ) CALL tra_mle_trp( kt, nit000, zuu, zvv, zww, 'TRA', Kmm ) ! add the mle transport (if necessary) 150 ! 151 ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 152 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 160 153 CALL iom_put( "uocetr_eff", zuu ) ! output effective transport 161 154 CALL iom_put( "vocetr_eff", zvv ) … … 163 156 ENDIF 164 157 ! 165 166 ! TEMP: [tiling] This c hange not necessary if using XIOS (subdomain support)158 !!gm ??? 159 ! TEMP: [tiling] This copy-in not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 167 160 CALL dia_ptr( kt, Kmm, zvv(A2D(nn_hls),:) ) ! diagnose the effective MSF 168 161 !!gm ??? 169 162 ! 170 163 … … 178 171 ! 179 172 CASE ( np_CEN ) ! Centered scheme : 2nd / 4th order 180 IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kmm), 'T', 1. )181 173 CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v ) 182 174 CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order 183 IF (nn_hls.EQ.2) THEN184 CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1., pts(:,:,:,:,Kmm), 'T', 1.)185 CALL lbc_lnk( 'traadv', zuu(:,:,:), 'U', -1., zvv(:,:,:), 'V', -1., zww(:,:,:), 'W', 1.)186 #if defined key_loop_fusion187 CALL tra_adv_fct_lf ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )188 #else189 175 CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v ) 190 #endif191 ELSE192 CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )193 END IF194 176 CASE ( np_MUS ) ! MUSCL 195 IF (nn_hls.EQ.2) THEN196 CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)197 #if defined key_loop_fusion198 CALL tra_adv_mus_lf ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )199 #else200 177 CALL tra_adv_mus ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups ) 201 #endif202 ELSE203 CALL tra_adv_mus ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )204 END IF205 178 CASE ( np_UBS ) ! UBS 206 IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)207 179 CALL tra_adv_ubs ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v ) 208 180 CASE ( np_QCK ) ! QUICKEST 209 IF (nn_hls.EQ.2) THEN210 CALL lbc_lnk( 'traadv', zuu(:,:,:), 'U', -1., zvv(:,:,:), 'V', -1.)211 CALL lbc_lnk( 'traadv', pts(:,:,:,:,Kbb), 'T', 1.)212 END IF213 181 CALL tra_adv_qck ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs ) 214 182 ! … … 225 193 ENDIF 226 194 227 ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_adv_*) and if XIOS has subdomain support (ldf_eiv_dia) 228 IF( ln_tile .AND. ntile == 0 ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 ) 229 195 ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 196 IF( ln_tile .AND. .NOT. l_istiled ) CALL dom_tile_start( ldhold=.TRUE. ) 230 197 ENDIF 231 198 ! ! print mean trends (used for debugging) … … 233 200 & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 234 201 235 ! TEMP: [tiling] This change not necessary if using XIOS (subdomain support)236 IF( ntile == 0.OR. ntile == nijtile ) THEN ! Do only for the full domain202 ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 203 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only for the full domain 237 204 DEALLOCATE( zuu, zvv, zww ) 238 205 ENDIF … … 306 273 CALL ctl_stop( 'tra_adv_init: FCT scheme, choose 2nd or 4th order' ) 307 274 ENDIF 275 ! TEMP: [tiling] This change not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 276 IF( ln_traadv_fct .AND. ln_tile ) THEN 277 CALL ctl_warn( 'tra_adv_init: FCT scheme does not yet work with tiling' ) 278 ENDIF 308 279 IF( ln_traadv_ubs .AND. ( nn_ubs_v /= 2 .AND. nn_ubs_v /= 4 ) ) THEN ! UBS 309 280 CALL ctl_stop( 'tra_adv_init: UBS scheme, choose 2nd or 4th order' ) -
NEMO/trunk/src/OCE/TRA/traadv_cen.F90
r14433 r14834 23 23 USE trc_oce ! share passive tracers/Ocean variables 24 24 USE lib_mpp ! MPP library 25 #if defined key_loop_fusion 26 USE traadv_cen_lf ! centered scheme (tra_adv_cen routine - loop fusion version) 27 #endif 25 28 26 29 IMPLICIT NONE … … 71 74 INTEGER , INTENT(in ) :: kn_cen_h ! =2/4 (2nd or 4th order scheme) 72 75 INTEGER , INTENT(in ) :: kn_cen_v ! =2/4 (2nd or 4th order scheme) 73 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)76 ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 74 77 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components 75 78 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation … … 82 85 !!---------------------------------------------------------------------- 83 86 ! 84 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 87 #if defined key_loop_fusion 88 CALL tra_adv_cen_lf ( kt, nit000, cdtype, pU, pV, pW, Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v ) 89 #else 90 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 85 91 IF( kt == kit000 ) THEN 86 92 IF(lwp) WRITE(numout,*) … … 119 125 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 120 126 END_3D 121 IF (nn_hls .EQ.1) CALL lbc_lnk( 'traadv_cen', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp) ! Lateral boundary cond.127 IF (nn_hls==1) CALL lbc_lnk( 'traadv_cen', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp, ld4only= .TRUE. ) ! Lateral boundary cond. 122 128 ! 123 129 DO_3D( nn_hls-1, 0, nn_hls-1, 0, 1, jpkm1 ) ! Horizontal advective fluxes … … 131 137 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v 132 138 END_3D 133 IF (nn_hls .EQ.1) CALL lbc_lnk( 'traadv_cen', zwx, 'U', -1. , zwy, 'V', -1. )139 IF (nn_hls==1) CALL lbc_lnk( 'traadv_cen', zwx, 'U', -1. , zwy, 'V', -1. ) 134 140 ! 135 141 CASE DEFAULT … … 184 190 END DO 185 191 ! 192 #endif 186 193 END SUBROUTINE tra_adv_cen 187 194 -
NEMO/trunk/src/OCE/TRA/traadv_fct.F90
r14820 r14834 34 34 PUBLIC tra_adv_fct ! called by traadv.F90 35 35 PUBLIC interp_4th_cpt ! called by traadv_cen.F90 36 PUBLIC tridia_solver ! called by traadv_fct_lf.F9037 PUBLIC nonosc ! called by traadv_fct_lf.F90 - key_agrif38 36 39 37 LOGICAL :: l_trd ! flag to compute trends … … 81 79 INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4) 82 80 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 83 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)81 ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case 84 82 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components 85 83 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation … … 95 93 !!---------------------------------------------------------------------- 96 94 ! 97 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 95 #if defined key_loop_fusion 96 CALL tra_adv_fct_lf ( kt, nit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) 97 #else 98 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 98 99 IF( kt == kit000 ) THEN 99 100 IF(lwp) WRITE(numout,*) … … 136 137 ! If adaptive vertical advection, check if it is needed on this PE at this time 137 138 IF( ln_zad_Aimp ) THEN 138 IF( MAXVAL( ABS( wi(A2D( nn_hls),:) ) ) > 0._wp ) ll_zAimp = .TRUE.139 IF( MAXVAL( ABS( wi(A2D(1),:) ) ) > 0._wp ) ll_zAimp = .TRUE. 139 140 END IF 140 141 ! If active adaptive vertical advection, build tridiagonal matrix … … 239 240 END DO 240 241 ! NOTE [ comm_cleanup ] : need to change sign to ensure halo 1 - halo 2 compatibility 241 CALL lbc_lnk( 'traadv_fct', zltu, 'T', -1.0_wp , zltv, 'T', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn)242 ! 243 DO_3D( 1, 0, 1, 0, 1, jpkm1 )242 CALL lbc_lnk( 'traadv_fct', zltu, 'T', -1.0_wp , zltv, 'T', -1.0_wp, ld4only= .TRUE. ) ! Lateral boundary cond. (unchanged sgn) 243 ! 244 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) 244 245 zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! 2 x C2 interpolation of T at u- & v-points 245 246 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) … … 254 255 & ) - zwy(ji,jj,jk) 255 256 END_3D 256 IF (nn_hls.EQ.2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp, zwy, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn)257 257 ! 258 258 CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested 259 259 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 260 260 ztv(:,:,jpk) = 0._wp 261 DO_3D( nn_hls , nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) ! 1st derivative (gradient)261 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) ! 1st derivative (gradient) 262 262 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 263 263 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 264 264 END_3D 265 IF (nn_hls .EQ.1) CALL lbc_lnk( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp) ! Lateral boundary cond. (unchanged sgn)265 IF (nn_hls==1) CALL lbc_lnk( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp, ld4only= .TRUE. ) ! Lateral boundary cond. (unchanged sgn) 266 266 ! 267 267 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Horizontal advective fluxes … … 275 275 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 276 276 END_3D 277 IF (nn_hls .EQ.2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn)277 IF (nn_hls==2) CALL lbc_lnk( 'traadv_fct', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 278 278 ! 279 279 END SELECT … … 298 298 ENDIF 299 299 ! 300 IF (nn_hls .EQ.1) THEN300 IF (nn_hls==1) THEN 301 301 CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'T', 1.0_wp ) 302 302 ELSE … … 381 381 ENDIF 382 382 ! 383 #endif 383 384 END SUBROUTINE tra_adv_fct 384 385 … … 456 457 END_2D 457 458 END DO 458 IF (nn_hls .EQ.1) CALL lbc_lnk( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp) ! lateral boundary cond. (unchanged sign)459 IF (nn_hls==1) CALL lbc_lnk( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp, ld4only= .TRUE. ) ! lateral boundary cond. (unchanged sign) 459 460 460 461 ! 3. monotonic flux in the i & j direction (paa & pbb) … … 677 678 END SUBROUTINE tridia_solver 678 679 680 #if defined key_loop_fusion 681 #define tracer_flux_i(out,zfp,zfm,ji,jj,jk) \ 682 zfp = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) ; \ 683 zfm = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) ) ; \ 684 out = 0.5 * ( zfp * pt(ji,jj,jk,jn,Kbb) + zfm * pt(ji+1,jj,jk,jn,Kbb) ) 685 686 #define tracer_flux_j(out,zfp,zfm,ji,jj,jk) \ 687 zfp = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) ) ; \ 688 zfm = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) ) ; \ 689 out = 0.5 * ( zfp * pt(ji,jj,jk,jn,Kbb) + zfm * pt(ji,jj+1,jk,jn,Kbb) ) 690 691 SUBROUTINE tra_adv_fct_lf( kt, kit000, cdtype, p2dt, pU, pV, pW, & 692 & Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) 693 !!---------------------------------------------------------------------- 694 !! *** ROUTINE tra_adv_fct *** 695 !! 696 !! ** Purpose : Compute the now trend due to total advection of tracers 697 !! and add it to the general trend of tracer equations 698 !! 699 !! ** Method : - 2nd or 4th FCT scheme on the horizontal direction 700 !! (choice through the value of kn_fct) 701 !! - on the vertical the 4th order is a compact scheme 702 !! - corrected flux (monotonic correction) 703 !! 704 !! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends 705 !! - send trends to trdtra module for further diagnostics (l_trdtra=T) 706 !! - poleward advective heat and salt transport (ln_diaptr=T) 707 !!---------------------------------------------------------------------- 708 INTEGER , INTENT(in ) :: kt ! ocean time-step index 709 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices 710 INTEGER , INTENT(in ) :: kit000 ! first time step index 711 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 712 INTEGER , INTENT(in ) :: kjpt ! number of tracers 713 INTEGER , INTENT(in ) :: kn_fct_h ! order of the FCT scheme (=2 or 4) 714 INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4) 715 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 716 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components 717 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 718 ! 719 INTEGER :: ji, jj, jk, jn ! dummy loop indices 720 REAL(wp) :: ztra ! local scalar 721 REAL(wp) :: zwx_im1, zfp_ui, zfp_ui_m1, zfp_vj, zfp_vj_m1, zfp_wk, zC2t_u, zC4t_u ! - - 722 REAL(wp) :: zwy_jm1, zfm_ui, zfm_ui_m1, zfm_vj, zfm_vj_m1, zfm_wk, zC2t_v, zC4t_v ! - - 723 REAL(wp) :: ztu, ztv, ztu_im1, ztu_ip1, ztv_jm1, ztv_jp1 724 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwx_3d, zwy_3d, zwz, ztw, zltu_3d, zltv_3d 725 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry 726 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup 727 LOGICAL :: ll_zAimp ! flag to apply adaptive implicit vertical advection 728 !!---------------------------------------------------------------------- 729 ! 730 IF( kt == kit000 ) THEN 731 IF(lwp) WRITE(numout,*) 732 IF(lwp) WRITE(numout,*) 'tra_adv_fct_lf : FCT advection scheme on ', cdtype 733 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 734 ENDIF 735 !! -- init to 0 736 zwx_3d(:,:,:) = 0._wp 737 zwy_3d(:,:,:) = 0._wp 738 zwz(:,:,:) = 0._wp 739 zwi(:,:,:) = 0._wp 740 ! 741 l_trd = .FALSE. ! set local switches 742 l_hst = .FALSE. 743 l_ptr = .FALSE. 744 ll_zAimp = .FALSE. 745 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 746 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 747 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 748 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 749 ! 750 IF( l_trd .OR. l_hst ) THEN 751 ALLOCATE( ztrdx(jpi,jpj,jpk), ztrdy(jpi,jpj,jpk), ztrdz(jpi,jpj,jpk) ) 752 ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp 753 ENDIF 754 ! 755 IF( l_ptr ) THEN 756 ALLOCATE( zptry(jpi,jpj,jpk) ) 757 zptry(:,:,:) = 0._wp 758 ENDIF 759 ! 760 ! If adaptive vertical advection, check if it is needed on this PE at this time 761 IF( ln_zad_Aimp ) THEN 762 IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE. 763 END IF 764 ! If active adaptive vertical advection, build tridiagonal matrix 765 IF( ll_zAimp ) THEN 766 ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 767 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 768 zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) & 769 & / e3t(ji,jj,jk,Krhs) 770 zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t(ji,jj,jk,Krhs) 771 zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs) 772 END_3D 773 END IF 774 ! 775 DO jn = 1, kjpt !== loop over the tracers ==! 776 ! 777 ! !== upstream advection with initial mass fluxes & intermediate update ==! 778 ! !* upstream tracer flux in the k direction *! 779 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 780 zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) ) 781 zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) ) 782 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk) 783 END_3D 784 IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked) 785 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 786 DO_2D( 1, 1, 1, 1 ) 787 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 788 END_2D 789 ELSE ! no cavities: only at the ocean surface 790 DO_2D( 1, 1, 1, 1 ) 791 zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) 792 END_2D 793 ENDIF 794 ENDIF 795 ! 796 ! !* upstream tracer flux in the i and j direction 797 DO jk = 1, jpkm1 798 DO jj = 1, jpj-1 799 tracer_flux_i(zwx_3d(1,jj,jk),zfp_ui,zfm_ui,1,jj,jk) 800 tracer_flux_j(zwy_3d(1,jj,jk),zfp_vj,zfm_vj,1,jj,jk) 801 END DO 802 DO ji = 1, jpi-1 803 tracer_flux_i(zwx_3d(ji,1,jk),zfp_ui,zfm_ui,ji,1,jk) 804 tracer_flux_j(zwy_3d(ji,1,jk),zfp_vj,zfm_vj,ji,1,jk) 805 END DO 806 DO_2D( 1, 1, 1, 1 ) 807 tracer_flux_i(zwx_3d(ji,jj,jk),zfp_ui,zfm_ui,ji,jj,jk) 808 tracer_flux_i(zwx_im1,zfp_ui_m1,zfm_ui_m1,ji-1,jj,jk) 809 tracer_flux_j(zwy_3d(ji,jj,jk),zfp_vj,zfm_vj,ji,jj,jk) 810 tracer_flux_j(zwy_jm1,zfp_vj_m1,zfm_vj_m1,ji,jj-1,jk) 811 ztra = - ( zwx_3d(ji,jj,jk) - zwx_im1 + zwy_3d(ji,jj,jk) - zwy_jm1 + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) 812 ! ! update and guess with monotonic sheme 813 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra & 814 & / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) 815 zwi(ji,jj,jk) = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & 816 & / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 817 END_2D 818 END DO 819 820 IF ( ll_zAimp ) THEN 821 CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 ) 822 ! 823 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 824 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 825 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 826 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 827 ztw(ji,jj,jk) = 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 828 zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes 829 END_3D 830 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 831 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 832 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 833 END_3D 834 ! 835 END IF 836 ! 837 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) 838 ztrdx(:,:,:) = zwx_3d(:,:,:) ; ztrdy(:,:,:) = zwy_3d(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 839 END IF 840 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 841 IF( l_ptr ) zptry(:,:,:) = zwy_3d(:,:,:) 842 ! 843 ! !== anti-diffusive flux : high order minus low order ==! 844 ! 845 SELECT CASE( kn_fct_h ) !* horizontal anti-diffusive fluxes 846 ! 847 CASE( 2 ) !- 2nd order centered 848 DO_3D( 2, 1, 2, 1, 1, jpkm1 ) 849 zwx_3d(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx_3d(ji,jj,jk) 850 zwy_3d(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy_3d(ji,jj,jk) 851 END_3D 852 ! 853 CASE( 4 ) !- 4th order centered 854 zltu_3d(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 855 zltv_3d(:,:,jpk) = 0._wp 856 ! ! Laplacian 857 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! 2nd derivative * 1/ 6 858 ! ! 1st derivative (gradient) 859 ztu = ( pt(ji+1,jj,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 860 ztu_im1 = ( pt(ji,jj,jk,jn,Kmm) - pt(ji-1,jj,jk,jn,Kmm) ) * umask(ji-1,jj,jk) 861 ztv = ( pt(ji,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 862 ztv_jm1 = ( pt(ji,jj,jk,jn,Kmm) - pt(ji,jj-1,jk,jn,Kmm) ) * vmask(ji,jj-1,jk) 863 ! ! 2nd derivative * 1/ 6 864 zltu_3d(ji,jj,jk) = ( ztu + ztu_im1 ) * r1_6 865 zltv_3d(ji,jj,jk) = ( ztv + ztv_jm1 ) * r1_6 866 END_2D 867 END DO 868 ! NOTE [ comm_cleanup ] : need to change sign to ensure halo 1 - halo 2 compatibility 869 CALL lbc_lnk( 'traadv_fct', zltu_3d, 'T', -1.0_wp , zltv_3d, 'T', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 870 ! 871 DO_3D( 2, 1, 2, 1, 1, jpkm1 ) 872 zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! 2 x C2 interpolation of T at u- & v-points 873 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 874 ! ! C4 minus upstream advective fluxes 875 ! round brackets added to fix the order of floating point operations 876 ! needed to ensure halo 1 - halo 2 compatibility 877 zwx_3d(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( zC2t_u + ( zltu_3d(ji,jj,jk) - zltu_3d(ji+1,jj,jk) & 878 & ) & ! bracket for halo 1 - halo 2 compatibility 879 & ) - zwx_3d(ji,jj,jk) 880 zwy_3d(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( zC2t_v + ( zltv_3d(ji,jj,jk) - zltv_3d(ji,jj+1,jk) & 881 & ) & ! bracket for halo 1 - halo 2 compatibility 882 & ) - zwy_3d(ji,jj,jk) 883 END_3D 884 ! 885 CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested 886 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Horizontal advective fluxes 887 ztu_im1 = ( pt(ji ,jj ,jk,jn,Kmm) - pt(ji-1,jj,jk,jn,Kmm) ) * umask(ji-1,jj,jk) 888 ztu_ip1 = ( pt(ji+2,jj ,jk,jn,Kmm) - pt(ji+1,jj,jk,jn,Kmm) ) * umask(ji+1,jj,jk) 889 890 ztv_jm1 = ( pt(ji,jj ,jk,jn,Kmm) - pt(ji,jj-1,jk,jn,Kmm) ) * vmask(ji,jj-1,jk) 891 ztv_jp1 = ( pt(ji,jj+2,jk,jn,Kmm) - pt(ji,jj+1,jk,jn,Kmm) ) * vmask(ji,jj+1,jk) 892 zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! 2 x C2 interpolation of T at u- & v-points (x2) 893 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 894 ! ! C4 interpolation of T at u- & v-points (x2) 895 zC4t_u = zC2t_u + r1_6 * ( ztu_im1 - ztu_ip1 ) 896 zC4t_v = zC2t_v + r1_6 * ( ztv_jm1 - ztv_jp1 ) 897 ! ! C4 minus upstream advective fluxes 898 zwx_3d(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx_3d(ji,jj,jk) 899 zwy_3d(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy_3d(ji,jj,jk) 900 END_3D 901 CALL lbc_lnk( 'traadv_fct', zwx_3d, 'U', -1.0_wp , zwy_3d, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 902 ! 903 END SELECT 904 ! 905 SELECT CASE( kn_fct_v ) !* vertical anti-diffusive fluxes (w-masked interior values) 906 ! 907 CASE( 2 ) !- 2nd order centered 908 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 909 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) & 910 & - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 911 END_3D 912 ! 913 CASE( 4 ) !- 4th order COMPACT 914 CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! zwt = COMPACT interpolation of T at w-point 915 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 916 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 917 END_3D 918 ! 919 END SELECT 920 IF( ln_linssh ) THEN ! top ocean value: high order = upstream ==>> zwz=0 921 zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked 922 ENDIF 923 ! 924 CALL lbc_lnk( 'traadv_fct', zwi, 'T', 1.0_wp) 925 ! 926 IF ( ll_zAimp ) THEN 927 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) !* trend and after field with monotonic scheme 928 ! ! total intermediate advective trends 929 ztra = - ( zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj ,jk ) & 930 & + zwy_3d(ji,jj,jk) - zwy_3d(ji ,jj-1,jk ) & 931 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 932 ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 933 END_3D 934 ! 935 CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 936 ! 937 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 938 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 939 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 940 zwz(ji,jj,jk) = zwz(ji,jj,jk) + 0.5 * e1e2t(ji,jj) * ( zfp_wk * ztw(ji,jj,jk) + zfm_wk * ztw(ji,jj,jk-1) ) * wmask(ji,jj,jk) 941 END_3D 942 END IF 943 ! 944 ! !== monotonicity algorithm ==! 945 ! 946 CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx_3d, zwy_3d, zwz, zwi, p2dt ) 947 ! 948 ! !== final trend with corrected fluxes ==! 949 ! 950 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 951 ztra = - ( zwx_3d(ji,jj,jk) - zwx_3d(ji-1,jj ,jk ) & 952 & + zwy_3d(ji,jj,jk) - zwy_3d(ji ,jj-1,jk ) & 953 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 954 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm) 955 zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 956 END_3D 957 ! 958 IF ( ll_zAimp ) THEN 959 ! 960 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 961 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! Interior value ( multiplied by wmask) 962 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 963 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 964 ztw(ji,jj,jk) = - 0.5 * e1e2t(ji,jj) * ( zfp_wk * zwi(ji,jj,jk) + zfm_wk * zwi(ji,jj,jk-1) ) * wmask(ji,jj,jk) 965 zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic 966 END_3D 967 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 968 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 969 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 970 END_3D 971 END IF 972 ! NOT TESTED - NEED l_trd OR l_hst TRUE 973 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport 974 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx_3d(:,:,:) ! <<< add anti-diffusive fluxes 975 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy_3d(:,:,:) ! to upstream fluxes 976 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! 977 ! 978 IF( l_trd ) THEN ! trend diagnostics 979 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) ) 980 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) ) 981 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, pt(:,:,:,jn,Kmm) ) 982 ENDIF 983 ! ! heat/salt transport 984 IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) ) 985 ! 986 ENDIF 987 ! NOT TESTED - NEED l_ptr TRUE 988 IF( l_ptr ) THEN ! "Poleward" transports 989 zptry(:,:,:) = zptry(:,:,:) + zwy_3d(:,:,:) ! <<< add anti-diffusive fluxes 990 CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) ) 991 ENDIF 992 ! 993 END DO ! end of tracer loop 994 ! 995 IF ( ll_zAimp ) THEN 996 DEALLOCATE( zwdia, zwinf, zwsup ) 997 ENDIF 998 IF( l_trd .OR. l_hst ) THEN 999 DEALLOCATE( ztrdx, ztrdy, ztrdz ) 1000 ENDIF 1001 IF( l_ptr ) THEN 1002 DEALLOCATE( zptry ) 1003 ENDIF 1004 ! 1005 END SUBROUTINE tra_adv_fct_lf 1006 #endif 679 1007 !!====================================================================== 680 1008 END MODULE traadv_fct -
NEMO/trunk/src/OCE/TRA/traadv_mus.F90
r14433 r14834 81 81 LOGICAL , INTENT(in ) :: ld_msc_ups ! use upstream scheme within muscl 82 82 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 83 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)83 ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 84 84 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components 85 85 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation … … 93 93 !!---------------------------------------------------------------------- 94 94 ! 95 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile95 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 96 96 IF( kt == kit000 ) THEN 97 97 IF(lwp) WRITE(numout,*) … … 139 139 zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( pt(ji,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) 140 140 END_3D 141 ! lateral boundary conditions (changed sign)142 IF ( nn_hls.EQ.1 ) CALL lbc_lnk( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp )143 141 ! !-- Slopes of tracer 144 142 zslpx(:,:,jpk) = 0._wp ! bottom values 145 143 zslpy(:,:,jpk) = 0._wp 146 DO_3D( nn_hls-1, 1, nn_hls-1,1, 1, jpkm1 )144 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 147 145 zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji-1,jj ,jk) ) & 148 146 & * ( 0.25 + SIGN( 0.25_wp, zwx(ji,jj,jk) * zwx(ji-1,jj ,jk) ) ) … … 151 149 END_3D 152 150 ! 153 DO_3D( nn_hls-1, 1, nn_hls-1,1, 1, jpkm1 ) !-- Slopes limitation151 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) !-- Slopes limitation 154 152 zslpx(ji,jj,jk) = SIGN( 1.0_wp, zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji ,jj,jk) ), & 155 153 & 2.*ABS( zwx (ji-1,jj,jk) ), & … … 159 157 & 2.*ABS( zwy (ji,jj ,jk) ) ) 160 158 END_3D 161 ! 162 DO_3D( nn_hls-1, 0, nn_hls-1, 0, 1, jpkm1 ) !-- MUSCL horizontal advective fluxes 159 ! NOTE [ comm_cleanup ] : need to change sign to ensure halo 1 - halo 2 compatibility 160 IF ( nn_hls==1 ) CALL lbc_lnk( 'traadv_mus', zslpx, 'T', -1.0_wp , zslpy, 'T', -1.0_wp ) ! lateral boundary conditions (changed sign) 161 ! 162 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) !-- MUSCL horizontal advective fluxes 163 163 ! MUSCL fluxes 164 164 z0u = SIGN( 0.5_wp, pU(ji,jj,jk) ) … … 176 176 zwy(ji,jj,jk) = pV(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) 177 177 END_3D 178 IF ( nn_hls.EQ.1 ) CALL lbc_lnk( 'traadv_mus', zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp ) ! lateral boundary conditions (changed sign)179 178 ! 180 179 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Tracer advective trend -
NEMO/trunk/src/OCE/TRA/traadv_qck.F90
r14433 r14834 27 27 USE lbclnk ! ocean lateral boundary condition (or mpp link) 28 28 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 29 #if defined key_loop_fusion 30 USE traadv_qck_lf ! QCK scheme (tra_adv_qck routine - loop fusion version) 31 #endif 29 32 30 33 IMPLICIT NONE … … 91 94 INTEGER , INTENT(in ) :: kjpt ! number of tracers 92 95 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 93 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)96 ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 94 97 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume transport components 95 98 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 96 99 !!---------------------------------------------------------------------- 97 100 ! 98 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 101 #if defined key_loop_fusion 102 CALL tra_adv_qck_lf ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs ) 103 #else 104 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 99 105 IF( kt == kit000 ) THEN 100 106 IF(lwp) WRITE(numout,*) … … 117 123 CALL tra_adv_cen2_k( kt, cdtype, pW, Kmm, pt, kjpt, Krhs ) 118 124 ! 125 #endif 119 126 END SUBROUTINE tra_adv_qck 120 127 … … 129 136 INTEGER , INTENT(in ) :: kjpt ! number of tracers 130 137 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 131 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)138 ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 132 139 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU ! i-velocity components 133 140 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! active tracers and RHS of tracer equation … … 149 156 zfd(ji,jj,jk) = pt(ji+1,jj,jk,jn,Kbb) ! Downstream in the x-direction for the tracer 150 157 END_3D 151 IF (nn_hls .EQ.1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp) ! Lateral boundary conditions158 IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. ) ! Lateral boundary conditions 152 159 153 160 ! … … 167 174 END_3D 168 175 !--- Lateral boundary conditions 169 IF (nn_hls .EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwx(:,:,:), 'T', 1.0_wp )176 IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwx(:,:,:), 'T', 1.0_wp ) 170 177 171 178 !--- QUICKEST scheme … … 176 183 zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2. 177 184 END_3D 178 IF (nn_hls .EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp) ! Lateral boundary conditions185 IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. ) ! Lateral boundary conditions 179 186 180 187 ! … … 214 221 INTEGER , INTENT(in ) :: kjpt ! number of tracers 215 222 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 216 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)223 ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 217 224 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pV ! j-velocity components 218 225 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! active tracers and RHS of tracer equation … … 229 236 zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0 230 237 ! 231 DO jk = 1, jpkm1 232 ! 233 !--- Computation of the ustream and downstream value of the tracer and the mask 234 DO_2D( 0, 0, nn_hls-1, nn_hls-1 ) 235 ! Upstream in the x-direction for the tracer 236 zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb) 237 ! Downstream in the x-direction for the tracer 238 zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb) 239 END_2D 240 END DO 241 IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp ) ! Lateral boundary conditions 238 !--- Computation of the ustream and downstream value of the tracer and the mask 239 DO_3D( 0, 0, nn_hls-1, nn_hls-1, 1, jpkm1 ) 240 ! Upstream in the x-direction for the tracer 241 zfc(ji,jj,jk) = pt(ji,jj-1,jk,jn,Kbb) 242 ! Downstream in the x-direction for the tracer 243 zfd(ji,jj,jk) = pt(ji,jj+1,jk,jn,Kbb) 244 END_3D 245 246 IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfc(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. ) ! Lateral boundary conditions 242 247 243 248 ! … … 259 264 260 265 !--- Lateral boundary conditions 261 IF (nn_hls .EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp )266 IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp , zfd(:,:,:), 'T', 1.0_wp, zfc(:,:,:), 'T', 1.0_wp, zwy(:,:,:), 'T', 1.0_wp ) 262 267 263 268 !--- QUICKEST scheme … … 268 273 zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2. 269 274 END_3D 270 IF (nn_hls .EQ.1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp) !--- Lateral boundary conditions275 IF (nn_hls==1) CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1.0_wp, ld4only= .TRUE. ) !--- Lateral boundary conditions 271 276 ! 272 277 ! Tracer flux on the x-direction … … 306 311 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 307 312 INTEGER , INTENT(in ) :: kjpt ! number of tracers 308 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)313 ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 309 314 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pW ! vertical velocity 310 315 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! active tracers and RHS of tracer equation … … 365 370 !---------------------------------------------------------------------- 366 371 ! 367 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )372 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 368 373 zc = puc(ji,jj,jk) ! Courant number 369 374 zcurv = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk) -
NEMO/trunk/src/OCE/TRA/traadv_ubs.F90
r14433 r14834 26 26 USE lbclnk ! ocean lateral boundary condition (or mpp link) 27 27 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 28 #if defined key_loop_fusion 29 USE traadv_ubs_lf ! UBS scheme (tra_adv_ubs routine - loop fusion version) 30 #endif 28 31 29 32 IMPLICIT NONE … … 92 95 INTEGER , INTENT(in ) :: kn_ubs_v ! number of tracers 93 96 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 94 ! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)97 ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 95 98 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume transport components 96 99 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation … … 103 106 !!---------------------------------------------------------------------- 104 107 ! 105 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile 108 #if defined key_loop_fusion 109 CALL tra_adv_ubs_lf ( kt, kit000, cdtype, p2dt, pU, pV, pW, Kbb, Kmm, pt, kjpt, Krhs, kn_ubs_v ) 110 #else 111 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 106 112 IF( kt == kit000 ) THEN 107 113 IF(lwp) WRITE(numout,*) … … 140 146 ! 141 147 END DO 142 IF (nn_hls .EQ.1) CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1.0_wp, zltv, 'T', 1.0_wp) ! Lateral boundary cond. (unchanged sgn)148 IF (nn_hls==1) CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1.0_wp, zltv, 'T', 1.0_wp, ld4only= .TRUE. ) ! Lateral boundary cond. (unchanged sgn) 143 149 ! 144 150 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) !== Horizontal advective fluxes ==! (UBS) … … 260 266 END DO 261 267 ! 268 #endif 262 269 END SUBROUTINE tra_adv_ubs 263 270 -
NEMO/trunk/src/OCE/TRA/trabbc.F90
r14072 r14834 102 102 ENDIF 103 103 ! 104 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain 105 CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) ) 106 ENDIF 104 CALL iom_put ( "hfgeou" , rho0_rcp * qgh_trd0(:,:) ) 105 107 106 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbc - Ta: ', mask1=tmask, clinfo3='tra-ta' ) 108 107 ! -
NEMO/trunk/src/OCE/TRA/trabbl.F90
r14820 r14834 126 126 CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_ldf - Ta: ', mask1=tmask, & 127 127 & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 128 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 129 CALL iom_put( "ahu_bbl", ahu_bbl ) ! bbl diffusive flux i-coef 130 CALL iom_put( "ahv_bbl", ahv_bbl ) ! bbl diffusive flux j-coef 131 ENDIF 128 CALL iom_put( "ahu_bbl", ahu_bbl ) ! bbl diffusive flux i-coef 129 CALL iom_put( "ahv_bbl", ahv_bbl ) ! bbl diffusive flux j-coef 132 130 ! 133 131 ENDIF … … 139 137 CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv - Ta: ', mask1=tmask, & 140 138 & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 141 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 142 ! lateral boundary conditions ; just need for outputs 143 CALL lbc_lnk( 'trabbl', utr_bbl, 'U', -1.0_wp , vtr_bbl, 'V', -1.0_wp ) 144 CALL iom_put( "uoce_bbl", utr_bbl ) ! bbl i-transport 145 CALL iom_put( "voce_bbl", vtr_bbl ) ! bbl j-transport 146 ENDIF 139 CALL iom_put( "uoce_bbl", utr_bbl ) ! bbl i-transport 140 CALL iom_put( "voce_bbl", vtr_bbl ) ! bbl j-transport 147 141 ! 148 142 ENDIF … … 215 209 216 210 211 ! NOTE: [tiling] tiling changes the results, but only the order of floating point operations is different 217 212 SUBROUTINE tra_bbl_adv( pt, pt_rhs, kjpt, Kmm ) 218 213 !!---------------------------------------------------------------------- … … 238 233 INTEGER :: iis , iid , ijs , ijd ! local integers 239 234 INTEGER :: ikus, ikud, ikvs, ikvd ! - - 240 INTEGER :: isi, isj ! - -241 235 REAL(wp) :: zbtr, ztra ! local scalars 242 236 REAL(wp) :: zu_bbl, zv_bbl ! - - 243 237 !!---------------------------------------------------------------------- 244 !245 IF( ntsi == Nis0 ) THEN ; isi = 1 ; ELSE ; isi = 0 ; ENDIF ! Avoid double-counting when using tiling246 IF( ntsj == Njs0 ) THEN ; isj = 1 ; ELSE ; isj = 0 ; ENDIF247 238 ! ! =========== 248 239 DO jn = 1, kjpt ! tracer loop 249 240 ! ! =========== 250 DO_2D ( isi, 0, isj, 0 ) ! CAUTION start from i=1 to update i=2 when cyclic east-west241 DO_2D_OVR( 1, 0, 1, 0 ) ! CAUTION start from i=1 to update i=2 when cyclic east-west 251 242 IF( utr_bbl(ji,jj) /= 0.e0 ) THEN ! non-zero i-direction bbl advection 252 243 ! down-slope i/k-indices (deep) & up-slope i/k indices (shelf) … … 340 331 !!---------------------------------------------------------------------- 341 332 ! 342 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile333 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 343 334 IF( kt == kit000 ) THEN 344 335 IF(lwp) WRITE(numout,*) … … 363 354 IF( nn_bbl_ldf == 1 ) THEN ! diffusive bbl ! 364 355 ! !-------------------! 365 DO_2D ( 1, 0, 1, 0 ) ! (criteria for non zero flux: grad(rho).grad(h) < 0 )356 DO_2D_OVR( 1, 0, 1, 0 ) ! (criteria for non zero flux: grad(rho).grad(h) < 0 ) 366 357 ! ! i-direction 367 358 za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point … … 393 384 ! 394 385 CASE( 1 ) != use of upper velocity 395 DO_2D ( 1, 0, 1, 0 ) ! criteria: grad(rho).grad(h)<0 and grad(rho).grad(h)<0386 DO_2D_OVR( 1, 0, 1, 0 ) ! criteria: grad(rho).grad(h)<0 and grad(rho).grad(h)<0 396 387 ! ! i-direction 397 388 za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem) ! 2*(alpha,beta) at u-point … … 422 413 CASE( 2 ) != bbl velocity = F( delta rho ) 423 414 zgbbl = grav * rn_gambbl 424 DO_2D ( 1, 0, 1, 0 ) ! criteria: rho_up > rho_down415 DO_2D_OVR( 1, 0, 1, 0 ) ! criteria: rho_up > rho_down 425 416 ! ! i-direction 426 417 ! down-slope T-point i/k-index (deep) & up-slope T-point i/k-index (shelf) -
NEMO/trunk/src/OCE/TRA/traisf.F90
r14072 r14834 47 47 IF( ln_timing ) CALL timing_start('tra_isf') 48 48 ! 49 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile49 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 50 50 IF( kt == nit000 ) THEN 51 51 IF(lwp) WRITE(numout,*) … … 79 79 ! 80 80 IF ( ln_isfdebug ) THEN 81 IF( ntile == 0.OR. ntile == nijtile ) THEN ! Do only for the full domain81 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only for the full domain 82 82 CALL debug('tra_isf: pts(:,:,:,:,Krhs) T', pts(:,:,:,1,Krhs)) 83 83 CALL debug('tra_isf: pts(:,:,:,:,Krhs) S', pts(:,:,:,2,Krhs)) -
NEMO/trunk/src/OCE/TRA/traldf.F90
r14189 r14834 17 17 USE oce ! ocean dynamics and tracers 18 18 USE dom_oce ! ocean space and time domain 19 ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*)20 USE domtile21 19 USE phycst ! physical constants 22 20 USE ldftra ! lateral diffusion: eddy diffusivity & EIV coeff. … … 58 56 !! 59 57 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds 60 ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*)61 LOGICAL :: lskip62 58 !!---------------------------------------------------------------------- 63 59 ! 64 60 IF( ln_timing ) CALL timing_start('tra_ldf') 65 61 ! 66 lskip = .FALSE.67 68 62 IF( l_trdtra ) THEN !* Save ta and sa trends 69 63 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) … … 71 65 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 72 66 ENDIF 73 74 ! TEMP: [tiling] These changes not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*) 75 IF( nldf_tra == np_blp .OR. nldf_tra == np_blp_i .OR. nldf_tra == np_blp_it ) THEN 76 IF( ln_tile ) THEN 77 IF( ntile == 1 ) THEN 78 CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) 79 ELSE 80 lskip = .TRUE. 81 ENDIF 82 ENDIF 83 ENDIF 84 IF( .NOT. lskip ) THEN 85 ! 86 SELECT CASE ( nldf_tra ) !* compute lateral mixing trend and add it to the general trend 87 CASE ( np_lap ) ! laplacian: iso-level operator 88 CALL tra_ldf_lap ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, 1 ) 89 CASE ( np_lap_i ) ! laplacian: standard iso-neutral operator (Madec) 90 CALL tra_ldf_iso ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, 1 ) 91 CASE ( np_lap_it ) ! laplacian: triad iso-neutral operator (griffies) 92 CALL tra_ldf_triad( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, 1 ) 93 CASE ( np_blp , np_blp_i , np_blp_it ) ! bilaplacian: iso-level & iso-neutral operators 94 IF(nn_hls.EQ.2) CALL lbc_lnk( 'tra_ldf', pts(:,:,:,:,Kbb), 'T',1.) 95 CALL tra_ldf_blp ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, nldf_tra ) 96 END SELECT 97 ! 98 IF( l_trdtra ) THEN !* save the horizontal diffusive trends for further diagnostics 99 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 100 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 101 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_ldf, ztrdt ) 102 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_ldf, ztrds ) 103 DEALLOCATE( ztrdt, ztrds ) 104 ENDIF 105 106 ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from tra_ldf_blp, zps_hde*) 107 IF( ln_tile .AND. ntile == 0 ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 1 ) 67 ! 68 SELECT CASE ( nldf_tra ) !* compute lateral mixing trend and add it to the general trend 69 CASE ( np_lap ) ! laplacian: iso-level operator 70 CALL tra_ldf_lap ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, 1 ) 71 CASE ( np_lap_i ) ! laplacian: standard iso-neutral operator (Madec) 72 CALL tra_ldf_iso ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, 1 ) 73 CASE ( np_lap_it ) ! laplacian: triad iso-neutral operator (griffies) 74 CALL tra_ldf_triad( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, 1 ) 75 CASE ( np_blp , np_blp_i , np_blp_it ) ! bilaplacian: iso-level & iso-neutral operators 76 CALL tra_ldf_blp ( kt, Kmm, nit000,'TRA', ahtu, ahtv, gtsu, gtsv, gtui, gtvi, pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, nldf_tra ) 77 END SELECT 78 ! 79 IF( l_trdtra ) THEN !* save the horizontal diffusive trends for further diagnostics 80 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 81 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 82 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_ldf, ztrdt ) 83 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_ldf, ztrds ) 84 DEALLOCATE( ztrdt, ztrds ) 108 85 ENDIF 109 86 ! !* print mean trends (used for debugging) -
NEMO/trunk/src/OCE/TRA/traldf_iso.F90
r14820 r14834 132 132 INTEGER :: ji, jj, jk, jn ! dummy loop indices 133 133 INTEGER :: ikt 134 INTEGER :: ierr 134 INTEGER :: ierr, iij ! local integer 135 135 REAL(wp) :: zmsku, zahu_w, zabe1, zcof1, zcoef3 ! local scalars 136 136 REAL(wp) :: zmskv, zahv_w, zabe2, zcof2, zcoef4 ! - - … … 141 141 ! 142 142 IF( kpass == 1 .AND. kt == kit000 ) THEN 143 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile143 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 144 144 IF(lwp) WRITE(numout,*) 145 145 IF(lwp) WRITE(numout,*) 'tra_ldf_iso : rotated laplacian diffusion operator on ', cdtype … … 147 147 ENDIF 148 148 ! 149 DO_3D ( 0, 0, 0, 0, 1, jpk )149 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 150 150 akz (ji,jj,jk) = 0._wp 151 151 ah_wslp2(ji,jj,jk) = 0._wp … … 153 153 ENDIF 154 154 ! 155 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile155 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 156 156 l_hst = .FALSE. 157 157 l_ptr = .FALSE. … … 161 161 ENDIF 162 162 ! 163 ! Define pt_rhs halo points for multi-point haloes in bilaplacian case 164 IF( nldf_tra == np_blp_i .AND. kpass == 1 ) THEN ; iij = nn_hls 165 ELSE ; iij = 1 166 ENDIF 167 163 168 ! 164 169 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) … … 172 177 IF( kpass == 1 ) THEN !== first pass only ==! 173 178 ! 174 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )179 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 175 180 ! 176 181 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & … … 185 190 & + ( pahu(ji-1,jj,jk-1) + pahu(ji ,jj,jk) & 186 191 & ) & ! bracket for halo 1 - halo 2 compatibility 187 & ) * zmsku 192 & ) * zmsku 188 193 zahv_w = ( ( pahv(ji,jj ,jk-1) + pahv(ji,jj-1,jk) & 189 194 & ) & ! bracket for halo 1 - halo 2 compatibility 190 & + ( pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) & 195 & + ( pahv(ji,jj-1,jk-1) + pahv(ji,jj ,jk) & 191 196 & ) & ! bracket for halo 1 - halo 2 compatibility 192 & ) * zmskv 197 & ) * zmskv 193 198 ! 194 199 ah_wslp2(ji,jj,jk) = zahu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & … … 197 202 ! 198 203 IF( ln_traldf_msc ) THEN ! stabilizing vertical diffusivity coefficient 199 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )204 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 200 205 ! round brackets added to fix the order of floating point operations 201 206 ! needed to ensure halo 1 - halo 2 compatibility … … 205 210 & ) & ! bracket for halo 1 - halo 2 compatibility 206 211 & + ( ( pahv(ji,jj ,jk) + pahv(ji,jj ,jk-1) ) / ( e2v(ji,jj ) * e2v(ji,jj ) ) & 207 & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) & 212 & + ( pahv(ji,jj-1,jk) + pahv(ji,jj-1,jk-1) ) / ( e2v(ji,jj-1) * e2v(ji,jj-1) ) & 208 213 & ) & ! bracket for halo 1 - halo 2 compatibility 209 & ) 214 & ) 210 215 END_3D 211 216 ! 212 217 IF( ln_traldf_blp ) THEN ! bilaplacian operator 213 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )218 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 214 219 akz(ji,jj,jk) = 16._wp & 215 220 & * ah_wslp2 (ji,jj,jk) & … … 219 224 END_3D 220 225 ELSEIF( ln_traldf_lap ) THEN ! laplacian operator 221 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )226 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 222 227 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 223 228 zcoef0 = rDt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) … … 227 232 ! 228 233 ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 229 DO_3D ( 0, 0, 0, 0, 1, jpk )234 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 230 235 akz(ji,jj,jk) = ah_wslp2(ji,jj,jk) 231 236 END_3D … … 240 245 !! I - masked horizontal derivative 241 246 !!---------------------------------------------------------------------- 242 !!gm : bug.... why (x,:,:)? (1,jpj,:) and (jpi,1,:) should be sufficient.... 243 zdit (ntsi-nn_hls,:,:) = 0._wp ; zdit (ntei+nn_hls,:,:) = 0._wp 244 zdjt (ntsi-nn_hls,:,:) = 0._wp ; zdjt (ntei+nn_hls,:,:) = 0._wp 245 !!end 247 zdit(:,:,:) = 0._wp 248 zdjt(:,:,:) = 0._wp 246 249 247 250 ! Horizontal tracer gradient 248 DO_3D( 1, 0, 1, 0, 1, jpkm1 )251 DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 ) 249 252 zdit(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 250 253 zdjt(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 251 254 END_3D 252 255 IF( ln_zps ) THEN ! botton and surface ocean correction of the horizontal gradient 253 DO_2D( 1, 0, 1, 0 )! bottom correction (partial bottom cell)256 DO_2D( iij, iij-1, iij, iij-1 ) ! bottom correction (partial bottom cell) 254 257 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 255 258 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 256 259 END_2D 257 260 IF( ln_isfcav ) THEN ! first wet level beneath a cavity 258 DO_2D( 1, 0, 1, 0)261 DO_2D( iij, iij-1, iij, iij-1 ) 259 262 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj)) = pgui(ji,jj,jn) 260 263 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj)) = pgvi(ji,jj,jn) … … 269 272 DO jk = 1, jpkm1 ! Horizontal slab 270 273 ! 271 DO_2D( 1, 1, 1, 1)274 DO_2D( iij, iij, iij, iij ) 272 275 ! !== Vertical tracer gradient 273 276 zdk1t(ji,jj) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * wmask(ji,jj,jk+1) ! level jk+1 … … 278 281 END_2D 279 282 ! 280 DO_2D( 1, 0, 1, 0) !== Horizontal fluxes283 DO_2D( iij, iij-1, iij, iij-1 ) !== Horizontal fluxes 281 284 zabe1 = pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) 282 285 zabe2 = pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) … … 296 299 & + zcof1 * ( ( zdkt (ji+1,jj) + zdk1t(ji,jj) & 297 300 & ) & ! bracket for halo 1 - halo 2 compatibility 298 & + ( zdk1t(ji+1,jj) + zdkt (ji,jj) & 301 & + ( zdk1t(ji+1,jj) + zdkt (ji,jj) & 299 302 & ) & ! bracket for halo 1 - halo 2 compatibility 300 303 & ) ) * umask(ji,jj,jk) … … 304 307 & + ( zdk1t(ji,jj+1) + zdkt (ji,jj) & 305 308 & ) & ! bracket for halo 1 - halo 2 compatibility 306 & ) ) * vmask(ji,jj,jk) 309 & ) ) * vmask(ji,jj,jk) 307 310 END_2D 308 311 ! 309 DO_2D( 0, 0, 0, 0) !== horizontal divergence and add to pta312 DO_2D( iij-1, iij-1, iij-1, iij-1 ) !== horizontal divergence and add to pta 310 313 ! round brackets added to fix the order of floating point operations 311 314 ! needed to ensure halo 1 - halo 2 compatibility … … 328 331 ztfw(:,:, 1 ) = 0._wp ; ztfw(:,:,jpk) = 0._wp 329 332 330 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! interior (2=<jk=<jpk-1)333 DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) ! interior (2=<jk=<jpk-1) 331 334 ! 332 335 zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & … … 348 351 & ) & ! bracket for halo 1 - halo 2 compatibility 349 352 & + ( zdit(ji-1,jj ,jk-1) + zdit(ji ,jj ,jk) & 350 & ) & ! bracket for halo 1 - halo 2 compatibility 353 & ) & ! bracket for halo 1 - halo 2 compatibility 351 354 & ) & 352 355 & + zcoef4 * ( ( zdjt(ji ,jj ,jk-1) + zdjt(ji ,jj-1,jk) & … … 358 361 ! !== add the vertical 33 flux ==! 359 362 IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz 360 DO_3D( 0, 0, 0, 0, 2, jpkm1 )363 DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 361 364 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) + e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * wmask(ji,jj,jk) & 362 365 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & … … 367 370 SELECT CASE( kpass ) 368 371 CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2 369 DO_3D( 0, 0, 0, 0, 2, jpkm1 )372 DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 370 373 ztfw(ji,jj,jk) = & 371 374 & ztfw(ji,jj,jk) + ah_wslp2(ji,jj,jk) * e1e2t(ji,jj) & … … 381 384 ENDIF 382 385 ! 383 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Divergence of vertical fluxes added to pta ==!386 DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 ) !== Divergence of vertical fluxes added to pta ==! 384 387 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + zsign * ( ztfw (ji,jj,jk) - ztfw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) & 385 388 & / e3t(ji,jj,jk,Kmm) -
NEMO/trunk/src/OCE/TRA/traldf_lap_blp.F90
r14820 r14834 103 103 ! 104 104 INTEGER :: ji, jj, jk, jn ! dummy loop indices 105 INTEGER :: i si, iei, isj, iej ! local integers105 INTEGER :: iij 106 106 REAL(wp) :: zsign ! local scalars 107 107 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: ztu, ztv, zaheeu, zaheev 108 108 !!---------------------------------------------------------------------- 109 109 ! 110 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile110 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 111 111 IF( kt == nit000 .AND. lwp ) THEN 112 112 WRITE(numout,*) … … 122 122 ENDIF 123 123 ! 124 ! Define pt_rhs halo points for multi-point haloes in bilaplacian case 125 IF( nldf_tra == np_blp .AND. kpass == 1 ) THEN ; iij = nn_hls 126 ELSE ; iij = 1 127 ENDIF 128 124 129 ! !== Initialization of metric arrays used for all tracers ==! 125 130 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) … … 127 132 ENDIF 128 133 129 IF( ntsi == Nis0 ) THEN ; isi = nn_hls - 1 ; ELSE ; isi = 0 ; ENDIF ! Avoid double-counting when using tiling 130 IF( ntsj == Njs0 ) THEN ; isj = nn_hls - 1 ; ELSE ; isj = 0 ; ENDIF 131 IF( ntei == Nie0 ) THEN ; iei = nn_hls - 1 ; ELSE ; iei = 0 ; ENDIF 132 IF( ntej == Nje0 ) THEN ; iej = nn_hls - 1 ; ELSE ; iej = 0 ; ENDIF 133 134 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) !== First derivative (gradient) ==! 134 DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 ) !== First derivative (gradient) ==! 135 135 zaheeu(ji,jj,jk) = zsign * pahu(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) !!gm * umask(ji,jj,jk) pah masked! 136 136 zaheev(ji,jj,jk) = zsign * pahv(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) !!gm * vmask(ji,jj,jk) … … 141 141 ! ! =========== ! 142 142 ! 143 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 ) !== First derivative (gradient) ==!143 DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 ) !== First derivative (gradient) ==! 144 144 ztu(ji,jj,jk) = zaheeu(ji,jj,jk) * ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) 145 145 ztv(ji,jj,jk) = zaheev(ji,jj,jk) * ( pt(ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) 146 146 END_3D 147 147 IF( ln_zps ) THEN ! set gradient at bottom/top ocean level 148 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! bottom148 DO_2D( iij, iij-1, iij, iij-1 ) ! bottom 149 149 ztu(ji,jj,mbku(ji,jj)) = zaheeu(ji,jj,mbku(ji,jj)) * pgu(ji,jj,jn) 150 150 ztv(ji,jj,mbkv(ji,jj)) = zaheev(ji,jj,mbkv(ji,jj)) * pgv(ji,jj,jn) 151 151 END_2D 152 152 IF( ln_isfcav ) THEN ! top in ocean cavities only 153 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )153 DO_2D( iij, iij-1, iij, iij-1 ) 154 154 IF( miku(ji,jj) > 1 ) ztu(ji,jj,miku(ji,jj)) = zaheeu(ji,jj,miku(ji,jj)) * pgui(ji,jj,jn) 155 155 IF( mikv(ji,jj) > 1 ) ztv(ji,jj,mikv(ji,jj)) = zaheev(ji,jj,mikv(ji,jj)) * pgvi(ji,jj,jn) … … 158 158 ENDIF 159 159 ! 160 DO_3D( i si, iei, isj, iej, 1, jpkm1 ) !== Second derivative (divergence) added to the general tracer trends ==!160 DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 ) !== Second derivative (divergence) added to the general tracer trends ==! 161 161 ! round brackets added to fix the order of floating point operations 162 162 ! needed to ensure halo 1 - halo 2 compatibility … … 215 215 !!--------------------------------------------------------------------- 216 216 ! 217 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile217 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 218 218 IF( kt == kit000 .AND. lwp ) THEN 219 219 WRITE(numout,*) … … 239 239 END SELECT 240 240 ! 241 CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp ) ! Lateral boundary conditions (unchanged sign)241 IF (nn_hls==1) CALL lbc_lnk( 'traldf_lap_blp', zlap(:,:,:,:) , 'T', 1.0_wp ) ! Lateral boundary conditions (unchanged sign) 242 242 ! ! Partial top/bottom cell: GRADh( zlap ) 243 243 IF( ln_isfcav .AND. ln_zps ) THEN ; CALL zps_hde_isf( kt, Kmm, kjpt, zlap, zglu, zglv, zgui, zgvi ) ! both top & bottom -
NEMO/trunk/src/OCE/TRA/traldf_triad.F90
r14820 r14834 13 13 USE oce ! ocean dynamics and active tracers 14 14 USE dom_oce ! ocean space and time domain 15 ! TEMP: [tiling] This change not necessary if XIOS has subdomain support16 USE domtile17 15 USE domutl, ONLY : is_tile 18 16 USE phycst ! physical constants … … 109 107 REAL(wp), DIMENSION(A2D_T(ktt_rhs),JPK,KJPT), INTENT(inout) :: pt_rhs ! tracer trend 110 108 ! 111 INTEGER :: ji, jj, jk, jn, kp ! dummy loop indices109 INTEGER :: ji, jj, jk, jn, kp, iij ! dummy loop indices 112 110 REAL(wp) :: zcoef0, ze3w_2, zsign ! - - 113 111 ! … … 115 113 REAL(wp) :: ze1ur, ze2vr, ze3wr, zdxt, zdyt, zdzt, zdyt_jp1, ze3wr_jp1, zdzt_jp1, zah_slp1, zah_slp_jp1, zaei_slp_jp1 116 114 REAL(wp) :: zah_slp, zaei_slp, zdxt_ip1, ze3wr_ip1, zdzt_ip1, zah_slp_ip1, zaei_slp_ip1, zaei_slp1 117 REAL(wp), DIMENSION(A2D(nn_hls),0:1) :: zdkt3d ! vertical tracer gradient at 2 levels 118 REAL(wp), DIMENSION(A2D(nn_hls) ) :: z2d ! 2D workspace 119 REAL(wp), DIMENSION(A2D(nn_hls) ,jpk) :: zdit, zdjt, zftu, zftv, ztfw ! 3D - 120 ! TEMP: [tiling] This can be A2D(nn_hls) if XIOS has subdomain support 121 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw 115 REAL(wp), DIMENSION(A2D(nn_hls),0:1) :: zdkt3d ! vertical tracer gradient at 2 levels 116 REAL(wp), DIMENSION(A2D(nn_hls) ) :: z2d ! 2D workspace 117 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdit, zdjt, zftu, zftv, ztfw, zpsi_uw, zpsi_vw ! 3D - 122 118 !!---------------------------------------------------------------------- 123 119 ! 124 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile120 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 125 121 IF( kpass == 1 .AND. kt == kit000 ) THEN 126 122 IF(lwp) WRITE(numout,*) … … 138 134 ENDIF 139 135 ! 136 ! Define pt_rhs halo points for multi-point haloes in bilaplacian case 137 IF( nldf_tra == np_blp_it .AND. kpass == 1 ) THEN ; iij = nn_hls 138 ELSE ; iij = 1 139 ENDIF 140 141 ! 140 142 IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign (eddy diffusivity >0) 141 143 ELSE ; zsign = -1._wp … … 148 150 IF( kpass == 1 ) THEN !== first pass only and whatever the tracer is ==! 149 151 ! 150 DO_3D ( 0, 0, 0, 0, 1, jpk )152 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 151 153 akz (ji,jj,jk) = 0._wp 152 154 ah_wslp2(ji,jj,jk) = 0._wp … … 154 156 ! 155 157 DO kp = 0, 1 ! i-k triads 156 DO_3D ( 0, 0, 0, 0, 1, jpkm1 )158 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 157 159 ze3wr = 1._wp / e3w(ji,jj,jk+kp,Kmm) 158 160 zbu = e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm) … … 177 179 ! 178 180 DO kp = 0, 1 ! j-k triads 179 DO_3D ( 0, 0, 0, 0, 1, jpkm1 )181 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 180 182 ze3wr = 1.0_wp / e3w(ji,jj,jk+kp,Kmm) 181 183 zbv = e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm) … … 203 205 ! 204 206 IF( ln_traldf_blp ) THEN ! bilaplacian operator 205 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )207 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 206 208 akz(ji,jj,jk) = 16._wp & 207 209 & * ah_wslp2 (ji,jj,jk) & … … 211 213 END_3D 212 214 ELSEIF( ln_traldf_lap ) THEN ! laplacian operator 213 DO_3D ( 0, 0, 0, 0, 2, jpkm1 )215 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 ) 214 216 ze3w_2 = e3w(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm) 215 217 zcoef0 = rDt * ( akz(ji,jj,jk) + ah_wslp2(ji,jj,jk) / ze3w_2 ) … … 219 221 ! 220 222 ELSE ! 33 flux set to zero with akz=ah_wslp2 ==>> computed in full implicit 221 DO_3D ( 0, 0, 0, 0, 1, jpk )223 DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk ) 222 224 akz(ji,jj,jk) = ah_wslp2(ji,jj,jk) 223 225 END_3D 224 226 ENDIF 225 227 ! 226 ! TEMP: [tiling] These changes not necessary if XIOS has subdomain support 227 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain 228 IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN 229 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 ) 230 231 zpsi_uw(:,:,:) = 0._wp 232 zpsi_vw(:,:,:) = 0._wp 233 234 DO kp = 0, 1 235 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 236 ! round brackets added to fix the order of floating point operations 237 ! needed to ensure halo 1 - halo 2 compatibility [ NOT TESTED ] 238 zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 239 & + ( 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji,jj,jk,1,kp) & 240 & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+1,jj,jk,0,kp) & 241 & ) ! bracket for halo 1 - halo 2 compatibility 242 zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) & 243 & + ( 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj,jk,1,kp) & 244 & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+1,jk,0,kp) & 245 & ) ! bracket for halo 1 - halo 2 compatibility 246 END_3D 247 END DO 248 CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 249 250 IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = nijtile ) 251 ENDIF 228 IF( ln_ldfeiv_dia .AND. cdtype == 'TRA' ) THEN 229 zpsi_uw(:,:,:) = 0._wp 230 zpsi_vw(:,:,:) = 0._wp 231 232 DO kp = 0, 1 233 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 234 ! round brackets added to fix the order of floating point operations 235 ! needed to ensure halo 1 - halo 2 compatibility 236 zpsi_uw(ji,jj,jk+kp) = zpsi_uw(ji,jj,jk+kp) & 237 & + ( 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji,jj,jk,1,kp) & 238 & + 0.25_wp * aeiu(ji,jj,jk) * e2u(ji,jj) * triadi_g(ji+1,jj,jk,0,kp) & 239 & ) ! bracket for halo 1 - halo 2 compatibility 240 zpsi_vw(ji,jj,jk+kp) = zpsi_vw(ji,jj,jk+kp) & 241 & + ( 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj,jk,1,kp) & 242 & + 0.25_wp * aeiv(ji,jj,jk) * e1v(ji,jj) * triadj_g(ji,jj+1,jk,0,kp) & 243 & ) ! bracket for halo 1 - halo 2 compatibility 244 END_3D 245 END DO 246 CALL ldf_eiv_dia( zpsi_uw, zpsi_vw, Kmm ) 252 247 ENDIF 253 248 ! … … 262 257 zftu(:,:,:) = 0._wp 263 258 zftv(:,:,:) = 0._wp 264 ! 265 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) !== before lateral T & S gradients at T-level jk ==! 259 zdit(:,:,:) = 0._wp 260 zdjt(:,:,:) = 0._wp 261 ! 262 DO_3D( iij, iij-1, iij, iij-1, 1, jpkm1 ) !== before lateral T & S gradients at T-level jk ==! 266 263 zdit(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn) - pt(ji,jj,jk,jn) ) * umask(ji,jj,jk) 267 264 zdjt(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn) - pt(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 268 265 END_3D 269 266 IF( ln_zps .AND. l_grad_zps ) THEN ! partial steps: correction at top/bottom ocean level 270 DO_2D( 1, 0, 1, 0) ! bottom level267 DO_2D( iij, iij-1, iij, iij-1 ) ! bottom level 271 268 zdit(ji,jj,mbku(ji,jj)) = pgu(ji,jj,jn) 272 269 zdjt(ji,jj,mbkv(ji,jj)) = pgv(ji,jj,jn) 273 270 END_2D 274 271 IF( ln_isfcav ) THEN ! top level (ocean cavities only) 275 DO_2D( 1, 0, 1, 0)272 DO_2D( iij, iij-1, iij, iij-1 ) 276 273 IF( miku(ji,jj) > 1 ) zdit(ji,jj,miku(ji,jj) ) = pgui(ji,jj,jn) 277 274 IF( mikv(ji,jj) > 1 ) zdjt(ji,jj,mikv(ji,jj) ) = pgvi(ji,jj,jn) … … 286 283 DO jk = 1, jpkm1 287 284 ! !== Vertical tracer gradient at level jk and jk+1 288 DO_2D( 1, 1, 1, 1)285 DO_2D( iij, iij, iij, iij ) 289 286 zdkt3d(ji,jj,1) = ( pt(ji,jj,jk,jn) - pt(ji,jj,jk+1,jn) ) * tmask(ji,jj,jk+1) 290 287 END_2D … … 293 290 IF( jk == 1 ) THEN ; zdkt3d(:,:,0) = zdkt3d(:,:,1) 294 291 ELSE 295 DO_2D( 1, 1, 1, 1)292 DO_2D( iij, iij, iij, iij ) 296 293 zdkt3d(ji,jj,0) = ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 297 294 END_2D … … 305 302 IF( ln_botmix_triad ) THEN 306 303 DO kp = 0, 1 !== Horizontal & vertical fluxes 307 DO_2D( 1, 0, 1, 0)304 DO_2D( iij, iij-1, iij, iij-1 ) 308 305 ze1ur = r1_e1u(ji,jj) 309 306 zdxt = zdit(ji,jj,jk) * ze1ur … … 317 314 zbu_ip1 = 0.25_wp * e1e2u(ji+1,jj) * e3u(ji+1,jj,jk,Kmm) 318 315 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked.... 319 zah = pahu(ji,jj,jk) 320 zah_ip1 = pahu(ji+1,jj,jk) 316 zah = pahu(ji,jj,jk) 317 zah_ip1 = pahu(ji+1,jj,jk) 321 318 zah_slp = zah * triadi(ji,jj,jk,1,kp) 322 319 zah_slp_ip1 = zah_ip1 * triadi(ji+1,jj,jk,1,kp) … … 333 330 & + ( zah * zdxt + zah_slp1 * zdzt_ip1 - zaei_slp1 * zdzt_ip1 ) * zbu * ze1ur & 334 331 & ) ! bracket for halo 1 - halo 2 compatibility 335 ztfw(ji+1,jj,jk+kp) = ztfw(ji+1,jj,jk+kp) & 332 ztfw(ji+1,jj,jk+kp) = ztfw(ji+1,jj,jk+kp) & 336 333 & - ( (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1 & 337 & + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1 & 338 & ) ! bracket for halo 1 - halo 2 compatibility 339 340 334 & + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1 & 335 & ) ! bracket for halo 1 - halo 2 compatibility 336 END_2D 337 END DO 341 338 ! 342 339 DO kp = 0, 1 343 DO_2D( 1, 0, 1, 0)340 DO_2D( iij, iij-1, iij, iij-1 ) 344 341 ze2vr = r1_e2v(ji,jj) 345 342 zdyt = zdjt(ji,jj,jk) * ze2vr … … 353 350 ! ln_botmix_triad is .T. don't mask zah for bottom half cells !!gm ????? ahu is masked.... 354 351 zah = pahv(ji,jj,jk) ! pahv(ji,jj+jp,jk) ???? 355 zah_jp1 = pahv(ji,jj+1,jk) 352 zah_jp1 = pahv(ji,jj+1,jk) 356 353 zah_slp = zah * triadj(ji,jj,jk,1,kp) 357 354 zah_slp1 = zah * triadj(ji,jj+1,jk,0,kp) 358 355 zah_slp_jp1 = zah_jp1 * triadj(ji,jj+1,jk,1,kp) 359 356 IF( ln_ldfeiv ) THEN 360 zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp) 361 zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp) 357 zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp) 358 zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp) 362 359 zaei_slp1 = aeiv(ji,jj,jk) * triadj_g(ji,jj+1,jk,0,kp) 363 360 ENDIF … … 377 374 ELSE 378 375 ! 379 DO kp = 0, 1 380 DO_2D( 1, 0, 1, 0)376 DO kp = 0, 1 !== Horizontal & vertical fluxes 377 DO_2D( iij, iij-1, iij, iij-1 ) 381 378 ze1ur = r1_e1u(ji,jj) 382 379 zdxt = zdit(ji,jj,jk) * ze1ur … … 391 388 ! ln_botmix_triad is .F. mask zah for bottom half cells 392 389 zah = pahu(ji,jj,jk) * umask(ji,jj,jk+kp) ! pahu(ji+ip,jj,jk) ===>> ???? 393 zah_ip1 = pahu(ji+1,jj,jk) * umask(ji+1,jj,jk+kp) 390 zah_ip1 = pahu(ji+1,jj,jk) * umask(ji+1,jj,jk+kp) 394 391 zah_slp = zah * triadi(ji,jj,jk,1,kp) 395 392 zah_slp_ip1 = zah_ip1 * triadi(ji+1,jj,jk,1,kp) 396 393 zah_slp1 = zah * triadi(ji+1,jj,jk,0,kp) 397 394 IF( ln_ldfeiv ) THEN 398 zaei_slp = aeiu(ji,jj,jk) * triadi_g(ji,jj,jk,1,kp) 395 zaei_slp = aeiu(ji,jj,jk) * triadi_g(ji,jj,jk,1,kp) 399 396 zaei_slp_ip1 = aeiu(ji+1,jj,jk) * triadi_g(ji+1,jj,jk,1,kp) 400 397 zaei_slp1 = aeiu(ji,jj,jk) * triadi_g(ji+1,jj,jk,0,kp) … … 406 403 & + ( zah * zdxt + zah_slp1 * zdzt_ip1 - zaei_slp1 * zdzt_ip1 ) * zbu * ze1ur & 407 404 & ) ! bracket for halo 1 - halo 2 compatibility 408 ztfw(ji+1,jj,jk+kp) = ztfw(ji+1,jj,jk+kp) & 405 ztfw(ji+1,jj,jk+kp) = ztfw(ji+1,jj,jk+kp) & 409 406 & - ( (zah_slp_ip1 + zaei_slp_ip1) * zdxt_ip1 * zbu_ip1 * ze3wr_ip1 & 410 407 & + ( zah_slp1 + zaei_slp1) * zdxt * zbu * ze3wr_ip1 & … … 414 411 ! 415 412 DO kp = 0, 1 416 DO_2D( 1, 0, 1, 0)413 DO_2D( iij, iij-1, iij, iij-1 ) 417 414 ze2vr = r1_e2v(ji,jj) 418 415 zdyt = zdjt(ji,jj,jk) * ze2vr … … 431 428 zah_slp_jp1 = zah_jp1 * triadj(ji,jj+1,jk,1,kp) 432 429 IF( ln_ldfeiv ) THEN 433 zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp) 430 zaei_slp = aeiv(ji,jj,jk) * triadj_g(ji,jj,jk,1,kp) 434 431 zaei_slp_jp1 = aeiv(ji,jj+1,jk) * triadj_g(ji,jj+1,jk,1,kp) 435 432 zaei_slp1 = aeiv(ji,jj,jk) * triadj_g(ji,jj+1,jk,0,kp) … … 449 446 ENDIF 450 447 ! !== horizontal divergence and add to the general trend ==! 451 DO_2D( 0, 0, 0, 0)448 DO_2D( iij-1, iij-1, iij-1, iij-1 ) 452 449 ! round brackets added to fix the order of floating point operations 453 450 ! needed to ensure halo 1 - halo 2 compatibility … … 464 461 ! !== add the vertical 33 flux ==! 465 462 IF( ln_traldf_lap ) THEN ! laplacian case: eddy coef = ah_wslp2 - akz 466 DO_3D( 0, 0, 1, 0, 2, jpkm1 )463 DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 467 464 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) & 468 465 & * ( ah_wslp2(ji,jj,jk) - akz(ji,jj,jk) ) & … … 472 469 SELECT CASE( kpass ) 473 470 CASE( 1 ) ! 1st pass : eddy coef = ah_wslp2 474 DO_3D( 0, 0, 1, 0, 2, jpkm1 )471 DO_3D( iij-1, iij-1, iij-1, iij-1, 2, jpkm1 ) 475 472 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) & 476 473 & * ah_wslp2(ji,jj,jk) * ( pt(ji,jj,jk-1,jn) - pt(ji,jj,jk,jn) ) 477 474 END_3D 478 475 CASE( 2 ) ! 2nd pass : eddy flux = ah_wslp2 and akz applied on pt and pt2 gradients, resp. 479 DO_3D( 0, 0, 1, 0, 2, jpkm1 )476 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 480 477 ztfw(ji,jj,jk) = ztfw(ji,jj,jk) - e1e2t(ji,jj) / e3w(ji,jj,jk,Kmm) * tmask(ji,jj,jk) & 481 478 & * ( ah_wslp2(ji,jj,jk) * ( pt (ji,jj,jk-1,jn) - pt (ji,jj,jk,jn) ) & … … 485 482 ENDIF 486 483 ! 487 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Divergence of vertical fluxes added to pta ==!484 DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 ) !== Divergence of vertical fluxes added to pta ==! 488 485 pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) & 489 486 & + zsign * ( ztfw(ji,jj,jk+1) - ztfw(ji,jj,jk) ) & -
NEMO/trunk/src/OCE/TRA/tramle.F90
r14433 r14834 87 87 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 88 88 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 89 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: pu ! in : 3 ocean transport components 90 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: pv ! out: same 3 transport components 91 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) :: pw ! increased by the MLE induced transport 89 ! TEMP: [tiling] Can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct 90 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu ! in : 3 ocean transport components 91 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pv ! out: same 3 transport components 92 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pw ! increased by the MLE induced transport 92 93 ! 93 94 INTEGER :: ji, jj, jk ! dummy loop indices … … 96 97 REAL(wp) :: zcvw, zmvw ! - - 97 98 INTEGER , DIMENSION(A2D(nn_hls)) :: inml_mle 98 REAL(wp), DIMENSION(A2D(nn_hls)) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_ MH99 REAL(wp), DIMENSION(A2D(nn_hls)) :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH 99 100 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpsi_uw, zpsi_vw 100 ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support)101 REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: zLf_NH102 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zpsiu_mle, zpsiv_mle103 101 !!---------------------------------------------------------------------- 104 102 ! … … 110 108 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 111 109 CASE ( 0 ) != min of the 2 neighbour MLDs 112 DO_2D( 1, 0, 1, 0)110 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 113 111 zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) ) 114 112 zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) ) 115 113 END_2D 116 114 CASE ( 1 ) != average of the 2 neighbour MLDs 117 DO_2D( 1, 0, 1, 0)115 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 118 116 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 119 117 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) 120 118 END_2D 121 119 CASE ( 2 ) != max of the 2 neighbour MLDs 122 DO_2D( 1, 0, 1, 0)120 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 123 121 zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) ) 124 122 zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) ) … … 126 124 END SELECT 127 125 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 128 DO_2D( 1, 0, 1, 0)126 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 129 127 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 130 128 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) & … … 137 135 ! 138 136 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 139 DO_2D( 1, 0, 1, 0)137 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 140 138 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2u(ji,jj) & 141 139 & * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) … … 149 147 ! !== MLD used for MLE ==! 150 148 ! ! compute from the 10m density to deal with the diurnal cycle 151 DO_2D( 1, 1, 1, 1)149 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 152 150 inml_mle(ji,jj) = mbkt(ji,jj) + 1 ! init. to number of ocean w-level (T-level + 1) 153 151 END_2D 154 152 IF ( nla10 > 0 ) THEN ! avoid case where first level is thicker than 10m 155 DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 ) ! from the bottom to nlb10 (10m)153 DO_3DS( nn_hls, nn_hls, nn_hls, nn_hls, jpkm1, nlb10, -1 ) ! from the bottom to nlb10 (10m) 156 154 IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle ) inml_mle(ji,jj) = jk ! Mixed layer 157 155 END_3D … … 163 161 zbm (:,:) = 0._wp 164 162 zn2 (:,:) = 0._wp 165 DO_3D( 1, 1, 1, 1, 1, ikmax ) ! MLD and mean buoyancy and N2 over the mixed layer163 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, ikmax ) ! MLD and mean buoyancy and N2 over the mixed layer 166 164 zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points 167 165 zmld(ji,jj) = zmld(ji,jj) + zc … … 172 170 SELECT CASE( nn_mld_uv ) ! MLD at u- & v-pts 173 171 CASE ( 0 ) != min of the 2 neighbour MLDs 174 DO_2D( 1, 0, 1, 0)172 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 175 173 zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) ) 176 174 zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) ) 177 175 END_2D 178 176 CASE ( 1 ) != average of the 2 neighbour MLDs 179 DO_2D( 1, 0, 1, 0)177 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 180 178 zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp 181 179 zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp 182 180 END_2D 183 181 CASE ( 2 ) != max of the 2 neighbour MLDs 184 DO_2D( 1, 0, 1, 0)182 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 185 183 zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) ) 186 184 zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) ) … … 188 186 END SELECT 189 187 ! ! convert density into buoyancy 190 DO_2D( 1, 1, 1, 1)188 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 191 189 zbm(ji,jj) = + grav * zbm(ji,jj) / MAX( e3t(ji,jj,1,Kmm), zmld(ji,jj) ) 192 190 END_2D … … 201 199 ! 202 200 IF( nn_mle == 0 ) THEN ! Fox-Kemper et al. 2010 formulation 203 DO_2D( 1, 0, 1, 0)201 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 204 202 zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 205 203 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) & … … 212 210 ! 213 211 ELSEIF( nn_mle == 1 ) THEN ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat) 214 DO_2D( 1, 0, 1, 0)212 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 215 213 zpsim_u(ji,jj) = rc_f * zhu(ji,jj) * zhu(ji,jj) * e2_e1u(ji,jj) & 216 214 & * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) ) … … 222 220 ! 223 221 IF( nn_conv == 1 ) THEN ! No MLE in case of convection 224 DO_2D( 1, 0, 1, 0)222 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 225 223 IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp ) zpsim_u(ji,jj) = 0._wp 226 224 IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp ) zpsim_v(ji,jj) = 0._wp … … 230 228 ENDIF ! end of ln_osm_mle conditional 231 229 ! !== structure function value at uw- and vw-points ==! 232 DO_2D( 1, 0, 1, 0)230 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 233 231 zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall) ! hu --> 1/hu 234 232 zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall) … … 238 236 zpsi_vw(:,:,:) = 0._wp 239 237 ! 240 DO_3D( 1, 0, 1, 0, 2, ikmax ) ! start from 2 : surface value = 0 238 DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 2, ikmax ) ! start from 2 : surface value = 0 239 241 240 zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj) 242 241 zcvw = 1._wp - ( gdepw(ji,jj+1,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhv(ji,jj) … … 252 251 ! !== transport increased by the MLE induced transport ==! 253 252 DO jk = 1, ikmax 254 DO_2D ( 1, 0, 1, 0 ) ! CAUTION pu,pv must be defined at row/column i=1 / j=1253 DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 255 254 pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) ) 256 255 pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) ) 257 256 END_2D 258 DO_2D ( 0, 0, 0, 0)257 DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 259 258 pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk) & 260 259 & + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) ) * wmask(ji,jj,1) … … 262 261 END DO 263 262 264 ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support)265 263 IF( cdtype == 'TRA') THEN !== outputs ==! 266 IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile267 ALLOCATE( zLf_NH(jpi,jpj), zpsiu_mle(jpi,jpj,jpk), zpsiv_mle(jpi,jpj,jpk) )268 zpsiu_mle(:,:,:) = 0._wp ; zpsiv_mle(:,:,:) = 0._wp269 ENDIF270 264 ! 271 265 IF (ln_osm_mle.and.ln_zdfosm) THEN … … 279 273 ENDIF 280 274 ! 275 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 276 ! 281 277 ! divide by cross distance to give streamfunction with dimensions m^2/s 282 278 DO_3D( 0, 0, 0, 0, 1, ikmax+1 ) 283 zpsi u_mle(ji,jj,jk) = zpsi_uw(ji,jj,jk) * r1_e2u(ji,jj)284 zpsi v_mle(ji,jj,jk) = zpsi_vw(ji,jj,jk) * r1_e1v(ji,jj)279 zpsi_uw(ji,jj,jk) = zpsi_uw(ji,jj,jk) * r1_e2u(ji,jj) 280 zpsi_vw(ji,jj,jk) = zpsi_vw(ji,jj,jk) * r1_e1v(ji,jj) 285 281 END_3D 286 287 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 288 CALL iom_put( "Lf_NHpf" , zLf_NH ) ! Lf = N H / f 289 CALL iom_put( "psiu_mle", zpsiu_mle ) ! i-mle streamfunction 290 CALL iom_put( "psiv_mle", zpsiv_mle ) ! j-mle streamfunction 291 DEALLOCATE( zLf_NH, zpsiu_mle, zpsiv_mle ) 292 ENDIF 282 CALL iom_put( "psiu_mle", zpsi_uw ) ! i-mle streamfunction 283 CALL iom_put( "psiv_mle", zpsi_vw ) ! j-mle streamfunction 293 284 ENDIF 294 285 ! … … 375 366 r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 ) 376 367 ! 368 ! Specifically, dbdx_mle, dbdy_mle and mld_prof need to be defined for nn_hls = 2 369 IF( nn_hls == 2 .AND. ln_osm_mle .AND. ln_zdfosm ) THEN 370 CALL ctl_stop('nn_hls = 2 cannot be used with ln_mle = ln_osm_mle = ln_zdfosm = T (zdfosm not updated for nn_hls = 2)') 371 ENDIF 377 372 ENDIF 378 373 ! -
NEMO/trunk/src/OCE/TRA/tranpc.F90
r14215 r14834 17 17 USE oce ! ocean dynamics and active tracers 18 18 USE dom_oce ! ocean space and time domain 19 ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed)20 USE domtile21 19 USE phycst ! physical constants 22 20 USE zdf_oce ! ocean vertical physics … … 81 79 LOGICAL, PARAMETER :: l_LB_debug = .FALSE. ! set to true if you want to follow what is 82 80 INTEGER :: ilc1, jlc1, klc1, nncpu ! actually happening in a water column at point "ilc1, jlc1" 83 INTEGER :: isi, isj, iei, iej84 81 LOGICAL :: lp_monitor_point = .FALSE. ! in CPU domain "nncpu" 85 82 !!---------------------------------------------------------------------- … … 105 102 CALL bn2 ( pts(:,:,:,:,Kaa), zab, zn2, Kmm ) ! after Brunt-Vaisala (given on W-points) 106 103 ! 107 IF( ntile == 0 .OR. ntile == 1 ) nnpcc = 0 ! Do only on the first tile 108 ! 109 IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF ! Avoid double-counting when using tiling 110 IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF 111 IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF 112 IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF 113 ! 114 DO_2D( isi, iei, isj, iej ) ! interior column only 104 IF( .NOT. l_istiled .OR. ntile == 1 ) nnpcc = 0 ! Do only on the first tile 105 ! 106 DO_2D_OVR( 0, 0, 0, 0 ) ! interior column only 115 107 ! 116 108 IF( tmask(ji,jj,2) == 1 ) THEN ! At least 2 ocean points … … 319 311 ENDIF 320 312 ! 321 IF( ntile == 0.OR. ntile == nijtile ) THEN ! Do only for the full domain313 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only for the full domain 322 314 IF( lwp .AND. l_LB_debug ) THEN 323 315 WRITE(numout,*) 'Exiting tra_npc , kt = ',kt,', => numb. of statically instable water-columns: ', nnpcc -
NEMO/trunk/src/OCE/TRA/traqsr.F90
r14215 r14834 108 108 ! 109 109 INTEGER :: ji, jj, jk ! dummy loop indices 110 INTEGER :: irgb , isi, iei, isj, iej! local integers110 INTEGER :: irgb ! local integers 111 111 REAL(wp) :: zchl, zcoef, z1_2 ! local scalars 112 112 REAL(wp) :: zc0 , zc1 , zc2 , zc3 ! - - … … 121 121 IF( ln_timing ) CALL timing_start('tra_qsr') 122 122 ! 123 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile123 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 124 124 IF( kt == nit000 ) THEN 125 125 IF(lwp) WRITE(numout,*) … … 137 137 ! ! before qsr induced heat content ! 138 138 ! !-----------------------------------! 139 IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF ! Avoid double-counting when using tiling140 IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF141 IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF142 IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF143 144 139 IF( kt == nit000 ) THEN !== 1st time step ==! 145 140 IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN ! read in restart 146 141 z1_2 = 0.5_wp 147 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile142 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 148 143 IF(lwp) WRITE(numout,*) ' nit000-1 qsr tracer content forcing field read in the restart file' 149 144 CALL iom_get( numror, jpdom_auto, 'qsr_hc_b', qsr_hc_b ) ! before heat content trend due to Qsr flux … … 151 146 ELSE ! No restart or Euler forward at 1st time step 152 147 z1_2 = 1._wp 153 DO_3D ( isi, iei, isj, iej, 1, jpk )148 DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 154 149 qsr_hc_b(ji,jj,jk) = 0._wp 155 150 END_3D … … 157 152 ELSE !== Swap of qsr heat content ==! 158 153 z1_2 = 0.5_wp 159 DO_3D ( isi, iei, isj, iej, 1, jpk )154 DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk ) 160 155 qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk) 161 156 END_3D … … 168 163 CASE( np_BIO ) !== bio-model fluxes ==! 169 164 ! 170 DO_3D ( isi, iei, isj, iej, 1, nksr )165 DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr ) 171 166 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( etot3(ji,jj,jk) - etot3(ji,jj,jk+1) ) 172 167 END_3D … … 179 174 ! 180 175 IF( nqsr == np_RGBc ) THEN !* Variable Chlorophyll 181 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only for the full domain182 IF( ln_tile ) CALL dom_tile ( ntsi, ntsj, ntei, ntej, ktile = 0 )! Use full domain176 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only for the full domain 177 IF( ln_tile ) CALL dom_tile_stop( ldhold=.TRUE. ) ! Use full domain 183 178 CALL fld_read( kt, 1, sf_chl ) ! Read Chl data and provides it at the current time step 184 IF( ln_tile ) CALL dom_tile ( ntsi, ntsj, ntei, ntej, ktile = 1) ! Revert to tile domain179 IF( ln_tile ) CALL dom_tile_start( ldhold=.TRUE. ) ! Revert to tile domain 185 180 ENDIF 186 181 ! … … 190 185 ! most expensive calculations) 191 186 ! 192 DO_2D ( isi, iei, isj, iej)187 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 193 188 ! zlogc = log(zchl) 194 189 zlogc = LOG ( MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) ) … … 209 204 210 205 ! 211 DO_3D ( isi, iei, isj, iej, 1, nksr + 1 )206 DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr + 1 ) 212 207 ! zchl = ALOG( ze0(ji,jj) ) 213 208 zlogc = ze0(ji,jj) … … 239 234 ! 240 235 zcoef = ( 1. - rn_abs ) / 3._wp !* surface equi-partition in R-G-B 241 DO_2D ( isi, iei, isj, iej)236 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 242 237 ze0(ji,jj) = rn_abs * qsr(ji,jj) 243 238 ze1(ji,jj) = zcoef * qsr(ji,jj) … … 250 245 ! 251 246 ! !* interior equi-partition in R-G-B depending on vertical profile of Chl 252 DO_3D ( isi, iei, isj, iej, 2, nksr + 1 )247 DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 2, nksr + 1 ) 253 248 ze3t = e3t(ji,jj,jk-1,Kmm) 254 249 irgb = NINT( ztmp3d(ji,jj,jk) ) … … 264 259 END_3D 265 260 ! 266 DO_3D ( isi, iei, isj, iej, 1, nksr ) !* now qsr induced heat content261 DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr ) !* now qsr induced heat content 267 262 qsr_hc(ji,jj,jk) = r1_rho0_rcp * ( ztmp3d(ji,jj,jk) - ztmp3d(ji,jj,jk+1) ) 268 263 END_3D … … 274 269 zz0 = rn_abs * r1_rho0_rcp ! surface equi-partition in 2-bands 275 270 zz1 = ( 1. - rn_abs ) * r1_rho0_rcp 276 DO_3D ( isi, iei, isj, iej, 1, nksr ) !* now qsr induced heat content271 DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, nksr ) !* now qsr induced heat content 277 272 zc0 = zz0 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk ,Kmm)*xsi1r ) 278 273 zc1 = zz0 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi0r ) + zz1 * EXP( -gdepw(ji,jj,jk+1,Kmm)*xsi1r ) … … 292 287 ! 293 288 ! sea-ice: store the 1st ocean level attenuation coefficient 294 DO_2D ( isi, iei, isj, iej)289 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 295 290 IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) ) 296 291 ELSE ; fraqsr_1lev(ji,jj) = 1._wp … … 298 293 END_2D 299 294 ! 300 ! TEMP: [tiling] This change not necessary and working array can use A2D(nn_hls) if using XIOS (subdomain support) 301 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only for the full domain 302 IF( iom_use('qsr3d') ) THEN ! output the shortwave Radiation distribution 303 ALLOCATE( zetot(jpi,jpj,jpk) ) 304 zetot(:,:,nksr+1:jpk) = 0._wp ! below ~400m set to zero 305 DO jk = nksr, 1, -1 306 zetot(:,:,jk) = zetot(:,:,jk+1) + qsr_hc(:,:,jk) * rho0_rcp 307 END DO 308 CALL iom_put( 'qsr3d', zetot ) ! 3D distribution of shortwave Radiation 309 DEALLOCATE( zetot ) 310 ENDIF 311 ENDIF 312 ! 313 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 295 IF( iom_use('qsr3d') ) THEN ! output the shortwave Radiation distribution 296 ALLOCATE( zetot(A2D(nn_hls),jpk) ) 297 zetot(:,:,nksr+1:jpk) = 0._wp ! below ~400m set to zero 298 DO_3DS(0, 0, 0, 0, nksr, 1, -1) 299 zetot(ji,jj,jk) = zetot(ji,jj,jk+1) + qsr_hc(ji,jj,jk) * rho0_rcp 300 END_3D 301 CALL iom_put( 'qsr3d', zetot ) ! 3D distribution of shortwave Radiation 302 DEALLOCATE( zetot ) 303 ENDIF 304 ! 305 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 314 306 IF( lrst_oce ) THEN ! write in the ocean restart file 315 307 CALL iom_rstput( kt, nitrst, numrow, 'qsr_hc_b' , qsr_hc ) -
NEMO/trunk/src/OCE/TRA/trasbc.F90
r14215 r14834 77 77 ! 78 78 INTEGER :: ji, jj, jk, jn ! dummy loop indices 79 INTEGER :: ikt, ikb , isi, iei, isj, iej! local integers79 INTEGER :: ikt, ikb ! local integers 80 80 REAL(wp) :: zfact, z1_e3t, zdep, ztim ! local scalar 81 81 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdt, ztrds … … 84 84 IF( ln_timing ) CALL timing_start('tra_sbc') 85 85 ! 86 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile86 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 87 87 IF( kt == nit000 ) THEN 88 88 IF(lwp) WRITE(numout,*) … … 98 98 ENDIF 99 99 ! 100 IF( ntsi == Nis0 ) THEN ; isi = nn_hls ; ELSE ; isi = 0 ; ENDIF ! Avoid double-counting when using tiling101 IF( ntsj == Njs0 ) THEN ; isj = nn_hls ; ELSE ; isj = 0 ; ENDIF102 IF( ntei == Nie0 ) THEN ; iei = nn_hls ; ELSE ; iei = 0 ; ENDIF103 IF( ntej == Nje0 ) THEN ; iej = nn_hls ; ELSE ; iej = 0 ; ENDIF104 105 100 !!gm This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist) 106 101 IF( .NOT.ln_traqsr ) THEN ! no solar radiation penetration 107 DO_2D ( isi, iei, isj, iej)102 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 108 103 qns(ji,jj) = qns(ji,jj) + qsr(ji,jj) ! total heat flux in qns 109 104 qsr(ji,jj) = 0._wp ! qsr set to zero … … 118 113 IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN ! Restart: read in restart file 119 114 zfact = 0.5_wp 120 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile115 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 121 116 IF(lwp) WRITE(numout,*) ' nit000-1 sbc tracer content field read in the restart file' 122 117 sbc_tsc(:,:,:) = 0._wp … … 126 121 ELSE ! No restart or restart not found: Euler forward time stepping 127 122 zfact = 1._wp 128 DO_2D ( isi, iei, isj, iej)123 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 129 124 sbc_tsc(ji,jj,:) = 0._wp 130 125 sbc_tsc_b(ji,jj,:) = 0._wp … … 133 128 ELSE !* other time-steps: swap of forcing fields 134 129 zfact = 0.5_wp 135 DO_2D ( isi, iei, isj, iej)130 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 136 131 sbc_tsc_b(ji,jj,:) = sbc_tsc(ji,jj,:) 137 132 END_2D 138 133 ENDIF 139 134 ! !== Now sbc tracer content fields ==! 140 DO_2D ( isi, iei, isj, iej)135 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) 141 136 sbc_tsc(ji,jj,jp_tem) = r1_rho0_rcp * qns(ji,jj) ! non solar heat flux 142 137 sbc_tsc(ji,jj,jp_sal) = r1_rho0 * sfx(ji,jj) ! salt flux due to freezing/melting 143 138 END_2D 144 139 IF( ln_linssh ) THEN !* linear free surface 145 DO_2D ( isi, iei, isj, iej) !==>> add concentration/dilution effect due to constant volume cell140 DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls ) !==>> add concentration/dilution effect due to constant volume cell 146 141 sbc_tsc(ji,jj,jp_tem) = sbc_tsc(ji,jj,jp_tem) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_tem,Kmm) 147 142 sbc_tsc(ji,jj,jp_sal) = sbc_tsc(ji,jj,jp_sal) + r1_rho0 * emp(ji,jj) * pts(ji,jj,1,jp_sal,Kmm) 148 143 END_2D !==>> output c./d. term 149 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 150 IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) ) 151 IF( iom_use('emp_x_sss') ) CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) ) 152 ENDIF 144 IF( iom_use('emp_x_sst') ) CALL iom_put( "emp_x_sst", emp (:,:) * pts(:,:,1,jp_tem,Kmm) ) 145 IF( iom_use('emp_x_sss') ) CALL iom_put( "emp_x_sss", emp (:,:) * pts(:,:,1,jp_sal,Kmm) ) 153 146 ENDIF 154 147 ! … … 160 153 END DO 161 154 ! 162 IF( ntile == 0.OR. ntile == nijtile ) THEN ! Do only on the last tile155 IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN ! Do only on the last tile 163 156 IF( lrst_oce ) THEN !== write sbc_tsc in the ocean restart file ==! 164 157 CALL iom_rstput( kt, nitrst, numrow, 'sbc_hc_b', sbc_tsc(:,:,jp_tem) ) … … 186 179 ENDIF 187 180 188 IF( ntile == 0 .OR. ntile == nijtile ) THEN ! Do only on the last tile 189 IF( iom_use('rnf_x_sst') ) CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) ) ! runoff term on sst 190 IF( iom_use('rnf_x_sss') ) CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) ) ! runoff term on sss 191 ENDIF 181 IF( iom_use('rnf_x_sst') ) CALL iom_put( "rnf_x_sst", rnf*pts(:,:,1,jp_tem,Kmm) ) ! runoff term on sst 182 IF( iom_use('rnf_x_sss') ) CALL iom_put( "rnf_x_sss", rnf*pts(:,:,1,jp_sal,Kmm) ) ! runoff term on sss 192 183 193 184 #if defined key_asminc -
NEMO/trunk/src/OCE/TRA/trazdf.F90
r14433 r14834 64 64 ! 65 65 IF( kt == nit000 ) THEN 66 IF( ntile == 0.OR. ntile == 1 ) THEN ! Do only on the first tile66 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 67 67 IF(lwp)WRITE(numout,*) 68 68 IF(lwp)WRITE(numout,*) 'tra_zdf : implicit vertical mixing on T & S' -
NEMO/trunk/src/OCE/TRA/zpshde.F90
r14433 r14834 47 47 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 48 48 INTEGER , INTENT(in ) :: kjpt ! number of tracers 49 REAL(wp), DIMENSION(:,:,:,:), INTENT(in out) :: pta ! 4D tracers fields49 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pta ! 4D tracers fields 50 50 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 51 REAL(wp), DIMENSION(:,:,:) , INTENT(in out), OPTIONAL :: prd ! 3D density anomaly fields51 REAL(wp), DIMENSION(:,:,:) , INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 52 52 REAL(wp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 53 53 ! … … 111 111 INTEGER , INTENT(in ) :: kjpt ! number of tracers 112 112 INTEGER , INTENT(in ) :: ktta, ktgt, ktrd, ktgr 113 REAL(wp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(in out) :: pta ! 4D tracers fields113 REAL(wp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(in ) :: pta ! 4D tracers fields 114 114 REAL(wp), DIMENSION(A2D_T(ktgt) ,KJPT), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 115 REAL(wp), DIMENSION(A2D_T(ktrd),JPK ), INTENT(in out), OPTIONAL :: prd ! 3D density anomaly fields115 REAL(wp), DIMENSION(A2D_T(ktrd),JPK ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 116 116 REAL(wp), DIMENSION(A2D_T(ktgr) ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 117 117 ! … … 124 124 ! 125 125 IF( ln_timing ) CALL timing_start( 'zps_hde') 126 IF (nn_hls.EQ.2) THEN127 CALL lbc_lnk( 'zpshde', pta, 'T', 1.0_wp)128 IF(PRESENT(prd)) CALL lbc_lnk( 'zpshde', prd, 'T', 1.0_wp)129 END IF130 126 ! 131 127 pgtu(:,:,:) = 0._wp ; zti (:,:,:) = 0._wp ; zhi (:,:) = 0._wp … … 134 130 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! 135 131 ! 136 DO_2D( nn_hls -1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! Gradient of density at the last level132 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! Gradient of density at the last level 137 133 iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points 138 134 ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1 … … 173 169 END DO 174 170 ! 175 IF (nn_hls .EQ.1) CALL lbc_lnk( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond.171 IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond. 176 172 ! 177 173 IF( PRESENT( prd ) ) THEN !== horizontal derivative of density anomalies (rd) ==! (optional part) … … 206 202 ENDIF 207 203 END_2D 208 IF (nn_hls .EQ.1) CALL lbc_lnk( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp ) ! Lateral boundary conditions204 IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp ) ! Lateral boundary conditions 209 205 ! 210 206 END IF … … 221 217 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 222 218 INTEGER , INTENT(in ) :: kjpt ! number of tracers 223 REAL(wp), DIMENSION(:,:,:,:), INTENT(in out) :: pta ! 4D tracers fields219 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pta ! 4D tracers fields 224 220 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 225 221 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 226 REAL(wp), DIMENSION(:,:,:) , INTENT(in out), OPTIONAL :: prd ! 3D density anomaly fields222 REAL(wp), DIMENSION(:,:,:) , INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 227 223 REAL(wp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 228 224 REAL(wp), DIMENSION(:,:) , INTENT( out), OPTIONAL :: pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) … … 291 287 INTEGER , INTENT(in ) :: kjpt ! number of tracers 292 288 INTEGER , INTENT(in ) :: ktta, ktgt, ktgti, ktrd, ktgr, ktgri 293 REAL(wp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(in out) :: pta ! 4D tracers fields289 REAL(wp), DIMENSION(A2D_T(ktta),JPK,KJPT), INTENT(in ) :: pta ! 4D tracers fields 294 290 REAL(wp), DIMENSION(A2D_T(ktgt) ,KJPT), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts 295 291 REAL(wp), DIMENSION(A2D_T(ktgti) ,KJPT), INTENT( out) :: pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF) 296 REAL(wp), DIMENSION(A2D_T(ktrd),JPK ), INTENT(in out), OPTIONAL :: prd ! 3D density anomaly fields292 REAL(wp), DIMENSION(A2D_T(ktrd),JPK ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields 297 293 REAL(wp), DIMENSION(A2D_T(ktgr) ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom) 298 294 REAL(wp), DIMENSION(A2D_T(ktgri) ), INTENT( out), OPTIONAL :: pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top) … … 307 303 IF( ln_timing ) CALL timing_start( 'zps_hde_isf') 308 304 ! 309 IF (nn_hls.EQ.2) THEN310 CALL lbc_lnk( 'zpshde', pta, 'T', 1.0_wp)311 IF (PRESENT(prd)) CALL lbc_lnk( 'zpshde', prd, 'T', 1.0_wp)312 END IF313 314 305 pgtu (:,:,:) = 0._wp ; pgtv (:,:,:) =0._wp 315 306 pgtui(:,:,:) = 0._wp ; pgtvi(:,:,:) =0._wp … … 319 310 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! 320 311 ! 321 DO_2D( nn_hls -1, nn_hls-1, nn_hls-1, nn_hls-1 )312 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 322 313 323 314 iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points … … 359 350 END DO 360 351 ! 361 IF (nn_hls .EQ.1) CALL lbc_lnk( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond.352 IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgtu(:,:,:), 'U', -1.0_wp , pgtv(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond. 362 353 363 354 ! horizontal derivative of density anomalies (rd) … … 401 392 END_2D 402 393 403 IF (nn_hls .EQ.1) CALL lbc_lnk( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp ) ! Lateral boundary conditions394 IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgru , 'U', -1.0_wp , pgrv , 'V', -1.0_wp ) ! Lateral boundary conditions 404 395 ! 405 396 END IF … … 408 399 ! 409 400 DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! ! 410 DO_2D( nn_hls -1, nn_hls-1, nn_hls-1, nn_hls-1 )401 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 411 402 iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1 412 403 ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1 … … 452 443 ! 453 444 END DO 454 IF (nn_hls .EQ.1) CALL lbc_lnk( 'zpshde', pgtui(:,:,:), 'U', -1.0_wp , pgtvi(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond.445 IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgtui(:,:,:), 'U', -1.0_wp , pgtvi(:,:,:), 'V', -1.0_wp ) ! Lateral boundary cond. 455 446 456 447 IF( PRESENT( prd ) ) THEN !== horizontal derivative of density anomalies (rd) ==! (optional part) … … 491 482 492 483 END_2D 493 IF (nn_hls .EQ.1) CALL lbc_lnk( 'zpshde', pgrui, 'U', -1.0_wp , pgrvi, 'V', -1.0_wp ) ! Lateral boundary conditions484 IF (nn_hls==1) CALL lbc_lnk( 'zpshde', pgrui, 'U', -1.0_wp , pgrvi, 'V', -1.0_wp ) ! Lateral boundary conditions 494 485 ! 495 486 END IF
Note: See TracChangeset
for help on using the changeset viewer.