- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/TRA/traadv_fct.F90 (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/TRA/traadv_fct.F90
r10425 r13463 20 20 USE diaptr ! poleward transport diagnostics 21 21 USE diaar5 ! AR5 diagnostics 22 USE phycst , ONLY : rau0_rcp 22 USE phycst , ONLY : rho0_rcp 23 USE zdf_oce , ONLY : ln_zad_Aimp 23 24 ! 24 25 USE in_out_manager ! I/O manager … … 44 45 45 46 !! * Substitutions 46 # include "vectopt_loop_substitute.h90" 47 # include "do_loop_substitute.h90" 48 # include "domzgr_substitute.h90" 47 49 !!---------------------------------------------------------------------- 48 50 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 52 54 CONTAINS 53 55 54 SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, p un, pvn, pwn, &55 & ptb, ptn, pta, kjpt, kn_fct_h, kn_fct_v )56 SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pU, pV, pW, & 57 & Kbb, Kmm, pt, kjpt, Krhs, kn_fct_h, kn_fct_v ) 56 58 !!---------------------------------------------------------------------- 57 59 !! *** ROUTINE tra_adv_fct *** … … 65 67 !! - corrected flux (monotonic correction) 66 68 !! 67 !! ** Action : - update pt awith the now advective tracer trends69 !! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends 68 70 !! - send trends to trdtra module for further diagnostics (l_trdtra=T) 69 !! - htr_adv, str_adv :poleward advective heat and salt transport (ln_diaptr=T)70 !!---------------------------------------------------------------------- 71 INTEGER , INTENT(in ) :: kt ! ocean time-step index72 INTEGER , INTENT(in ) :: kit000 ! first time step index73 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)74 INTEGER , INTENT(in ) :: kjpt ! number of tracers75 INTEGER , INTENT(in ) :: kn_fct_h ! order of the FCT scheme (=2 or 4)76 INTEGER , INTENT(in ) :: kn_fct_v! order of the FCT scheme (=2 or 4)77 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step78 REAL(wp) , DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components79 REAL(wp), DIMENSION(jpi,jpj,jpk ,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields80 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt ), INTENT(inout) :: pta ! tracer trend71 !! - poleward advective heat and salt transport (ln_diaptr=T) 72 !!---------------------------------------------------------------------- 73 INTEGER , INTENT(in ) :: kt ! ocean time-step index 74 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices 75 INTEGER , INTENT(in ) :: kit000 ! first time step index 76 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 77 INTEGER , INTENT(in ) :: kjpt ! number of tracers 78 INTEGER , INTENT(in ) :: kn_fct_h ! order of the FCT scheme (=2 or 4) 79 INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4) 80 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 81 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components 82 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 81 83 ! 82 84 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 86 88 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw 87 89 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz, zptry 90 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zwinf, zwdia, zwsup 91 LOGICAL :: ll_zAimp ! flag to apply adaptive implicit vertical advection 88 92 !!---------------------------------------------------------------------- 89 93 ! … … 93 97 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 94 98 ENDIF 99 !! -- init to 0 100 zwi(:,:,:) = 0._wp 101 zwx(:,:,:) = 0._wp 102 zwy(:,:,:) = 0._wp 103 zwz(:,:,:) = 0._wp 104 ztu(:,:,:) = 0._wp 105 ztv(:,:,:) = 0._wp 106 zltu(:,:,:) = 0._wp 107 zltv(:,:,:) = 0._wp 108 ztw(:,:,:) = 0._wp 95 109 ! 96 110 l_trd = .FALSE. ! set local switches 97 111 l_hst = .FALSE. 98 112 l_ptr = .FALSE. 99 IF( ( cdtype =='TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 100 IF( cdtype =='TRA' .AND. ln_diaptr ) l_ptr = .TRUE. 101 IF( cdtype =='TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 113 ll_zAimp = .FALSE. 114 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype =='TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 115 IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE. 116 IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. & 102 117 & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE. 103 118 ! … … 117 132 zwi(:,:,:) = 0._wp 118 133 ! 134 ! If adaptive vertical advection, check if it is needed on this PE at this time 135 IF( ln_zad_Aimp ) THEN 136 IF( MAXVAL( ABS( wi(:,:,:) ) ) > 0._wp ) ll_zAimp = .TRUE. 137 END IF 138 ! If active adaptive vertical advection, build tridiagonal matrix 139 IF( ll_zAimp ) THEN 140 ALLOCATE(zwdia(jpi,jpj,jpk), zwinf(jpi,jpj,jpk),zwsup(jpi,jpj,jpk)) 141 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 142 zwdia(ji,jj,jk) = 1._wp + p2dt * ( MAX( wi(ji,jj,jk) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) & 143 & / e3t(ji,jj,jk,Krhs) 144 zwinf(ji,jj,jk) = p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) / e3t(ji,jj,jk,Krhs) 145 zwsup(ji,jj,jk) = -p2dt * MAX( wi(ji,jj,jk+1) , 0._wp ) / e3t(ji,jj,jk,Krhs) 146 END_3D 147 END IF 148 ! 119 149 DO jn = 1, kjpt !== loop over the tracers ==! 120 150 ! 121 151 ! !== upstream advection with initial mass fluxes & intermediate update ==! 122 152 ! !* upstream tracer flux in the i and j direction 123 DO jk = 1, jpkm1 124 DO jj = 1, jpjm1 125 DO ji = 1, fs_jpim1 ! vector opt. 126 ! upstream scheme 127 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 128 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 129 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 130 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 131 zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) ) 132 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) ) 133 END DO 134 END DO 135 END DO 153 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 154 ! upstream scheme 155 zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) 156 zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) ) 157 zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) ) 158 zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) ) 159 zwx(ji,jj,jk) = 0.5 * ( zfp_ui * pt(ji,jj,jk,jn,Kbb) + zfm_ui * pt(ji+1,jj ,jk,jn,Kbb) ) 160 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * pt(ji,jj,jk,jn,Kbb) + zfm_vj * pt(ji ,jj+1,jk,jn,Kbb) ) 161 END_3D 136 162 ! !* upstream tracer flux in the k direction *! 137 DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask) 138 DO jj = 1, jpj 139 DO ji = 1, jpi 140 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 141 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 142 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 143 END DO 144 END DO 145 END DO 163 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 164 zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) ) 165 zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) ) 166 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) 167 END_3D 146 168 IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked) 147 169 IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface 148 DO jj = 1, jpj 149 DO ji = 1, jpi 150 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface 151 END DO 152 END DO 170 DO_2D( 1, 1, 1, 1 ) 171 zwz(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface 172 END_2D 153 173 ELSE ! no cavities: only at the ocean surface 154 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) 174 DO_2D( 1, 1, 1, 1 ) 175 zwz(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) 176 END_2D 155 177 ENDIF 156 178 ENDIF 157 179 ! 158 DO jk = 1, jpkm1 !* trend and after field with monotonic scheme 159 DO jj = 2, jpjm1 160 DO ji = fs_2, fs_jpim1 ! vector opt. 161 ! ! total intermediate advective trends 162 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 163 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 164 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 165 ! ! update and guess with monotonic sheme 166 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 167 zwi(ji,jj,jk) = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 168 END DO 169 END DO 170 END DO 180 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 181 ! ! total intermediate advective trends 182 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 183 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 184 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 185 ! ! update and guess with monotonic sheme 186 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra & 187 & / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk) 188 zwi(ji,jj,jk) = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) & 189 & / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 190 END_3D 191 192 IF ( ll_zAimp ) THEN 193 CALL tridia_solver( zwdia, zwsup, zwinf, zwi, zwi , 0 ) 194 ! 195 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp ; 196 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 197 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 198 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 199 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) 200 zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! update vertical fluxes 201 END_3D 202 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 203 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 204 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 205 END_3D 206 ! 207 END IF 171 208 ! 172 209 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes) … … 181 218 ! 182 219 CASE( 2 ) !- 2nd order centered 183 DO jk = 1, jpkm1 184 DO jj = 1, jpjm1 185 DO ji = 1, fs_jpim1 ! vector opt. 186 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk) 187 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk) 188 END DO 189 END DO 190 END DO 220 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 221 zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj,jk,jn,Kmm) ) - zwx(ji,jj,jk) 222 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj+1,jk,jn,Kmm) ) - zwy(ji,jj,jk) 223 END_3D 191 224 ! 192 225 CASE( 4 ) !- 4th order centered … … 194 227 zltv(:,:,jpk) = 0._wp 195 228 DO jk = 1, jpkm1 ! Laplacian 196 DO jj = 1, jpjm1 ! 1st derivative (gradient) 197 DO ji = 1, fs_jpim1 ! vector opt. 198 ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 199 ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 200 END DO 201 END DO 202 DO jj = 2, jpjm1 ! 2nd derivative * 1/ 6 203 DO ji = fs_2, fs_jpim1 ! vector opt. 204 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6 205 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6 206 END DO 207 END DO 229 DO_2D( 1, 0, 1, 0 ) 230 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 231 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 232 END_2D 233 DO_2D( 0, 0, 0, 0 ) 234 zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6 235 zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6 236 END_2D 208 237 END DO 209 CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1. , zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn) 210 ! 211 DO jk = 1, jpkm1 ! Horizontal advective fluxes 212 DO jj = 1, jpjm1 213 DO ji = 1, fs_jpim1 ! vector opt. 214 zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points 215 zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) 216 ! ! C4 minus upstream advective fluxes 217 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk) 218 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk) 219 END DO 220 END DO 221 END DO 238 CALL lbc_lnk_multi( 'traadv_fct', zltu, 'T', 1.0_wp , zltv, 'T', 1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 239 ! 240 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 241 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 242 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 243 ! ! C4 minus upstream advective fluxes 244 zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk) 245 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk) 246 END_3D 222 247 ! 223 248 CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested 224 249 ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero 225 250 ztv(:,:,jpk) = 0._wp 226 DO jk = 1, jpkm1 ! 1st derivative (gradient) 227 DO jj = 1, jpjm1 228 DO ji = 1, fs_jpim1 ! vector opt. 229 ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk) 230 ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk) 231 END DO 232 END DO 233 END DO 234 CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1. , ztv, 'V', -1. ) ! Lateral boundary cond. (unchanged sgn) 235 ! 236 DO jk = 1, jpkm1 ! Horizontal advective fluxes 237 DO jj = 2, jpjm1 238 DO ji = 2, fs_jpim1 ! vector opt. 239 zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points (x2) 240 zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) 241 ! ! C4 interpolation of T at u- & v-points (x2) 242 zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) ) 243 zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) ) 244 ! ! C4 minus upstream advective fluxes 245 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 246 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 247 END DO 248 END DO 249 END DO 251 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 252 ztu(ji,jj,jk) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk) 253 ztv(ji,jj,jk) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk) 254 END_3D 255 CALL lbc_lnk_multi( 'traadv_fct', ztu, 'U', -1.0_wp , ztv, 'V', -1.0_wp ) ! Lateral boundary cond. (unchanged sgn) 256 ! 257 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 258 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) 259 zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) 260 ! ! C4 interpolation of T at u- & v-points (x2) 261 zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) ) 262 zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) ) 263 ! ! C4 minus upstream advective fluxes 264 zwx(ji,jj,jk) = 0.5_wp * pU(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk) 265 zwy(ji,jj,jk) = 0.5_wp * pV(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk) 266 END_3D 250 267 ! 251 268 END SELECT … … 254 271 ! 255 272 CASE( 2 ) !- 2nd order centered 256 DO jk = 2, jpkm1 257 DO jj = 2, jpjm1 258 DO ji = fs_2, fs_jpim1 259 zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) & 260 & - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 261 END DO 262 END DO 263 END DO 273 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 274 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * 0.5_wp * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) & 275 & - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 276 END_3D 264 277 ! 265 278 CASE( 4 ) !- 4th order COMPACT 266 CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! zwt = COMPACT interpolation of T at w-point 267 DO jk = 2, jpkm1 268 DO jj = 2, jpjm1 269 DO ji = fs_2, fs_jpim1 270 zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 271 END DO 272 END DO 273 END DO 279 CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! zwt = COMPACT interpolation of T at w-point 280 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 281 zwz(ji,jj,jk) = ( pW(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk) 282 END_3D 274 283 ! 275 284 END SELECT … … 277 286 zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked 278 287 ENDIF 279 ! 280 CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1., zwx, 'U', -1. , zwy, 'V', -1., zwz, 'W', 1. ) 288 ! 289 IF ( ll_zAimp ) THEN 290 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 291 ! ! total intermediate advective trends 292 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 293 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 294 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 295 ztw(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 296 END_3D 297 ! 298 CALL tridia_solver( zwdia, zwsup, zwinf, ztw, ztw , 0 ) 299 ! 300 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 301 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 302 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 303 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) 304 END_3D 305 END IF 306 ! 307 CALL lbc_lnk_multi( 'traadv_fct', zwi, 'T', 1.0_wp, zwx, 'U', -1.0_wp , zwy, 'V', -1.0_wp, zwz, 'W', 1.0_wp ) 281 308 ! 282 309 ! !== monotonicity algorithm ==! 283 310 ! 284 CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )311 CALL nonosc( Kmm, pt(:,:,:,jn,Kbb), zwx, zwy, zwz, zwi, p2dt ) 285 312 ! 286 313 ! !== final trend with corrected fluxes ==! 287 314 ! 288 DO jk = 1, jpkm1 289 DO jj = 2, jpjm1 290 DO ji = fs_2, fs_jpim1 ! vector opt. 291 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 292 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 293 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) & 294 & * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 295 END DO 296 END DO 297 END DO 315 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 316 ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 317 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 318 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj) 319 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztra / e3t(ji,jj,jk,Kmm) 320 zwi(ji,jj,jk) = zwi(ji,jj,jk) + p2dt * ztra / e3t(ji,jj,jk,Krhs) * tmask(ji,jj,jk) 321 END_3D 322 ! 323 IF ( ll_zAimp ) THEN 324 ! 325 ztw(:,:,1) = 0._wp ; ztw(:,:,jpk) = 0._wp 326 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 327 zfp_wk = wi(ji,jj,jk) + ABS( wi(ji,jj,jk) ) 328 zfm_wk = wi(ji,jj,jk) - ABS( wi(ji,jj,jk) ) 329 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) 330 zwz(ji,jj,jk) = zwz(ji,jj,jk) + ztw(ji,jj,jk) ! Update vertical fluxes for trend diagnostic 331 END_3D 332 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 333 pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) & 334 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 335 END_3D 336 END IF 298 337 ! 299 338 IF( l_trd .OR. l_hst ) THEN ! trend diagnostics // heat/salt transport … … 303 342 ! 304 343 IF( l_trd ) THEN ! trend diagnostics 305 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )306 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )307 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )344 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztrdx, pU, pt(:,:,:,jn,Kmm) ) 345 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztrdy, pV, pt(:,:,:,jn,Kmm) ) 346 CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, ztrdz, pW, pt(:,:,:,jn,Kmm) ) 308 347 ENDIF 309 348 ! ! heat/salt transport … … 318 357 END DO ! end of tracer loop 319 358 ! 359 IF ( ll_zAimp ) THEN 360 DEALLOCATE( zwdia, zwinf, zwsup ) 361 ENDIF 320 362 IF( l_trd .OR. l_hst ) THEN 321 363 DEALLOCATE( ztrdx, ztrdy, ztrdz ) … … 328 370 329 371 330 SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )372 SUBROUTINE nonosc( Kmm, pbef, paa, pbb, pcc, paft, p2dt ) 331 373 !!--------------------------------------------------------------------- 332 374 !! *** ROUTINE nonosc *** … … 341 383 !! in-space based differencing for fluid 342 384 !!---------------------------------------------------------------------- 385 INTEGER , INTENT(in ) :: Kmm ! time level index 343 386 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 344 387 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field … … 347 390 INTEGER :: ji, jj, jk ! dummy loop indices 348 391 INTEGER :: ikm1 ! local integer 349 REAL( wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars350 REAL( wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -351 REAL( wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo352 !!---------------------------------------------------------------------- 353 ! 354 zbig = 1.e+40_ wp355 zrtrn = 1.e-15_ wp356 zbetup(:,:,:) = 0._ wp ; zbetdo(:,:,:) = 0._wp392 REAL(dp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars 393 REAL(dp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 394 REAL(dp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo 395 !!---------------------------------------------------------------------- 396 ! 397 zbig = 1.e+40_dp 398 zrtrn = 1.e-15_dp 399 zbetup(:,:,:) = 0._dp ; zbetdo(:,:,:) = 0._dp 357 400 358 401 ! Search local extrema … … 366 409 DO jk = 1, jpkm1 367 410 ikm1 = MAX(jk-1,1) 368 DO jj = 2, jpjm1 369 DO ji = fs_2, fs_jpim1 ! vector opt. 370 371 ! search maximum in neighbourhood 372 zup = MAX( zbup(ji ,jj ,jk ), & 373 & zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), & 374 & zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), & 375 & zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) ) 376 377 ! search minimum in neighbourhood 378 zdo = MIN( zbdo(ji ,jj ,jk ), & 379 & zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), & 380 & zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), & 381 & zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) ) 382 383 ! positive part of the flux 384 zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) & 385 & + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) & 386 & + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) ) 387 388 ! negative part of the flux 389 zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) & 390 & + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) & 391 & + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) ) 392 393 ! up & down beta terms 394 zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt 395 zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 396 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt 397 END DO 398 END DO 411 DO_2D( 0, 0, 0, 0 ) 412 413 ! search maximum in neighbourhood 414 zup = MAX( zbup(ji ,jj ,jk ), & 415 & zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), & 416 & zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), & 417 & zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) ) 418 419 ! search minimum in neighbourhood 420 zdo = MIN( zbdo(ji ,jj ,jk ), & 421 & zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), & 422 & zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), & 423 & zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) ) 424 425 ! positive part of the flux 426 zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) & 427 & + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) & 428 & + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) ) 429 430 ! negative part of the flux 431 zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) & 432 & + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) & 433 & + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) ) 434 435 ! up & down beta terms 436 zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt 437 zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 438 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt 439 END_2D 399 440 END DO 400 CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1. , zbetdo, 'T', 1.) ! lateral boundary cond. (unchanged sign)441 CALL lbc_lnk_multi( 'traadv_fct', zbetup, 'T', 1.0_wp , zbetdo, 'T', 1.0_wp ) ! lateral boundary cond. (unchanged sign) 401 442 402 443 ! 3. monotonic flux in the i & j direction (paa & pbb) 403 444 ! ---------------------------------------- 404 DO jk = 1, jpkm1 405 DO jj = 2, jpjm1 406 DO ji = fs_2, fs_jpim1 ! vector opt. 407 zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 408 zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 409 zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) ) 410 paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 411 412 zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 413 zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 414 zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) ) 415 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 416 417 ! monotonic flux in the k direction, i.e. pcc 418 ! ------------------------------------------- 419 za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 420 zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 421 zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) ) 422 pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 423 END DO 424 END DO 425 END DO 426 CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1. , pbb, 'V', -1. ) ! lateral boundary condition (changed sign) 445 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 446 zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 447 zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 448 zcu = ( 0.5 + SIGN( 0.5_wp , paa(ji,jj,jk) ) ) 449 paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 450 451 zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 452 zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 453 zcv = ( 0.5 + SIGN( 0.5_wp , pbb(ji,jj,jk) ) ) 454 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 455 456 ! monotonic flux in the k direction, i.e. pcc 457 ! ------------------------------------------- 458 za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) ) 459 zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 460 zc = ( 0.5 + SIGN( 0.5_wp , pcc(ji,jj,jk+1) ) ) 461 pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 462 END_3D 463 CALL lbc_lnk_multi( 'traadv_fct', paa, 'U', -1.0_wp , pbb, 'V', -1.0_wp ) ! lateral boundary condition (changed sign) 427 464 ! 428 465 END SUBROUTINE nonosc … … 444 481 !!---------------------------------------------------------------------- 445 482 446 DO jk = 3, jpkm1 !== build the three diagonal matrix ==! 447 DO jj = 1, jpj 448 DO ji = 1, jpi 449 zwd (ji,jj,jk) = 4._wp 450 zwi (ji,jj,jk) = 1._wp 451 zws (ji,jj,jk) = 1._wp 452 zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 453 ! 454 IF( tmask(ji,jj,jk+1) == 0._wp) THEN ! Switch to second order centered at bottom 455 zwd (ji,jj,jk) = 1._wp 456 zwi (ji,jj,jk) = 0._wp 457 zws (ji,jj,jk) = 0._wp 458 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 459 ENDIF 460 END DO 461 END DO 462 END DO 463 ! 464 jk = 2 ! Switch to second order centered at top 465 DO jj = 1, jpj 466 DO ji = 1, jpi 483 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 484 zwd (ji,jj,jk) = 4._wp 485 zwi (ji,jj,jk) = 1._wp 486 zws (ji,jj,jk) = 1._wp 487 zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 488 ! 489 IF( tmask(ji,jj,jk+1) == 0._wp) THEN ! Switch to second order centered at bottom 467 490 zwd (ji,jj,jk) = 1._wp 468 491 zwi (ji,jj,jk) = 0._wp 469 492 zws (ji,jj,jk) = 0._wp 470 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 471 END DO 472 END DO 493 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 494 ENDIF 495 END_3D 496 ! 497 jk = 2 ! Switch to second order centered at top 498 DO_2D( 1, 1, 1, 1 ) 499 zwd (ji,jj,jk) = 1._wp 500 zwi (ji,jj,jk) = 0._wp 501 zws (ji,jj,jk) = 0._wp 502 zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) ) 503 END_2D 473 504 ! 474 505 ! !== tridiagonal solve ==! 475 DO jj = 1, jpj ! first recurrence 476 DO ji = 1, jpi 477 zwt(ji,jj,2) = zwd(ji,jj,2) 478 END DO 479 END DO 480 DO jk = 3, jpkm1 481 DO jj = 1, jpj 482 DO ji = 1, jpi 483 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 484 END DO 485 END DO 486 END DO 487 ! 488 DO jj = 1, jpj ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 489 DO ji = 1, jpi 490 pt_out(ji,jj,2) = zwrm(ji,jj,2) 491 END DO 492 END DO 493 DO jk = 3, jpkm1 494 DO jj = 1, jpj 495 DO ji = 1, jpi 496 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 497 END DO 498 END DO 499 END DO 500 501 DO jj = 1, jpj ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 502 DO ji = 1, jpi 503 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 504 END DO 505 END DO 506 DO jk = jpk-2, 2, -1 507 DO jj = 1, jpj 508 DO ji = 1, jpi 509 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 510 END DO 511 END DO 512 END DO 506 DO_2D( 1, 1, 1, 1 ) 507 zwt(ji,jj,2) = zwd(ji,jj,2) 508 END_2D 509 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 510 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 511 END_3D 512 ! 513 DO_2D( 1, 1, 1, 1 ) 514 pt_out(ji,jj,2) = zwrm(ji,jj,2) 515 END_2D 516 DO_3D( 1, 1, 1, 1, 3, jpkm1 ) 517 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 518 END_3D 519 520 DO_2D( 1, 1, 1, 1 ) 521 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 522 END_2D 523 DO_3DS( 1, 1, 1, 1, jpk-2, 2, -1 ) 524 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 525 END_3D 513 526 ! 514 527 END SUBROUTINE interp_4th_cpt_org … … 533 546 ! !== build the three diagonal matrix & the RHS ==! 534 547 ! 535 DO jk = 3, jpkm1 ! interior (from jk=3 to jpk-1) 536 DO jj = 2, jpjm1 537 DO ji = fs_2, fs_jpim1 538 zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal 539 zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal 540 zws (ji,jj,jk) = wmask(ji,jj,jk) ! upper diagonal 541 zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk) & ! RHS 542 & * ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) ) 543 END DO 544 END DO 545 END DO 548 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 549 zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal 550 zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal 551 zws (ji,jj,jk) = wmask(ji,jj,jk) ! upper diagonal 552 zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk) & ! RHS 553 & * ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) ) 554 END_3D 546 555 ! 547 556 !!gm … … 556 565 END IF 557 566 ! 558 DO jj = 2, jpjm1 ! 2nd order centered at top & bottom 559 DO ji = fs_2, fs_jpim1 560 ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point 561 ikb = mbkt(ji,jj) ! - above the last wet point 562 ! 563 zwd (ji,jj,ikt) = 1._wp ! top 564 zwi (ji,jj,ikt) = 0._wp 565 zws (ji,jj,ikt) = 0._wp 566 zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) ) 567 ! 568 zwd (ji,jj,ikb) = 1._wp ! bottom 569 zwi (ji,jj,ikb) = 0._wp 570 zws (ji,jj,ikb) = 0._wp 571 zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 572 END DO 573 END DO 567 DO_2D( 0, 0, 0, 0 ) 568 ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point 569 ikb = MAX(mbkt(ji,jj), 2) ! - above the last wet point 570 ! 571 zwd (ji,jj,ikt) = 1._wp ! top 572 zwi (ji,jj,ikt) = 0._wp 573 zws (ji,jj,ikt) = 0._wp 574 zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,ikt-1) + pt_in(ji,jj,ikt) ) 575 ! 576 zwd (ji,jj,ikb) = 1._wp ! bottom 577 zwi (ji,jj,ikb) = 0._wp 578 zws (ji,jj,ikb) = 0._wp 579 zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,ikb-1) + pt_in(ji,jj,ikb) ) 580 END_2D 574 581 ! 575 582 ! !== tridiagonal solver ==! 576 583 ! 577 DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 578 DO ji = fs_2, fs_jpim1 579 zwt(ji,jj,2) = zwd(ji,jj,2) 580 END DO 581 END DO 582 DO jk = 3, jpkm1 583 DO jj = 2, jpjm1 584 DO ji = fs_2, fs_jpim1 585 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 586 END DO 587 END DO 588 END DO 589 ! 590 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 591 DO ji = fs_2, fs_jpim1 592 pt_out(ji,jj,2) = zwrm(ji,jj,2) 593 END DO 594 END DO 595 DO jk = 3, jpkm1 596 DO jj = 2, jpjm1 597 DO ji = fs_2, fs_jpim1 598 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 599 END DO 600 END DO 601 END DO 602 603 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 604 DO ji = fs_2, fs_jpim1 605 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 606 END DO 607 END DO 608 DO jk = jpk-2, 2, -1 609 DO jj = 2, jpjm1 610 DO ji = fs_2, fs_jpim1 611 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 612 END DO 613 END DO 614 END DO 584 DO_2D( 0, 0, 0, 0 ) 585 zwt(ji,jj,2) = zwd(ji,jj,2) 586 END_2D 587 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 588 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 589 END_3D 590 ! 591 DO_2D( 0, 0, 0, 0 ) 592 pt_out(ji,jj,2) = zwrm(ji,jj,2) 593 END_2D 594 DO_3D( 0, 0, 0, 0, 3, jpkm1 ) 595 pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 596 END_3D 597 598 DO_2D( 0, 0, 0, 0 ) 599 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 600 END_2D 601 DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 ) 602 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 603 END_3D 615 604 ! 616 605 END SUBROUTINE interp_4th_cpt … … 649 638 kstart = 1 + klev 650 639 ! 651 DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 652 DO ji = fs_2, fs_jpim1 653 zwt(ji,jj,kstart) = pD(ji,jj,kstart) 654 END DO 655 END DO 656 DO jk = kstart+1, jpkm1 657 DO jj = 2, jpjm1 658 DO ji = fs_2, fs_jpim1 659 zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 660 END DO 661 END DO 662 END DO 663 ! 664 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 665 DO ji = fs_2, fs_jpim1 666 pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 667 END DO 668 END DO 669 DO jk = kstart+1, jpkm1 670 DO jj = 2, jpjm1 671 DO ji = fs_2, fs_jpim1 672 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 673 END DO 674 END DO 675 END DO 676 677 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 678 DO ji = fs_2, fs_jpim1 679 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 680 END DO 681 END DO 682 DO jk = jpk-2, kstart, -1 683 DO jj = 2, jpjm1 684 DO ji = fs_2, fs_jpim1 685 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 686 END DO 687 END DO 688 END DO 640 DO_2D( 0, 0, 0, 0 ) 641 zwt(ji,jj,kstart) = pD(ji,jj,kstart) 642 END_2D 643 DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 ) 644 zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1) 645 END_3D 646 ! 647 DO_2D( 0, 0, 0, 0 ) 648 pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart) 649 END_2D 650 DO_3D( 0, 0, 0, 0, kstart+1, jpkm1 ) 651 pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1) 652 END_3D 653 654 DO_2D( 0, 0, 0, 0 ) 655 pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) 656 END_2D 657 DO_3DS( 0, 0, 0, 0, jpk-2, kstart, -1 ) 658 pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk) 659 END_3D 689 660 ! 690 661 END SUBROUTINE tridia_solver
Note: See TracChangeset
for help on using the changeset viewer.
