Changeset 6989 for branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90
- Timestamp:
- 2016-10-05T09:43:42+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limadv_umx.F90
r6860 r6989 29 29 30 30 PUBLIC lim_adv_umx ! routine called by limtrp.F90 31 32 INTEGER :: nn_ult_order = 4 ! order of the ULTIMATE scheme 33 34 REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6 31 32 REAL(wp) :: z1_6 = 1._wp / 6._wp ! =1/6 33 REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120 35 34 36 35 !! * Substitutions … … 43 42 CONTAINS 44 43 45 SUBROUTINE lim_adv_umx( lcon, kt, pdt, pu, pv, puc, pvc, pt, ptc, pfu, pfv)44 SUBROUTINE lim_adv_umx( kt, pdt, puc, pvc, pubox, pvbox, ptc ) 46 45 !!---------------------------------------------------------------------- 47 46 !! *** ROUTINE lim_adv_umx *** … … 56 55 !! ** Action : - pt the after advective tracer 57 56 !!---------------------------------------------------------------------- 58 LOGICAL , INTENT(in ) :: lcon ! "continuity" equations for a and H59 57 INTEGER , INTENT(in ) :: kt ! number of iteration 60 58 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 61 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu, pv ! 2 ice velocity components 62 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity components (for a or H) 63 ! 2 ice transport components (for tracers q) 64 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer field 59 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity components => u*e2 60 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity 65 61 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ptc ! tracer content field 66 REAL(wp), DIMENSION(jpi,jpj), INTENT( out), OPTIONAL :: pfu, pfv ! advective fluxes at u- and v-points67 62 ! 68 63 INTEGER :: ji, jj ! dummy loop indices … … 77 72 CALL wrk_alloc( jpi,jpj, zt_ups, zfu_ups, zfv_ups, ztrd, zfu_ho, zfv_ho, zt_u, zt_v ) 78 73 ! 79 !80 !clem zt_ups(:,:) = 0._wp81 74 ! 82 75 ! upstream advection with initial mass fluxes & intermediate update … … 88 81 zfp_vj = pvc(ji,jj) + ABS( pvc(ji,jj) ) 89 82 zfm_vj = pvc(ji,jj) - ABS( pvc(ji,jj) ) 90 zfu_ups(ji,jj) = 0.5_wp * ( zfp_ui * pt (ji,jj) + zfm_ui * pt(ji+1,jj ) )91 zfv_ups(ji,jj) = 0.5_wp * ( zfp_vj * pt (ji,jj) + zfm_vj * pt(ji ,jj+1) )83 zfu_ups(ji,jj) = 0.5_wp * ( zfp_ui * ptc(ji,jj) + zfm_ui * ptc(ji+1,jj ) ) 84 zfv_ups(ji,jj) = 0.5_wp * ( zfp_vj * ptc(ji,jj) + zfm_vj * ptc(ji ,jj+1) ) 92 85 END DO 93 86 END DO … … 106 99 ! High order (_ho) fluxes 107 100 ! ----------------------- 108 SELECT CASE( nn_ ult_order)101 SELECT CASE( nn_limadv_ord ) 109 102 CASE ( 20 ) ! centered second order 110 103 DO jj = 2, jpjm1 111 104 DO ji = fs_2, fs_jpim1 ! vector opt. 112 zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( pt(ji,jj) + pt(ji+1,jj) ) 113 zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( pt(ji,jj) + pt(ji,jj+1) ) 114 END DO 115 END DO 116 ! 117 CASE ( 1 , 4 ) ! 1st to 4th order ULTIMATE-MACHO scheme 118 CALL macho( lcon, kt, nn_ult_order, pdt, pt, pu, pv, zt_u, zt_v ) 119 !! CALL macho( lcon, kt, nn_ult_order, pdt, ptc, pu, pv, zt_u, zt_v ) 105 zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( ptc(ji,jj) + ptc(ji+1,jj) ) 106 zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj) + ptc(ji,jj+1) ) 107 END DO 108 END DO 109 ! 110 CASE ( 1:5 ) ! 1st to 5th order ULTIMATE-MACHO scheme 111 CALL macho( kt, nn_limadv_ord, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 120 112 ! 121 113 DO jj = 2, jpjm1 … … 153 145 CALL lbc_lnk( ptc(:,:) , 'T', 1. ) 154 146 ! 155 IF( PRESENT( pfu ) ) THEN156 DO jj = 2, jpjm1157 DO ji = fs_2, fs_jpim1 ! vector opt.158 pfu(ji,jj) = zfu_ups(ji,jj) + zfu_ho(ji,jj)159 pfv(ji,jj) = zfv_ups(ji,jj) + zfv_ho(ji,jj)160 END DO161 END DO162 CALL lbc_lnk_multi( pfu, 'U', -1., pfv, 'V', -1. ) ! Lateral bondary conditions163 ENDIF164 147 ! 165 148 CALL wrk_dealloc( jpi,jpj, zt_ups, zfu_ups, zfv_ups, ztrd, zfu_ho, zfv_ho, zt_u, zt_v ) … … 170 153 171 154 172 SUBROUTINE macho( lcon, kt, k_order, pdt, pt, pu, pv, pt_u, pt_v )155 SUBROUTINE macho( kt, k_order, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 173 156 !!--------------------------------------------------------------------- 174 157 !! *** ROUTINE ultimate_x *** … … 181 164 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 182 165 !!---------------------------------------------------------------------- 183 LOGICAL , INTENT(in ) :: lcon ! "continuity" equations for a and ah184 166 INTEGER , INTENT(in ) :: kt ! number of iteration 185 167 INTEGER , INTENT(in ) :: k_order ! order of the ULTIMATE scheme 186 168 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 187 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu, pv ! 2 ice velocity components 188 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 169 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer fields 170 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity components 171 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity 189 172 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points 190 173 ! … … 201 184 ! 202 185 ! !-- ultimate interpolation of pt at u-point --! 203 CALL ultimate_x( lcon, k_order, pdt, pt, pu, pt_u ) 204 ! 205 ! !-- advective form update in zzt --! 206 DO jj = 1, jpj 207 DO ji = fs_2, fs_jpim1 ! vector opt. 208 ! 209 IF( pu(ji,jj) * pu(ji-1,jj) <= 0._wp ) THEN ; zc_box = 0._wp 210 ELSEIF( pu(ji ,jj) > 0._wp ) THEN ; zc_box = pu(ji-1,jj) 211 ELSE ; zc_box = pu(ji ,jj) 212 ENDIF 213 ! 214 zzt(ji,jj) = pt(ji,jj) - zc_box * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e12t(ji,jj) 215 IF( lcon .eqv. .TRUE. ) THEN ! clem => for a.div(u) ?? 216 zzt(ji,jj) = zzt(ji,jj) - pt(ji,jj) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e12t(ji,jj) 217 END IF 218 END DO 219 END DO 220 ! 221 ! !-- ultimate interpolation of pt at v-point --! 222 CALL ultimate_y( lcon, k_order, pdt, zzt, pv, pt_v ) 223 ! 224 ELSE !== even ice time step: adv_y then adv_x ==! 225 ! 226 ! !-- ultimate interpolation of pt at v-point --! 227 CALL ultimate_y( lcon, k_order, pdt, pt, pv, pt_v ) 186 CALL ultimate_x( k_order, pdt, ptc, puc, pt_u ) 228 187 ! 229 188 ! !-- advective form update in zzt --! 230 189 DO jj = 2, jpjm1 231 DO ji = 1, jpi 232 ! 233 IF( pv(ji,jj) * pv(ji,jj-1) <= 0._wp ) THEN ; zc_box = 0._wp 234 ELSEIF( pv(ji,jj ) > 0._wp ) THEN ; zc_box = pv(ji,jj-1) 235 ELSE ; zc_box = pv(ji,jj ) 236 ENDIF 237 ! 238 zzt(ji,jj) = pt(ji,jj) - zc_box * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e12t(ji,jj) 239 IF( lcon .eqv. .TRUE. ) THEN ! clem => for a.div(u) ?? 240 zzt(ji,jj) = zzt(ji,jj) - pt(ji,jj) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e12t(ji,jj) 241 END IF 242 END DO 243 END DO 190 DO ji = fs_2, fs_jpim1 ! vector opt. 191 zzt(ji,jj) = ptc(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj) & 192 & - ptc (ji,jj) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e12t(ji,jj) 193 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 194 END DO 195 END DO 196 CALL lbc_lnk( zzt, 'T', 1. ) 197 ! 198 ! !-- ultimate interpolation of pt at v-point --! 199 CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v ) 200 ! 201 ELSE !== even ice time step: adv_y then adv_x ==! 202 ! 203 ! !-- ultimate interpolation of pt at v-point --! 204 CALL ultimate_y( k_order, pdt, ptc, pvc, pt_v ) 205 ! 206 ! !-- advective form update in zzt --! 207 DO jj = 2, jpjm1 208 DO ji = fs_2, fs_jpim1 209 zzt(ji,jj) = ptc(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj) & 210 & - ptc (ji,jj) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e12t(ji,jj) 211 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 212 END DO 213 END DO 214 CALL lbc_lnk( zzt, 'T', 1. ) 244 215 ! 245 216 ! !-- ultimate interpolation of pt at u-point --! 246 CALL ultimate_x( lcon, k_order, pdt, zzt, pu, pt_u )217 CALL ultimate_x( k_order, pdt, zzt, puc, pt_u ) 247 218 ! 248 219 ENDIF … … 255 226 256 227 257 SUBROUTINE ultimate_x( lcon, k_order, pdt, pt, pu, pt_u )228 SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u ) 258 229 !!--------------------------------------------------------------------- 259 230 !! *** ROUTINE ultimate_x *** … … 266 237 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 267 238 !!---------------------------------------------------------------------- 268 LOGICAL , INTENT(in ) :: lcon ! "continuity" equations for a and ah269 239 INTEGER , INTENT(in ) :: k_order ! ocean time-step index 270 240 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 271 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pu 241 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc ! ice i-velocity component 272 242 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 273 243 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u ! tracer at u-point 274 244 ! 275 245 INTEGER :: ji, jj ! dummy loop indices 276 REAL(wp) :: zcu, zdx2 ! - -277 REAL(wp), POINTER, DIMENSION(:,:) :: ztu , zltu246 REAL(wp) :: zcu, zdx2, zdx4 ! - - 247 REAL(wp), POINTER, DIMENSION(:,:) :: ztu1, ztu2, ztu3, ztu4 278 248 !!---------------------------------------------------------------------- 279 249 ! 280 250 IF( nn_timing == 1 ) CALL timing_start('ultimate_x') 281 251 ! 282 CALL wrk_alloc( jpi,jpj, ztu, zltu ) 252 CALL wrk_alloc( jpi,jpj, ztu1, ztu2, ztu3, ztu4 ) 253 ! 254 ! !-- Laplacian in i-direction --! 255 DO jj = 2, jpjm1 ! First derivative (gradient) 256 DO ji = 1, fs_jpim1 257 ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 258 END DO 259 ! ! Second derivative (Laplacian) 260 DO ji = fs_2, fs_jpim1 261 ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj) 262 END DO 263 END DO 264 CALL lbc_lnk( ztu2, 'T', 1. ) 265 ! 266 ! !-- BiLaplacian in i-direction --! 267 DO jj = 2, jpjm1 ! Third derivative 268 DO ji = 1, fs_jpim1 269 ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 270 END DO 271 ! ! Fourth derivative 272 DO ji = fs_2, fs_jpim1 273 ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj) 274 END DO 275 END DO 276 CALL lbc_lnk( ztu4, 'T', 1. ) 277 ! 283 278 ! 284 279 SELECT CASE (k_order ) … … 288 283 DO jj = 1, jpj 289 284 DO ji = 1, fs_jpim1 ! vector opt. 290 pt_u(ji,jj) = 0.5_wp * ( ( pt(ji+1,jj) + pt(ji,jj) )&291 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) * umask(ji,jj,1)285 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & 286 & - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 292 287 END DO 293 288 END DO 294 289 ! 295 290 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 291 ! 296 292 DO jj = 1, jpj 297 293 DO ji = 1, fs_jpim1 ! vector opt. 298 pt_u(ji,jj) = 0.5_wp * ( ( pt(ji+1,jj) + pt(ji,jj) ) &299 !! & - pdt * pu(ji,jj) * ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) ) * umask(ji,jj,1) 300 & - pdt * pu(ji,jj) * ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e12u(ji,jj) ) * umask(ji,jj,1)294 zcu = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 295 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & 296 & - zcu * ( pt(ji+1,jj) - pt(ji,jj) ) ) 301 297 END DO 302 298 END DO … … 305 301 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 306 302 ! 307 ! !-- Laplacian in i- and in j-directions --!308 DO jj = 2, jpjm1 ! First derivative (gradient)309 DO ji = 1, fs_jpim1 ! vector opt.310 ztu(ji,jj) = ( pt(ji+1,jj ) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)311 END DO312 ! ! Second derivative (divergence)313 DO ji = fs_2, fs_jpim1 ! vector opt.314 zltu(ji,jj) = ( ztu(ji,jj) - ztu(ji-1,jj) ) * r1_e1t(ji,jj)315 END DO316 END DO317 CALL lbc_lnk( zltu, 'T', 1. ) ! Lateral boundary cond. (unchanged sign)318 !319 303 DO jj = 1, jpj 320 304 DO ji = 1, fs_jpim1 ! vector opt. 321 !! zcu = pu(ji,jj) * pdt * r1_e1u(ji,jj) 322 zcu = pu(ji,jj) * pdt * r1_e12u(ji,jj) 305 zcu = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 323 306 zdx2 = e1u(ji,jj) * e1u(ji,jj) 307 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 324 308 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) & 325 309 & - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) & 326 & - r1_6 * zdx2 * ( 1._wp - zcu*zcu ) * ( zltu(ji+1,jj) + zltu(ji,jj) &327 & - SIGN( 1._wp, zcu ) * ( z ltu(ji+1,jj) - zltu(ji,jj) ) ) )310 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj) + ztu2(ji,jj) & 311 & - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) ) 328 312 END DO 329 313 END DO 330 314 ! 331 315 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 332 !333 ! !-- Laplacian in i- and in j-directions --!334 DO jj = 1, jpjm1 ! First derivative (gradient)335 DO ji = 1, fs_jpim1 ! vector opt.336 ztu(ji,jj) = ( pt(ji+1,jj ) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1)337 END DO338 ! ! Second derivative (divergence)339 DO ji = fs_2, fs_jpim1 ! vector opt.340 zltu(ji,jj) = ( ztu(ji,jj) - ztu(ji-1,jj) ) * r1_e1t(ji,jj)341 END DO342 END DO343 CALL lbc_lnk( zltu, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn)344 316 ! 345 317 DO jj = 1, jpj 346 318 DO ji = 1, fs_jpim1 ! vector opt. 347 !! zcu = pu(ji,jj) * pdt * r1_e1u(ji,jj) 348 zcu = pu(ji,jj) * pdt * r1_e12u(ji,jj) 319 zcu = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 349 320 zdx2 = e1u(ji,jj) * e1u(ji,jj) 321 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 350 322 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) & 351 323 & - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) & 352 & - r1_6 * zdx2 * ( 1._wp - zcu*zcu ) * ( zltu(ji+1,jj) + zltu(ji,jj) & 353 & - 0.5_wp * zcu * ( zltu(ji+1,jj) - zltu(ji,jj) ) ) ) 324 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj) + ztu2(ji,jj) & 325 & - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) ) 326 END DO 327 END DO 328 ! 329 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 330 ! 331 DO jj = 1, jpj 332 DO ji = 1, fs_jpim1 ! vector opt. 333 zcu = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 334 zdx2 = e1u(ji,jj) * e1u(ji,jj) 335 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 336 zdx4 = zdx2 * zdx2 337 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) & 338 & - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) & 339 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj) + ztu2(ji,jj) & 340 & - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) & 341 & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj) & 342 & - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) ) 354 343 END DO 355 344 END DO … … 357 346 END SELECT 358 347 ! 359 CALL wrk_dealloc( jpi,jpj, ztu , zltu)348 CALL wrk_dealloc( jpi,jpj, ztu1, ztu2, ztu3, ztu4 ) 360 349 ! 361 350 IF( nn_timing == 1 ) CALL timing_stop('ultimate_x') … … 364 353 365 354 366 SUBROUTINE ultimate_y( lcon, k_order, pdt, pt, pv, pt_v )355 SUBROUTINE ultimate_y( k_order, pdt, pt, pvc, pt_v ) 367 356 !!--------------------------------------------------------------------- 368 357 !! *** ROUTINE ultimate_y *** … … 375 364 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 376 365 !!---------------------------------------------------------------------- 377 LOGICAL , INTENT(in ) :: lcon ! "continuity" equations for a and ah378 366 INTEGER , INTENT(in ) :: k_order ! ocean time-step index 379 367 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 380 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pv 368 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvc ! ice j-velocity component 381 369 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 382 370 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_v ! tracer at v-point 383 371 ! 384 372 INTEGER :: ji, jj ! dummy loop indices 385 REAL(wp) :: zcv, zdy2 ! - -386 REAL(wp), POINTER, DIMENSION(:,:) :: ztv , zltv373 REAL(wp) :: zcv, zdy2, zdy4 ! - - 374 REAL(wp), POINTER, DIMENSION(:,:) :: ztv1, ztv2, ztv3, ztv4 387 375 !!---------------------------------------------------------------------- 388 376 ! 389 377 IF( nn_timing == 1 ) CALL timing_start('ultimate_y') 390 378 ! 391 CALL wrk_alloc( jpi,jpj, ztv, zltv ) 379 CALL wrk_alloc( jpi,jpj, ztv1, ztv2, ztv3, ztv4 ) 380 ! 381 ! !-- Laplacian in j-direction --! 382 DO jj = 1, jpjm1 ! First derivative (gradient) 383 DO ji = fs_2, fs_jpim1 384 ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 385 END DO 386 END DO 387 DO jj = 2, jpjm1 ! Second derivative (Laplacian) 388 DO ji = fs_2, fs_jpim1 389 ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj) 390 END DO 391 END DO 392 CALL lbc_lnk( ztv2, 'T', 1. ) 393 ! 394 ! !-- BiLaplacian in j-direction --! 395 DO jj = 1, jpjm1 ! First derivative 396 DO ji = fs_2, fs_jpim1 397 ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 398 END DO 399 END DO 400 DO jj = 2, jpjm1 ! Second derivative 401 DO ji = fs_2, fs_jpim1 402 ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj) 403 END DO 404 END DO 405 CALL lbc_lnk( ztv4, 'T', 1. ) 406 ! 392 407 ! 393 408 SELECT CASE (k_order ) … … 397 412 DO jj = 1, jpjm1 398 413 DO ji = 1, jpi 399 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) &400 & - SIGN( 1._wp, pv (ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ))414 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & 415 & - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 401 416 END DO 402 417 END DO … … 405 420 DO jj = 1, jpjm1 406 421 DO ji = 1, jpi 407 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) &408 !! & - pdt * pv(ji,jj) * ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) ) 409 & - pdt * pv(ji,jj) * ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e12v(ji,jj) )422 zcv = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 423 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & 424 & - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) ) 410 425 END DO 411 426 END DO … … 413 428 ! 414 429 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 415 !416 ! !-- Laplacian in i- and in j-directions --!417 DO jj = 1, jpjm1 ! First derivative (gradient)418 DO ji = fs_2, fs_jpim1 ! vector opt.419 ztv(ji,jj) = ( pt(ji ,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)420 END DO421 END DO422 DO jj = 2, jpjm1 ! Second derivative (divergence)423 DO ji = fs_2, fs_jpim1 ! vector opt.424 zltv(ji,jj) = ( ztv(ji,jj) - ztv(ji,jj-1) ) * r1_e2t(ji,jj)425 END DO426 END DO427 CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sign)428 430 ! 429 431 DO jj = 1, jpjm1 430 432 DO ji = 1, jpi 431 !! zcv = pv(ji,jj) * pdt * r1_e2v(ji,jj) 432 zcv = pv(ji,jj) * pdt * r1_e12v(ji,jj) 433 zcv = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 433 434 zdy2 = e2v(ji,jj) * e2v(ji,jj) 434 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & 435 & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & 436 & - r1_6 * zdy2 * ( 1._wp - zcv*zcv ) * ( zltv(ji,jj+1) + zltv(ji,jj) & 437 & - SIGN( 1._wp, zcv ) * ( zltv(ji,jj+1) - zltv(ji,jj) ) ) ) 435 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 436 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & 437 & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & 438 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1) + ztv2(ji,jj) & 439 & - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) ) 438 440 END DO 439 441 END DO 440 442 ! 441 443 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 442 !443 ! !-- Laplacian in i- and in j-directions --!444 DO jj = 1, jpjm1 ! First derivative (gradient)445 DO ji = fs_2, fs_jpim1 ! vector opt.446 ztv(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1)447 END DO448 END DO449 DO jj = 2, jpjm1 ! Second derivative (divergence)450 DO ji = fs_2, fs_jpim1 ! vector opt.451 zltv(ji,jj) = ( ztv(ji,jj) - ztv(ji,jj-1) ) * r1_e2t(ji,jj)452 END DO453 END DO454 CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn)455 444 ! 456 445 DO jj = 1, jpjm1 457 446 DO ji = 1, jpi 458 !! zcv = pv(ji,jj) * pdt * r1_e2v(ji,jj) 459 zcv = pv(ji,jj) * pdt * r1_e12v(ji,jj) 447 zcv = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 460 448 zdy2 = e2v(ji,jj) * e2v(ji,jj) 461 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & 462 & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & 463 & - r1_6 * zdy2 * ( 1._wp - zcv*zcv ) * ( zltv(ji,jj+1) + zltv(ji,jj) & 464 & - 0.5_wp * zcv * ( zltv(ji,jj+1) - zltv(ji,jj) ) ) ) 449 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 450 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & 451 & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & 452 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1) + ztv2(ji,jj) & 453 & - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) ) 454 END DO 455 END DO 456 ! 457 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 458 ! 459 DO jj = 1, jpjm1 460 DO ji = 1, jpi 461 zcv = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 462 zdy2 = e2v(ji,jj) * e2v(ji,jj) 463 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 464 zdy4 = zdy2 * zdy2 465 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & 466 & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & 467 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1) + ztv2(ji,jj) & 468 & - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) & 469 & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj) & 470 & - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) ) 465 471 END DO 466 472 END DO … … 468 474 END SELECT 469 475 ! 470 CALL wrk_dealloc( jpi,jpj, ztv , zltv)476 CALL wrk_dealloc( jpi,jpj, ztv1, ztv2, ztv3, ztv4 ) 471 477 ! 472 478 IF( nn_timing == 1 ) CALL timing_stop('ultimate_y')
Note: See TracChangeset
for help on using the changeset viewer.
