Changeset 10180
- Timestamp:
- 2018-10-08T17:06:04+02:00 (6 years ago)
- Location:
- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src
- Files:
-
- 4 edited
-
ICE/icedyn_adv_umx.F90 (modified) (29 diffs)
-
ICE/icethd.F90 (modified) (2 diffs)
-
ICE/icethd_zdf_bl99.F90 (modified) (13 diffs)
-
OCE/LBC/lib_mpp.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icedyn_adv_umx.F90
r10170 r10180 71 71 INTEGER :: ji, jj, jk, jl, jt ! dummy loop indices 72 72 INTEGER :: initad ! number of sub-timestep for the advection 73 INTEGER :: ipl ! third dimention of tracer array 74 73 75 REAL(wp) :: zcfl , zusnit, zdt ! - - 74 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zudy, zvdx, zcu_box, zcv_box 76 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zudy, zvdx, zcu_box, zcv_box 77 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zpato 75 78 !!---------------------------------------------------------------------- 76 79 ! … … 78 81 ! 79 82 ALLOCATE( zudy(jpi,jpj) , zvdx(jpi,jpj) , zcu_box(jpi,jpj) , zcv_box(jpi,jpj) ) 83 ALLOCATE( zpato(jpi,jpj,1) ) 80 84 ! 81 85 ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- ! … … 109 113 END DO 110 114 115 zpato (:,:,1) = pato_i(:,:) 116 111 117 !---------------! 112 118 !== advection ==! 113 119 !---------------! 114 120 DO jt = 1, initad 115 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pato_i(:,:) )! Open water area116 DO jl = 1, jpl117 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,jl) ) ! Ice area118 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,jl) ) ! Ice volume119 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,jl) ) ! Saltcontent120 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,jl) ) ! Age content121 DO jk = 1, nlay_i122 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,jl) ) ! Ice heat content123 END DO124 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) ) ! Snow volume125 DO jk = 1, nlay_s126 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,jl) ) ! Snow heat content127 END DO128 IF ( ln_pnd_H12 ) THEN129 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) ) ! Melt pond fraction130 CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) ) ! Melt pond volume131 ENDIF132 END DO133 END DO134 ! 135 DEALLOCATE( zudy, zvdx, zcu_box, zcv_box )121 CALL adv_umx( k_order, kt, 1, zdt, zudy, zvdx, zcu_box, zcv_box, zpato(:,:,1) ) ! Open water area 122 CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pa_i(:,:,:) ) ! Ice area 123 CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_i(:,:,:) ) ! Ice volume 124 CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, psv_i(:,:,:) ) ! Salt content 125 CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, poa_i(:,:,:) ) ! Age content 126 DO jk = 1, nlay_i 127 CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pe_i(:,:,jk,:) ) ! Ice heat content 128 END DO 129 CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,:) ) ! Snow volume 130 DO jk = 1, nlay_s 131 CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,jk,:) ) ! Snow heat content 132 END DO 133 IF ( ln_pnd_H12 ) THEN 134 CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,:) ) ! Melt pond fraction 135 CALL adv_umx( k_order, kt, jpl, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,:) ) ! Melt pond volume 136 ENDIF 137 END DO 138 ! 139 pato_i(:,:) = zpato (:,:,1) 140 ! 141 DEALLOCATE( zudy, zvdx, zcu_box, zcv_box, zpato ) 136 142 ! 137 143 END SUBROUTINE ice_dyn_adv_umx 138 144 139 145 140 SUBROUTINE adv_umx( k_order, kt, pdt, puc, pvc, pubox, pvbox, ptc )146 SUBROUTINE adv_umx( k_order, kt, ipl, pdt, puc, pvc, pubox, pvbox, ptc ) 141 147 !!---------------------------------------------------------------------- 142 148 !! *** ROUTINE adv_umx *** … … 151 157 !! ** Action : - pt the after advective tracer 152 158 !!---------------------------------------------------------------------- 153 INTEGER , INTENT(in ) :: k_order ! order of the ULTIMATE scheme 154 INTEGER , INTENT(in ) :: kt ! number of iteration 155 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 156 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 157 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity 158 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: ptc ! tracer content field 159 ! 160 INTEGER :: ji, jj ! dummy loop indices 159 INTEGER , INTENT(in ) :: k_order ! order of the ULTIMATE scheme 160 INTEGER , INTENT(in ) :: kt ! number of iteration 161 INTEGER , INTENT(in ) :: ipl ! third dimension of tracer array 162 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 163 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: puc , pvc ! 2 ice velocity components => u*e2 164 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pubox, pvbox ! upstream velocity 165 REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(inout) :: ptc ! tracer content field 166 ! 167 INTEGER :: ji, jj, jl ! dummy loop indices 168 161 169 REAL(wp) :: ztra ! local scalar 162 REAL(wp), DIMENSION(jpi,jpj) :: zfu_ups, zfu_ho, zt_u, zt_ups 163 REAL(wp), DIMENSION(jpi,jpj) :: zfv_ups, zfv_ho, zt_v, ztrd 170 REAL(wp), DIMENSION(jpi,jpj,ipl) :: zfu_ups, zfu_ho, zt_u, zt_ups 171 REAL(wp), DIMENSION(jpi,jpj,ipl) :: zfv_ups, zfv_ho, zt_v, ztrd 172 173 DO jl = 1, ipl 164 174 !!---------------------------------------------------------------------- 165 175 ! … … 168 178 DO jj = 1, jpjm1 ! upstream tracer flux in the i and j direction 169 179 DO ji = 1, fs_jpim1 ! vector opt. 170 zfu_ups(ji,jj ) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj)171 zfv_ups(ji,jj ) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1)180 zfu_ups(ji,jj,jl) = MAX( puc(ji,jj), 0._wp ) * ptc(ji,jj,jl) + MIN( puc(ji,jj), 0._wp ) * ptc(ji+1,jj,jl) 181 zfv_ups(ji,jj,jl) = MAX( pvc(ji,jj), 0._wp ) * ptc(ji,jj,jl) + MIN( pvc(ji,jj), 0._wp ) * ptc(ji,jj+1,jl) 172 182 END DO 173 183 END DO 174 184 175 DO jj = 2, jpjm1 ! total intermediate advective trends176 DO ji = fs_2, fs_jpim1 ! vector opt.177 ztra = - ( zfu_ups(ji,jj) - zfu_ups(ji-1,jj) &178 & + zfv_ups(ji,jj) - zfv_ups(ji ,jj-1) ) * r1_e1e2t(ji,jj)185 DO jj = 2, jpjm1 ! total intermediate advective trends 186 DO ji = fs_2, fs_jpim1 ! vector opt. 187 ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj ,jl) & 188 & + zfv_ups(ji,jj,jl) - zfv_ups(ji ,jj-1,jl) ) * r1_e1e2t(ji,jj) 179 189 ! 180 ztrd(ji,jj) = ztra ! upstream trend [ -div(uh) or -div(uhT) ]181 zt_ups (ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1)! guess after content field with monotonic scheme182 END DO183 END DO184 CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. ) ! Lateral boundary conditions (unchanged sign)190 ztrd(ji,jj,jl) = ztra ! upstream trend [ -div(uh) or -div(uhT) ] 191 zt_ups (ji,jj,jl) = ( ptc(ji,jj,jl) + pdt * ztra ) * tmask(ji,jj,1) ! guess after content field with monotonic scheme 192 END DO 193 END DO 194 END DO 185 195 186 196 ! High order (_ho) fluxes … … 188 198 SELECT CASE( k_order ) 189 199 CASE ( 20 ) ! centered second order 190 DO jj = 1, jpjm1 191 DO ji = 1, fs_jpim1 ! vector opt. 192 zfu_ho(ji,jj) = 0.5 * puc(ji,jj) * ( ptc(ji,jj) + ptc(ji+1,jj) ) 193 zfv_ho(ji,jj) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj) + ptc(ji,jj+1) ) 200 DO jl = 1, ipl 201 DO jj = 1, jpjm1 202 DO ji = 1, fs_jpim1 ! vector opt. 203 zfu_ho(ji,jj,jl) = 0.5 * puc(ji,jj) * ( ptc(ji,jj,jl) + ptc(ji+1,jj,jl) ) 204 zfv_ho(ji,jj,jl) = 0.5 * pvc(ji,jj) * ( ptc(ji,jj,jl) + ptc(ji,jj+1,jl) ) 205 END DO 194 206 END DO 195 207 END DO 196 208 ! 197 209 CASE ( 1:5 ) ! 1st to 5th order ULTIMATE-MACHO scheme 198 CALL macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 199 ! 200 DO jj = 1, jpjm1 201 DO ji = 1, fs_jpim1 ! vector opt. 202 zfu_ho(ji,jj) = puc(ji,jj) * zt_u(ji,jj) 203 zfv_ho(ji,jj) = pvc(ji,jj) * zt_v(ji,jj) 210 CALL macho( k_order, kt, ipl, pdt, ptc, puc, pvc, pubox, pvbox, zt_u, zt_v ) 211 ! 212 DO jl = 1, ipl 213 DO jj = 1, jpjm1 214 DO ji = 1, fs_jpim1 ! vector opt. 215 zfu_ho(ji,jj,jl) = puc(ji,jj) * zt_u(ji,jj,jl) 216 zfv_ho(ji,jj,jl) = pvc(ji,jj) * zt_v(ji,jj,jl) 217 END DO 204 218 END DO 205 219 END DO … … 209 223 ! antidiffusive flux : high order minus low order 210 224 ! -------------------------------------------------- 211 DO jj = 1, jpjm1 212 DO ji = 1, fs_jpim1 ! vector opt. 213 zfu_ho(ji,jj) = zfu_ho(ji,jj) - zfu_ups(ji,jj) 214 zfv_ho(ji,jj) = zfv_ho(ji,jj) - zfv_ups(ji,jj) 215 END DO 216 END DO 225 DO jl = 1, ipl 226 DO jj = 1, jpjm1 227 DO ji = 1, fs_jpim1 ! vector opt. 228 zfu_ho(ji,jj,jl) = zfu_ho(ji,jj,jl) - zfu_ups(ji,jj,jl) 229 zfv_ho(ji,jj,jl) = zfv_ho(ji,jj,jl) - zfv_ups(ji,jj,jl) 230 END DO 231 END DO 232 END DO 233 234 CALL lbc_lnk("icedyn_adv_umx",zt_ups, 'T', 1. ) ! Lateral boundary conditions (unchanged sign) 217 235 218 236 ! monotonicity algorithm 219 237 ! ------------------------- 220 CALL nonosc_2d( ptc, zfu_ho, zfv_ho, zt_ups, pdt )238 CALL nonosc_2d( ipl, ptc, zfu_ho, zfv_ho, zt_ups, pdt ) 221 239 222 240 ! final trend with corrected fluxes 223 241 ! ------------------------------------ 224 DO jj = 2, jpjm1 225 DO ji = fs_2, fs_jpim1 ! vector opt. 226 ztra = ztrd(ji,jj) - ( zfu_ho(ji,jj) - zfu_ho(ji-1,jj ) & 227 & + zfv_ho(ji,jj) - zfv_ho(ji ,jj-1) ) * r1_e1e2t(ji,jj) 228 ptc(ji,jj) = ( ptc(ji,jj) + pdt * ztra ) * tmask(ji,jj,1) 242 DO jl = 1, ipl 243 DO jj = 2, jpjm1 244 DO ji = fs_2, fs_jpim1 ! vector opt. 245 ztra = ztrd(ji,jj,jl) - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj ,jl) & 246 & + zfv_ho(ji,jj,jl) - zfv_ho(ji ,jj-1,jl) ) * r1_e1e2t(ji,jj) 247 ptc(ji,jj,jl) = ptc(ji,jj,jl) + pdt * ztra * tmask(ji,jj,1) 248 END DO 229 249 END DO 230 250 END DO … … 233 253 END SUBROUTINE adv_umx 234 254 235 236 SUBROUTINE macho( k_order, kt, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 255 SUBROUTINE macho( k_order, kt, ipl, pdt, ptc, puc, pvc, pubox, pvbox, pt_u, pt_v ) 237 256 !!--------------------------------------------------------------------- 238 257 !! *** ROUTINE ultimate_x *** … … 245 264 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 246 265 !!---------------------------------------------------------------------- 247 INTEGER , INTENT(in ) :: k_order ! order of the ULTIMATE scheme 248 INTEGER , INTENT(in ) :: kt ! number of iteration 249 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 250 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ptc ! tracer fields 251 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc, pvc ! 2 ice velocity components 252 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pubox, pvbox ! upstream velocity 253 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points 254 ! 255 INTEGER :: ji, jj ! dummy loop indices 266 INTEGER , INTENT(in ) :: k_order ! order of the ULTIMATE scheme 267 INTEGER , INTENT(in ) :: kt ! number of iteration 268 INTEGER , INTENT(in ) :: ipl ! third dimension of tracer array 269 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 270 REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in ) :: ptc ! tracer fields 271 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: puc, pvc ! 2 ice velocity components 272 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pubox, pvbox ! upstream velocity 273 REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points 274 ! 275 INTEGER :: ji, jj, jl ! dummy loop indices 256 276 REAL(wp) :: zc_box ! - - 257 REAL(wp), DIMENSION(jpi,jpj ) :: zzt277 REAL(wp), DIMENSION(jpi,jpj,ipl) :: zzt 258 278 !!---------------------------------------------------------------------- 259 279 ! … … 261 281 ! 262 282 ! !-- ultimate interpolation of pt at u-point --! 263 CALL ultimate_x( k_order, pdt, ptc, puc, pt_u )283 CALL ultimate_x( k_order, ipl, pdt, ptc, puc, pt_u ) 264 284 ! 265 285 ! !-- advective form update in zzt --! 266 DO jj = 2, jpjm1 267 DO ji = fs_2, fs_jpim1 ! vector opt. 268 zzt(ji,jj) = ptc(ji,jj) - pubox(ji,jj) * pdt * ( pt_u(ji,jj) - pt_u(ji-1,jj) ) * r1_e1t(ji,jj) & 269 & - ptc (ji,jj) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj) 270 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 286 DO jl = 1, ipl 287 DO jj = 2, jpjm1 288 DO ji = fs_2, fs_jpim1 ! vector opt. 289 zzt(ji,jj,jl) = ptc(ji,jj,jl) - pubox(ji,jj) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj) & 290 & - ptc(ji,jj,jl) * pdt * ( puc (ji,jj) - puc (ji-1,jj) ) * r1_e1e2t(ji,jj) 291 zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 292 END DO 271 293 END DO 272 294 END DO … … 274 296 ! 275 297 ! !-- ultimate interpolation of pt at v-point --! 276 CALL ultimate_y( k_order, pdt, zzt, pvc, pt_v )298 CALL ultimate_y( k_order, ipl, pdt, zzt, pvc, pt_v ) 277 299 ! 278 300 ELSE !== even ice time step: adv_y then adv_x ==! 279 301 ! 280 302 ! !-- ultimate interpolation of pt at v-point --! 281 CALL ultimate_y( k_order, pdt, ptc, pvc, pt_v )303 CALL ultimate_y( k_order, ipl, pdt, ptc, pvc, pt_v ) 282 304 ! 283 305 ! !-- advective form update in zzt --! 284 DO jj = 2, jpjm1 285 DO ji = fs_2, fs_jpim1 286 zzt(ji,jj) = ptc(ji,jj) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj) - pt_v(ji,jj-1) ) * r1_e2t(ji,jj) & 287 & - ptc (ji,jj) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj) 288 zzt(ji,jj) = zzt(ji,jj) * tmask(ji,jj,1) 306 DO jl = 1, ipl 307 DO jj = 2, jpjm1 308 DO ji = fs_2, fs_jpim1 309 zzt(ji,jj,jl) = ptc(ji,jj,jl) - pvbox(ji,jj) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj) & 310 & - ptc (ji,jj,jl) * pdt * ( pvc (ji,jj) - pvc (ji,jj-1) ) * r1_e1e2t(ji,jj) 311 zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 312 END DO 289 313 END DO 290 314 END DO … … 292 316 ! 293 317 ! !-- ultimate interpolation of pt at u-point --! 294 CALL ultimate_x( k_order, pdt, zzt, puc, pt_u )318 CALL ultimate_x( k_order, ipl, pdt, zzt, puc, pt_u ) 295 319 ! 296 320 ENDIF … … 299 323 300 324 301 SUBROUTINE ultimate_x( k_order, pdt, pt, puc, pt_u )325 SUBROUTINE ultimate_x( k_order, ipl, pdt, pt, puc, pt_u ) 302 326 !!--------------------------------------------------------------------- 303 327 !! *** ROUTINE ultimate_x *** … … 310 334 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 311 335 !!---------------------------------------------------------------------- 312 INTEGER , INTENT(in ) :: k_order ! ocean time-step index 313 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 314 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: puc ! ice i-velocity component 315 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 316 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_u ! tracer at u-point 317 ! 318 INTEGER :: ji, jj ! dummy loop indices 336 INTEGER , INTENT(in ) :: k_order ! ocean time-step index 337 INTEGER , INTENT(in ) :: ipl ! third dimension of tracer array 338 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 339 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: puc ! ice i-velocity component 340 REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in ) :: pt ! tracer fields 341 REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT( out) :: pt_u ! tracer at u-point 342 ! 343 INTEGER :: ji, jj, jl ! dummy loop indices 319 344 REAL(wp) :: zcu, zdx2, zdx4 ! - - 320 REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu2, ztu3, ztu4 345 REAL(wp), DIMENSION(jpi,jpj) :: ztu1, ztu3 346 REAL(wp), DIMENSION(jpi,jpj,ipl) :: ztu2, ztu4 321 347 !!---------------------------------------------------------------------- 322 348 ! 323 349 ! !-- Laplacian in i-direction --! 324 DO jj = 2, jpjm1 ! First derivative (gradient) 325 DO ji = 1, fs_jpim1 326 ztu1(ji,jj) = ( pt(ji+1,jj) - pt(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 327 END DO 328 ! ! Second derivative (Laplacian) 329 DO ji = fs_2, fs_jpim1 330 ztu2(ji,jj) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj) 350 DO jl = 1, ipl 351 DO jj = 2, jpjm1 ! First derivative (gradient) 352 DO ji = 1, fs_jpim1 353 ztu1(ji,jj) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 354 END DO 355 ! ! Second derivative (Laplacian) 356 DO ji = fs_2, fs_jpim1 357 ztu2(ji,jj,jl) = ( ztu1(ji,jj) - ztu1(ji-1,jj) ) * r1_e1t(ji,jj) 358 END DO 331 359 END DO 332 360 END DO … … 334 362 ! 335 363 ! !-- BiLaplacian in i-direction --! 336 DO jj = 2, jpjm1 ! Third derivative 337 DO ji = 1, fs_jpim1 338 ztu3(ji,jj) = ( ztu2(ji+1,jj) - ztu2(ji,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 339 END DO 340 ! ! Fourth derivative 341 DO ji = fs_2, fs_jpim1 342 ztu4(ji,jj) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj) 364 DO jl = 1, ipl 365 DO jj = 2, jpjm1 ! Third derivative 366 DO ji = 1, fs_jpim1 367 ztu3(ji,jj) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 368 END DO 369 ! ! Fourth derivative 370 DO ji = fs_2, fs_jpim1 371 ztu4(ji,jj,jl) = ( ztu3(ji,jj) - ztu3(ji-1,jj) ) * r1_e1t(ji,jj) 372 END DO 343 373 END DO 344 374 END DO … … 350 380 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 351 381 ! 352 DO jj = 2, jpjm1 353 DO ji = 1, fs_jpim1 ! vector opt. 354 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & 355 & - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj) - pt(ji,jj) ) ) 382 DO jl = 1, ipl 383 DO jj = 2, jpjm1 384 DO ji = 1, fs_jpim1 ! vector opt. 385 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 386 & - SIGN( 1._wp, puc(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 387 END DO 356 388 END DO 357 389 END DO … … 359 391 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 360 392 ! 361 DO jj = 2, jpjm1 362 DO ji = 1, fs_jpim1 ! vector opt. 363 zcu = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 364 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj) + pt(ji,jj) & 365 & - zcu * ( pt(ji+1,jj) - pt(ji,jj) ) ) 393 DO jl = 1, ipl 394 DO jj = 2, jpjm1 395 DO ji = 1, fs_jpim1 ! vector opt. 396 zcu = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 397 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 398 & - zcu * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) 399 END DO 366 400 END DO 367 401 END DO … … 369 403 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 370 404 ! 371 DO jj = 2, jpjm1 372 DO ji = 1, fs_jpim1 ! vector opt. 373 zcu = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 374 zdx2 = e1u(ji,jj) * e1u(ji,jj) 405 DO jl = 1, ipl 406 DO jj = 2, jpjm1 407 DO ji = 1, fs_jpim1 ! vector opt. 408 zcu = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 409 zdx2 = e1u(ji,jj) * e1u(ji,jj) 375 410 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 376 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) & 377 & - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) & 378 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj) + ztu2(ji,jj) & 379 & - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) ) 411 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 412 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 413 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 414 & - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 415 END DO 380 416 END DO 381 417 END DO … … 383 419 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 384 420 ! 385 DO jj = 2, jpjm1 386 DO ji = 1, fs_jpim1 ! vector opt. 387 zcu = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 388 zdx2 = e1u(ji,jj) * e1u(ji,jj) 421 DO jl = 1, ipl 422 DO jj = 2, jpjm1 423 DO ji = 1, fs_jpim1 ! vector opt. 424 zcu = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 425 zdx2 = e1u(ji,jj) * e1u(ji,jj) 389 426 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 390 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) & 391 & - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) & 392 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj) + ztu2(ji,jj) & 393 & - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) ) 427 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 428 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 429 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 430 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 431 END DO 394 432 END DO 395 433 END DO … … 397 435 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 398 436 ! 399 DO jj = 2, jpjm1 400 DO ji = 1, fs_jpim1 ! vector opt. 401 zcu = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 402 zdx2 = e1u(ji,jj) * e1u(ji,jj) 437 DO jl = 1, ipl 438 DO jj = 2, jpjm1 439 DO ji = 1, fs_jpim1 ! vector opt. 440 zcu = puc(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 441 zdx2 = e1u(ji,jj) * e1u(ji,jj) 403 442 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 404 zdx4 = zdx2 * zdx2 405 pt_u(ji,jj) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj) + pt (ji,jj) & 406 & - zcu * ( pt (ji+1,jj) - pt (ji,jj) ) ) & 407 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj) + ztu2(ji,jj) & 408 & - 0.5_wp * zcu * ( ztu2(ji+1,jj) - ztu2(ji,jj) ) ) & 409 & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj) + ztu4(ji,jj) & 410 & - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj) - ztu4(ji,jj) ) ) ) 443 zdx4 = zdx2 * zdx2 444 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 445 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 446 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 447 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) & 448 & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) & 449 & - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 450 END DO 411 451 END DO 412 452 END DO … … 417 457 418 458 419 SUBROUTINE ultimate_y( k_order, pdt, pt, pvc, pt_v )459 SUBROUTINE ultimate_y( k_order, ipl, pdt, pt, pvc, pt_v ) 420 460 !!--------------------------------------------------------------------- 421 461 !! *** ROUTINE ultimate_y *** … … 428 468 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 429 469 !!---------------------------------------------------------------------- 430 INTEGER , INTENT(in ) :: k_order ! ocean time-step index 431 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 432 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pvc ! ice j-velocity component 433 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pt ! tracer fields 434 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pt_v ! tracer at v-point 435 ! 436 INTEGER :: ji, jj ! dummy loop indices 470 INTEGER , INTENT(in ) :: k_order ! ocean time-step index 471 INTEGER , INTENT(in ) :: ipl ! third dimension of tracer array 472 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 473 REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pvc ! ice j-velocity component 474 REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT(in ) :: pt ! tracer fields 475 REAL(wp), DIMENSION(jpi,jpj,ipl), INTENT( out) :: pt_v ! tracer at v-point 476 ! 477 INTEGER :: ji, jj, jl ! dummy loop indices 437 478 REAL(wp) :: zcv, zdy2, zdy4 ! - - 438 REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv2, ztv3, ztv4 479 REAL(wp), DIMENSION(jpi,jpj) :: ztv1, ztv3 480 REAL(wp), DIMENSION(jpi,jpj,ipl) :: ztv2, ztv4 439 481 !!---------------------------------------------------------------------- 440 482 ! 441 483 ! !-- Laplacian in j-direction --! 442 DO jj = 1, jpjm1 ! First derivative (gradient) 443 DO ji = fs_2, fs_jpim1 444 ztv1(ji,jj) = ( pt(ji,jj+1) - pt(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 445 END DO 446 END DO 447 DO jj = 2, jpjm1 ! Second derivative (Laplacian) 448 DO ji = fs_2, fs_jpim1 449 ztv2(ji,jj) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj) 484 DO jl = 1, ipl 485 DO jj = 1, jpjm1 ! First derivative (gradient) 486 DO ji = fs_2, fs_jpim1 487 ztv1(ji,jj) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 488 END DO 489 END DO 490 DO jj = 2, jpjm1 ! Second derivative (Laplacian) 491 DO ji = fs_2, fs_jpim1 492 ztv2(ji,jj,jl) = ( ztv1(ji,jj) - ztv1(ji,jj-1) ) * r1_e2t(ji,jj) 493 END DO 450 494 END DO 451 495 END DO … … 453 497 ! 454 498 ! !-- BiLaplacian in j-direction --! 455 DO jj = 1, jpjm1 ! First derivative 456 DO ji = fs_2, fs_jpim1 457 ztv3(ji,jj) = ( ztv2(ji,jj+1) - ztv2(ji,jj) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 458 END DO 459 END DO 460 DO jj = 2, jpjm1 ! Second derivative 461 DO ji = fs_2, fs_jpim1 462 ztv4(ji,jj) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj) 499 DO jl = 1, ipl 500 DO jj = 1, jpjm1 ! First derivative 501 DO ji = fs_2, fs_jpim1 502 ztv3(ji,jj) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 503 END DO 504 END DO 505 DO jj = 2, jpjm1 ! Second derivative 506 DO ji = fs_2, fs_jpim1 507 ztv4(ji,jj,jl) = ( ztv3(ji,jj) - ztv3(ji,jj-1) ) * r1_e2t(ji,jj) 508 END DO 463 509 END DO 464 510 END DO … … 469 515 ! 470 516 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 471 DO jj = 1, jpjm1 472 DO ji = fs_2, fs_jpim1 473 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & 474 & - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1) - pt(ji,jj) ) ) 517 DO jl = 1, ipl 518 DO jj = 1, jpjm1 519 DO ji = fs_2, fs_jpim1 520 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 521 & - SIGN( 1._wp, pvc(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 522 END DO 475 523 END DO 476 524 END DO 477 525 ! 478 526 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 479 DO jj = 1, jpjm1 480 DO ji = fs_2, fs_jpim1 481 zcv = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 482 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1) + pt(ji,jj) ) & 483 & - zcv * ( pt(ji,jj+1) - pt(ji,jj) ) ) 527 DO jl = 1, ipl 528 DO jj = 1, jpjm1 529 DO ji = fs_2, fs_jpim1 530 zcv = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 531 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & 532 & - zcv * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) 533 END DO 484 534 END DO 485 535 END DO … … 487 537 ! 488 538 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 489 DO jj = 1, jpjm1 490 DO ji = fs_2, fs_jpim1 491 zcv = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 492 zdy2 = e2v(ji,jj) * e2v(ji,jj) 539 DO jl = 1, ipl 540 DO jj = 1, jpjm1 541 DO ji = fs_2, fs_jpim1 542 zcv = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 543 zdy2 = e2v(ji,jj) * e2v(ji,jj) 493 544 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 494 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & 495 & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & 496 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1) + ztv2(ji,jj) & 497 & - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) ) 545 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 546 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 547 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 548 & - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 549 END DO 498 550 END DO 499 551 END DO 500 552 ! 501 553 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 502 DO jj = 1, jpjm1 503 DO ji = fs_2, fs_jpim1 504 zcv = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 505 zdy2 = e2v(ji,jj) * e2v(ji,jj) 554 DO jl = 1, ipl 555 DO jj = 1, jpjm1 556 DO ji = fs_2, fs_jpim1 557 zcv = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 558 zdy2 = e2v(ji,jj) * e2v(ji,jj) 506 559 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 507 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & 508 & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & 509 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1) + ztv2(ji,jj) & 510 & - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) ) 560 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 561 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 562 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 563 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 564 END DO 511 565 END DO 512 566 END DO 513 567 ! 514 568 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 515 DO jj = 1, jpjm1 516 DO ji = fs_2, fs_jpim1 517 zcv = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 518 zdy2 = e2v(ji,jj) * e2v(ji,jj) 569 DO jl = 1, ipl 570 DO jj = 1, jpjm1 571 DO ji = fs_2, fs_jpim1 572 zcv = pvc(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 573 zdy2 = e2v(ji,jj) * e2v(ji,jj) 519 574 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 520 zdy4 = zdy2 * zdy2 521 pt_v(ji,jj) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1) + pt (ji,jj) & 522 & - zcv * ( pt (ji,jj+1) - pt (ji,jj) ) ) & 523 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1) + ztv2(ji,jj) & 524 & - 0.5_wp * zcv * ( ztv2(ji,jj+1) - ztv2(ji,jj) ) ) & 525 & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1) + ztv4(ji,jj) & 526 & - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1) - ztv4(ji,jj) ) ) ) 575 zdy4 = zdy2 * zdy2 576 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 577 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 578 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 579 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) & 580 & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) & 581 & - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) ) 582 END DO 527 583 END DO 528 584 END DO … … 533 589 534 590 535 SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, pdt )591 SUBROUTINE nonosc_2d( ipl, pbef, paa, pbb, paft, pdt ) 536 592 !!--------------------------------------------------------------------- 537 593 !! *** ROUTINE nonosc *** … … 546 602 !! in-space based differencing for fluid 547 603 !!---------------------------------------------------------------------- 548 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 549 REAL(wp), DIMENSION (jpi,jpj), INTENT(in ) :: pbef, paft ! before & after field 550 REAL(wp), DIMENSION (jpi,jpj), INTENT(inout) :: paa, pbb ! monotonic fluxes in the 2 directions 551 ! 552 INTEGER :: ji, jj ! dummy loop indices 604 INTEGER , INTENT(in ) :: ipl ! third dimension of tracer array 605 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 606 REAL(wp), DIMENSION (jpi,jpj,ipl), INTENT(in ) :: pbef, paft ! before & after field 607 REAL(wp), DIMENSION (jpi,jpj,ipl), INTENT(inout) :: paa, pbb ! monotonic fluxes in the 2 directions 608 ! 609 INTEGER :: ji, jj, jl ! dummy loop indices 553 610 INTEGER :: ikm1 ! local integer 554 611 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zsml, z1_dt ! local scalars 555 612 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 556 REAL(wp), DIMENSION(jpi,jpj) :: zbetup, zbetdo, zbup, zbdo, zmsk, zdiv 613 REAL(wp), DIMENSION(jpi,jpj) :: zbup, zbdo, zmsk 614 REAL(wp), DIMENSION(jpi,jpj,ipl) :: zbetup, zbetdo, zdiv 557 615 !!---------------------------------------------------------------------- 558 616 ! … … 561 619 562 620 ! test on divergence 563 DO jj = 2, jpjm1 564 DO ji = fs_2, fs_jpim1 ! vector opt. 565 zdiv(ji,jj) = - ( paa(ji,jj) - paa(ji-1,jj ) & 566 & + pbb(ji,jj) - pbb(ji ,jj-1) ) 621 DO jl = 1, ipl 622 DO jj = 2, jpjm1 623 DO ji = fs_2, fs_jpim1 ! vector opt. 624 zdiv(ji,jj,jl) = - ( paa(ji,jj,jl) - paa(ji-1,jj ,jl) & 625 & + pbb(ji,jj,jl) - pbb(ji ,jj-1,jl) ) 626 END DO 567 627 END DO 568 628 END DO 569 629 CALL lbc_lnk( 'icedyn_adv_umx', zdiv, 'T', 1. ) ! Lateral boundary conditions (unchanged sign) 570 630 571 ! Determine ice masks for before and after tracers 572 WHERE( pbef(:,:) == 0._wp .AND. paft(:,:) == 0._wp .AND. zdiv(:,:) == 0._wp ) ; zmsk(:,:) = 0._wp 573 ELSEWHERE ; zmsk(:,:) = 1._wp * tmask(:,:,1) 574 END WHERE 575 576 ! Search local extrema 577 ! -------------------- 578 ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 579 ! zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ), & 580 ! & paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ) ) 581 ! zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ), & 582 ! & paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ) ) 583 zbup(:,:) = MAX( pbef(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ), & 584 & paft(:,:) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ) ) 585 zbdo(:,:) = MIN( pbef(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ), & 586 & paft(:,:) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ) ) 631 DO jl = 1, ipl 632 ! Determine ice masks for before and after tracers 633 WHERE( pbef(:,:,jl) == 0._wp .AND. paft(:,:,jl) == 0._wp .AND. zdiv(:,:,jl) == 0._wp ) 634 zmsk(:,:) = 0._wp 635 ELSEWHERE 636 zmsk(:,:) = 1._wp * tmask(:,:,1) 637 END WHERE 638 639 ! Search local extrema 640 ! -------------------- 641 ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 642 ! zbup(:,:) = MAX( pbef(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ), & 643 ! & paft(:,:) * tmask(:,:,1) - zbig * ( 1.e0 - tmask(:,:,1) ) ) 644 ! zbdo(:,:) = MIN( pbef(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ), & 645 ! & paft(:,:) * tmask(:,:,1) + zbig * ( 1.e0 - tmask(:,:,1) ) ) 646 zbup(:,:) = MAX( pbef(:,:,jl) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ), & 647 & paft(:,:,jl) * zmsk(:,:) - zbig * ( 1.e0 - zmsk(:,:) ) ) 648 zbdo(:,:) = MIN( pbef(:,:,jl) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ), & 649 & paft(:,:,jl) * zmsk(:,:) + zbig * ( 1.e0 - zmsk(:,:) ) ) 587 650 588 651 z1_dt = 1._wp / pdt … … 595 658 & zbdo(ji ,jj-1), zbdo(ji ,jj+1) ) 596 659 ! 597 zpos = MAX( 0., paa(ji-1,jj ) ) - MIN( 0., paa(ji ,jj) ) & ! positive/negative part of the flux598 & + MAX( 0., pbb(ji ,jj-1) ) - MIN( 0., pbb(ji ,jj) )599 zneg = MAX( 0., paa(ji ,jj ) ) - MIN( 0., paa(ji-1,jj) ) &600 & + MAX( 0., pbb(ji ,jj ) ) - MIN( 0., pbb(ji ,jj-1) )660 zpos = MAX( 0., paa(ji-1,jj ,jl) ) - MIN( 0., paa(ji ,jj ,jl) ) & ! positive/negative part of the flux 661 & + MAX( 0., pbb(ji ,jj-1,jl) ) - MIN( 0., pbb(ji ,jj ,jl) ) 662 zneg = MAX( 0., paa(ji ,jj ,jl) ) - MIN( 0., paa(ji-1,jj ,jl) ) & 663 & + MAX( 0., pbb(ji ,jj ,jl) ) - MIN( 0., pbb(ji ,jj-1,jl) ) 601 664 ! 602 zbt = e1e2t(ji,jj) * z1_dt ! up & down beta terms 603 zbetup(ji,jj) = ( zup - paft(ji,jj) ) / ( zpos + zsml ) * zbt 604 zbetdo(ji,jj) = ( paft(ji,jj) - zdo ) / ( zneg + zsml ) * zbt 665 zbt = e1e2t(ji,jj) * z1_dt ! up & down beta terms 666 zbetup(ji,jj,jl) = ( zup - paft(ji,jj,jl) ) / ( zpos + zsml ) * zbt 667 zbetdo(ji,jj,jl) = ( paft(ji,jj,jl) - zdo ) / ( zneg + zsml ) * zbt 668 END DO 605 669 END DO 606 670 END DO 607 671 CALL lbc_lnk_multi( 'icedyn_adv_umx', zbetup, 'T', 1., zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign) 608 672 609 ! monotonic flux in the i & j direction (paa & pbb) 610 ! ------------------------------------- 611 DO jj = 2, jpjm1 612 DO ji = 1, fs_jpim1 ! vector opt. 613 zau = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji+1,jj) ) 614 zbu = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji+1,jj) ) 615 zcu = 0.5 + SIGN( 0.5 , paa(ji,jj) ) 616 ! 617 paa(ji,jj) = paa(ji,jj) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 618 END DO 619 END DO 620 ! 621 DO jj = 1, jpjm1 622 DO ji = fs_2, fs_jpim1 ! vector opt. 623 zav = MIN( 1._wp , zbetdo(ji,jj) , zbetup(ji,jj+1) ) 624 zbv = MIN( 1._wp , zbetup(ji,jj) , zbetdo(ji,jj+1) ) 625 zcv = 0.5 + SIGN( 0.5 , pbb(ji,jj) ) 626 ! 627 pbb(ji,jj) = pbb(ji,jj) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 673 DO jl = 1, ipl 674 ! monotonic flux in the i & j direction (paa & pbb) 675 ! ------------------------------------- 676 DO jj = 2, jpjm1 677 DO ji = 1, fs_jpim1 ! vector opt. 678 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 679 zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) 680 zcu = 0.5 + SIGN( 0.5 , paa(ji,jj,jl) ) 681 ! 682 paa(ji,jj,jl) = paa(ji,jj,jl) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 683 END DO 684 END DO 685 ! 686 DO jj = 1, jpjm1 687 DO ji = fs_2, fs_jpim1 ! vector opt. 688 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 689 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) 690 zcv = 0.5 + SIGN( 0.5 , pbb(ji,jj,jl) ) 691 ! 692 pbb(ji,jj,jl) = pbb(ji,jj,jl) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 693 END DO 628 694 END DO 629 695 END DO -
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icethd.F90
r10170 r10180 212 212 END DO 213 213 214 IF( lk_mpp ) CALL mpp_ini_ice( npti , numout )215 216 214 IF( npti > 0 ) THEN ! If there is no ice, do nothing. 217 215 ! … … 250 248 ! ! --- & Move to 2D arrays --- ! 251 249 ! 252 IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ??253 250 ENDIF 254 251 ! -
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/ICE/icethd_zdf_bl99.F90
r10069 r10180 83 83 INTEGER, DIMENSION(jpij) :: jm_min ! reference number of top equation 84 84 INTEGER, DIMENSION(jpij) :: jm_max ! reference number of bottom equation 85 ! 85 86 LOGICAL, DIMENSION(jpij) :: l_T_converged ! true when T converges (per grid point) 87 ! 86 88 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system 87 89 REAL(wp) :: zg1 = 2._wp ! … … 113 115 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence 114 116 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity 117 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i_cp ! copy 115 118 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice 116 119 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice … … 201 204 ! 202 205 iconv = 0 ! number of iterations 203 zdti_max = 1000._wp ! maximal value of error on all points204 !206 ! 207 l_T_converged(:) = .FALSE. 205 208 ! !============================! 206 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) ! Iterative procedure begins ! 209 ! Convergence calculated until all sub-domain grid points have converged 210 ! Calculations keep going for all grid points until sub-domain convergence (vectorisation optimisation) 211 ! but values are not taken into account (results independant of MPI partitioning) 212 ! 213 DO WHILE ( ( .NOT. ALL (l_T_converged(1:npti)) ) .AND. iconv < iconv_max ) ! Iterative procedure begins ! 207 214 ! !============================! 208 215 iconv = iconv + 1 … … 217 224 ! 218 225 DO ji = 1, npti 219 ztcond_i (ji,0) = rcnd_i + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 )220 ztcond_i (ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 )226 ztcond_i_cp(ji,0) = rcnd_i + zbeta * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 227 ztcond_i_cp(ji,nlay_i) = rcnd_i + zbeta * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) 221 228 END DO 222 229 DO jk = 1, nlay_i-1 223 230 DO ji = 1, npti 224 ztcond_i (ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &225 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 )231 ztcond_i_cp(ji,jk) = rcnd_i + zbeta * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 232 & MIN( -epsi10, 0.5_wp * (t_i_1d(ji,jk) + t_i_1d(ji,jk+1)) - rt0 ) 226 233 END DO 227 234 END DO … … 230 237 ! 231 238 DO ji = 1, npti 232 ztcond_i (ji,0) = rcnd_i + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) &233 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 )234 ztcond_i (ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) &235 & - 0.011_wp * ( t_bo_1d(ji) - rt0 )239 ztcond_i_cp(ji,0) = rcnd_i + 0.09_wp * sz_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 240 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 241 ztcond_i_cp(ji,nlay_i) = rcnd_i + 0.09_wp * sz_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 242 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 236 243 END DO 237 244 DO jk = 1, nlay_i-1 238 245 DO ji = 1, npti 239 ztcond_i (ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / &240 & MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) &241 & - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 )246 ztcond_i_cp(ji,jk) = rcnd_i + 0.09_wp * 0.5_wp * ( sz_i_1d(ji,jk) + sz_i_1d(ji,jk+1) ) / & 247 & MIN( -epsi10, 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) & 248 & - 0.011_wp * ( 0.5_wp * ( t_i_1d (ji,jk) + t_i_1d (ji,jk+1) ) - rt0 ) 242 249 END DO 243 250 END DO 244 251 ! 245 252 ENDIF 246 ztcond_i(1:npti,:) = MAX( zkimin, ztcond_i(1:npti,:) ) 253 254 ! Variable used after iterations 255 ! Value must be frozen after convergence for MPP independance reason 256 DO ji = 1, npti 257 IF ( .NOT. l_T_converged(ji) ) & 258 ztcond_i(ji,:) = MAX( zkimin, ztcond_i_cp(ji,:) ) 259 END DO 247 260 ! 248 261 !--- G(he) : enhancement of thermal conductivity in mono-category case … … 270 283 !----------------- 271 284 !--- Snow 285 ! Variable used after iterations 286 ! Value must be frozen after convergence for MPP independance reason 272 287 DO jk = 0, nlay_s-1 273 288 DO ji = 1, npti 274 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 289 IF ( .NOT. l_T_converged(ji) ) & 290 zkappa_s(ji,jk) = zghe(ji) * rn_cnd_s * z1_h_s(ji) 275 291 END DO 276 292 END DO 277 293 DO ji = 1, npti ! Snow-ice interface 278 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 279 IF( zfac > epsi10 ) THEN 280 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 281 ELSE 282 zkappa_s(ji,nlay_s) = 0._wp 294 IF ( .NOT. l_T_converged(ji) ) THEN 295 zfac = 0.5_wp * ( ztcond_i(ji,0) * zh_s(ji) + rn_cnd_s * zh_i(ji) ) 296 IF( zfac > epsi10 ) THEN 297 zkappa_s(ji,nlay_s) = zghe(ji) * rn_cnd_s * ztcond_i(ji,0) / zfac 298 ELSE 299 zkappa_s(ji,nlay_s) = 0._wp 300 ENDIF 283 301 ENDIF 284 302 END DO 285 303 286 304 !--- Ice 305 ! Variable used after iterations 306 ! Value must be frozen after convergence for MPP independance reason 287 307 DO jk = 0, nlay_i 288 308 DO ji = 1, npti 289 zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 309 IF ( .NOT. l_T_converged(ji) ) & 310 zkappa_i(ji,jk) = zghe(ji) * ztcond_i(ji,jk) * z1_h_i(ji) 290 311 END DO 291 312 END DO 292 313 DO ji = 1, npti ! Snow-ice interface 293 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 314 IF ( .NOT. l_T_converged(ji) ) & 315 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 294 316 END DO 295 317 ! … … 326 348 ! update of the non solar flux according to the update in T_su 327 349 DO ji = 1, npti 328 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 350 ! Variable used after iterations 351 ! Value must be frozen after convergence for MPP independance reason 352 IF ( .NOT. l_T_converged(ji) ) & 353 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 329 354 END DO 330 355 … … 496 521 ! ice temperatures 497 522 DO ji = 1, npti 498 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 523 ! Variable used after iterations 524 ! Value must be frozen after convergence for MPP independance reason 525 IF ( .NOT. l_T_converged(ji) ) & 526 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 499 527 END DO 500 528 … … 502 530 DO ji = 1, npti 503 531 jk = jm - nlay_s - 1 504 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 532 IF ( .NOT. l_T_converged(ji) ) & 533 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 505 534 END DO 506 535 END DO 507 536 508 537 DO ji = 1, npti 509 ! snow temperatures 510 IF( h_s_1d(ji) > 0._wp ) THEN 511 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 512 ENDIF 513 ! surface temperature 514 ztsub(ji) = t_su_1d(ji) 515 IF( t_su_1d(ji) < rt0 ) THEN 516 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 517 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 538 ! Variables used after iterations 539 ! Value must be frozen after convergence for MPP independance reason 540 IF ( .NOT. l_T_converged(ji) ) THEN 541 ! snow temperatures 542 IF( h_s_1d(ji) > 0._wp ) THEN 543 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 544 ENDIF 545 ! surface temperature 546 ztsub(ji) = t_su_1d(ji) 547 IF( t_su_1d(ji) < rt0 ) THEN 548 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 549 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 550 ENDIF 518 551 ENDIF 519 552 END DO … … 524 557 ! check that nowhere it has started to melt 525 558 ! zdti_max is a measure of error, it has to be under zdti_bnd 526 zdti_max = 0._wp 527 DO ji = 1, npti 528 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 529 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 530 END DO 531 532 DO jk = 1, nlay_s 533 DO ji = 1, npti 534 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 535 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 536 END DO 537 END DO 538 539 DO jk = 1, nlay_i 540 DO ji = 1, npti 541 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 542 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 543 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 544 END DO 545 END DO 546 547 ! Compute spatial maximum over all errors 548 ! note that this could be optimized substantially by iterating only the non-converging points 549 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 550 ! 559 560 DO ji = 1, npti 561 562 zdti_max = 0._wp 563 564 IF ( .NOT. l_T_converged(ji) ) THEN 565 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 566 zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) ) 567 568 t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 569 zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 570 571 DO jk = 1, nlay_i 572 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 573 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 574 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 575 END DO 576 577 IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 578 579 ENDIF 580 581 END DO 582 551 583 !----------------------------------------! 552 584 ! ! … … 670 702 671 703 ! ice temperatures 672 DO ji = 1, npti 673 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 704 DO ji = 1, npti 705 ! Variable used after iterations 706 ! Value must be frozen after convergence for MPP independance reason 707 IF ( .NOT. l_T_converged(ji) ) & 708 t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 674 709 END DO 675 710 676 711 DO jm = nlay_i + nlay_s, nlay_s + 2, -1 677 712 DO ji = 1, npti 678 jk = jm - nlay_s - 1 679 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 713 IF ( .NOT. l_T_converged(ji) ) THEN 714 jk = jm - nlay_s - 1 715 t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 716 ENDIF 680 717 END DO 681 718 END DO … … 683 720 ! snow temperatures 684 721 DO ji = 1, npti 685 IF( h_s_1d(ji) > 0._wp ) THEN 686 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 722 ! Variable used after iterations 723 ! Value must be frozen after convergence for MPP independance reason 724 IF ( .NOT. l_T_converged(ji) ) THEN 725 IF( h_s_1d(ji) > 0._wp ) THEN 726 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 727 ENDIF 687 728 ENDIF 688 729 END DO … … 693 734 ! check that nowhere it has started to melt 694 735 ! zdti_max is a measure of error, it has to be under zdti_bnd 695 zdti_max = 0._wp 696 697 DO jk = 1, nlay_s 698 DO ji = 1, npti 699 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 700 zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 701 END DO 702 END DO 703 704 DO jk = 1, nlay_i 705 DO ji = 1, npti 706 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 707 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 708 zdti_max = MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 709 END DO 710 END DO 711 712 ! Compute spatial maximum over all errors 713 ! note that this could be optimized substantially by iterating only the non-converging points 714 IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 736 737 DO ji = 1, npti 738 739 zdti_max = 0._wp 740 741 IF ( .NOT. l_T_converged(ji) ) THEN 742 ! t_s 743 t_s_1d(ji,1:nlay_s) = MAX( MIN( t_s_1d(ji,1:nlay_s), rt0 ), rt0 - 100._wp ) 744 zdti_max = MAX ( zdti_max , MAXVAL( ABS( t_s_1d(ji,1:nlay_s) - ztsb(ji,1:nlay_s) ) ) ) 745 ! t_i 746 DO jk = 0, nlay_i 747 ztmelts = -rTmlt * sz_i_1d(ji,jk) + rt0 748 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelts ), rt0 - 100._wp ) 749 zdti_max = MAX ( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 750 END DO 751 752 IF ( zdti_max < zdti_bnd ) l_T_converged(ji) = .TRUE. 753 754 ENDIF 755 756 END DO 715 757 716 758 ENDIF ! k_jules -
NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/LBC/lib_mpp.F90
r10172 r10180 85 85 PUBLIC mpp_max_multiple 86 86 PUBLIC mppscatter, mppgather 87 PUBLIC mpp_ini_ ice, mpp_ini_znl87 PUBLIC mpp_ini_znl 88 88 PUBLIC mppsize 89 89 PUBLIC mppsend, mpprecv ! needed by TAM and ICB routines … … 133 133 134 134 INTEGER :: MPI_SUMDD 135 136 ! variables used in case of sea-ice137 INTEGER, PUBLIC :: ncomm_ice !: communicator made by the processors with sea-ice (public so that it can be freed in icethd)138 INTEGER :: ngrp_iworld ! group ID for the world processors (for rheology)139 INTEGER :: ngrp_ice ! group ID for the ice processors (for rheology)140 INTEGER :: ndim_rank_ice ! number of 'ice' processors141 INTEGER :: n_ice_root ! number (in the comm_ice) of proc 0 in the ice comm142 INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: nrank_ice ! dimension ndim_rank_ice143 135 144 136 ! variables used for zonal integration … … 1011 1003 1012 1004 1013 SUBROUTINE mpp_ini_ice( pindic, kumout )1014 !!----------------------------------------------------------------------1015 !! *** routine mpp_ini_ice ***1016 !!1017 !! ** Purpose : Initialize special communicator for ice areas1018 !! condition together with global variables needed in the ddmpp folding1019 !!1020 !! ** Method : - Look for ice processors in ice routines1021 !! - Put their number in nrank_ice1022 !! - Create groups for the world processors and the ice processors1023 !! - Create a communicator for ice processors1024 !!1025 !! ** output1026 !! njmppmax = njmpp for northern procs1027 !! ndim_rank_ice = number of processors with ice1028 !! nrank_ice (ndim_rank_ice) = ice processors1029 !! ngrp_iworld = group ID for the world processors1030 !! ngrp_ice = group ID for the ice processors1031 !! ncomm_ice = communicator for the ice procs.1032 !! n_ice_root = number (in the world) of proc 0 in the ice comm.1033 !!1034 !!----------------------------------------------------------------------1035 INTEGER, INTENT(in) :: pindic1036 INTEGER, INTENT(in) :: kumout ! ocean.output logical unit1037 !!1038 INTEGER :: jjproc1039 INTEGER :: ii, ierr1040 INTEGER, ALLOCATABLE, DIMENSION(:) :: kice1041 INTEGER, ALLOCATABLE, DIMENSION(:) :: zwork1042 !!----------------------------------------------------------------------1043 !1044 ALLOCATE( kice(jpnij), zwork(jpnij), STAT=ierr )1045 IF( ierr /= 0 ) THEN1046 WRITE(kumout, cform_err)1047 WRITE(kumout,*) 'mpp_ini_ice : failed to allocate 2, 1D arrays (jpnij in length)'1048 CALL mppstop1049 ENDIF1050 1051 ! Look for how many procs with sea-ice1052 !1053 kice = 01054 DO jjproc = 1, jpnij1055 IF( jjproc == narea .AND. pindic .GT. 0 ) kice(jjproc) = 11056 END DO1057 !1058 zwork = 01059 CALL MPI_ALLREDUCE( kice, zwork, jpnij, mpi_integer, mpi_sum, mpi_comm_oce, ierr )1060 ndim_rank_ice = SUM( zwork )1061 1062 ! Allocate the right size to nrank_north1063 IF( ALLOCATED ( nrank_ice ) ) DEALLOCATE( nrank_ice )1064 ALLOCATE( nrank_ice(ndim_rank_ice) )1065 !1066 ii = 01067 nrank_ice = 01068 DO jjproc = 1, jpnij1069 IF( zwork(jjproc) == 1) THEN1070 ii = ii + 11071 nrank_ice(ii) = jjproc -11072 ENDIF1073 END DO1074 1075 ! Create the world group1076 CALL MPI_COMM_GROUP( mpi_comm_oce, ngrp_iworld, ierr )1077 1078 ! Create the ice group from the world group1079 CALL MPI_GROUP_INCL( ngrp_iworld, ndim_rank_ice, nrank_ice, ngrp_ice, ierr )1080 1081 ! Create the ice communicator , ie the pool of procs with sea-ice1082 CALL MPI_COMM_CREATE( mpi_comm_oce, ngrp_ice, ncomm_ice, ierr )1083 1084 ! Find proc number in the world of proc 0 in the north1085 ! The following line seems to be useless, we just comment & keep it as reminder1086 ! CALL MPI_GROUP_TRANSLATE_RANKS(ngrp_ice,1,0,ngrp_iworld,n_ice_root,ierr)1087 !1088 CALL MPI_GROUP_FREE(ngrp_ice, ierr)1089 CALL MPI_GROUP_FREE(ngrp_iworld, ierr)1090 1091 DEALLOCATE(kice, zwork)1092 !1093 END SUBROUTINE mpp_ini_ice1094 1095 1096 1005 SUBROUTINE mpp_ini_znl( kumout ) 1097 1006 !!---------------------------------------------------------------------- … … 1662 1571 LOGICAL, PUBLIC, PARAMETER :: lk_mpp = .FALSE. !: mpp flag 1663 1572 LOGICAL, PUBLIC :: ln_nnogather !: namelist control of northfold comms (needed here in case "key_mpp_mpi" is not used) 1664 INTEGER :: ncomm_ice1665 1573 INTEGER, PUBLIC :: mpi_comm_oce ! opa local communicator 1666 1574 !!---------------------------------------------------------------------- … … 1815 1723 STOP ! non MPP case, just stop the run 1816 1724 END SUBROUTINE mppstop 1817 1818 SUBROUTINE mpp_ini_ice( kcom, knum )1819 INTEGER :: kcom, knum1820 WRITE(*,*) 'mpp_ini_ice: You should not have seen this print! error?', kcom, knum1821 END SUBROUTINE mpp_ini_ice1822 1725 1823 1726 SUBROUTINE mpp_ini_znl( knum )
Note: See TracChangeset
for help on using the changeset viewer.
