Changeset 10679 for NEMO/branches
- Timestamp:
- 2019-02-14T14:11:43+01:00 (5 years ago)
- Location:
- NEMO/branches/2019/fix_ticket2238_solution2
- Files:
-
- 7 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icb_oce.F90
r10425 r10679 63 63 REAL(wp) :: uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi ! properties of iceberg environment 64 64 REAL(wp) :: mass_of_bits, heat_density 65 REAL(wp), DIMENSION(4) :: xRK1, xRK2, xRK3, xRK4 66 REAL(wp), DIMENSION(4) :: yRK1, yRK2, yRK3, yRK4 65 67 END TYPE point 66 68 … … 85 87 ! Extra arrays with bigger halo, needed when interpolating forcing onto iceberg position 86 88 ! particularly for MPP when iceberg can lie inside T grid but outside U, V, or f grid 87 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: uo_e, vo_e 88 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ff_e, tt_e, fr_e, hicth 89 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ua_e, va_e 90 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ssh_e 91 #if defined key_si3 || defined key_cice 92 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ui_e, vi_e 93 #endif 89 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hicth 94 90 95 91 !!gm almost all those PARAM ARE defined in NEMO … … 168 164 icb_alloc = icb_alloc + ill 169 165 ! 170 ! expanded arrays for bilinear interpolation 171 ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) , & 172 & vo_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) , & 173 #if defined key_si3 || defined key_cice 174 & ui_e(0:jpi+1,0:jpj+1) , & 175 & vi_e(0:jpi+1,0:jpj+1) , & 176 #endif 177 & ff_e(0:jpi+1,0:jpj+1) , fr_e(0:jpi+1,0:jpj+1) , & 178 & tt_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) , & 179 & hicth(0:jpi+1,0:jpj+1), & 180 & first_width(nclasses) , first_length(nclasses) , & 166 ALLOCATE( first_width(nclasses) , first_length(nclasses) , & 181 167 & src_calving (jpi,jpj) , & 182 168 & src_calving_hflx(jpi,jpj) , STAT=ill) -
NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbdyn.F90
r10570 r10679 18 18 USE icbutl ! iceberg utility routines 19 19 USE icbdia ! iceberg budget routines 20 USE icblbc ! iceberg: lateral boundary routines (including mpp) 21 ! 22 USE in_out_manager 23 USE lib_mpp ! massively parallel library 20 24 21 25 IMPLICIT NONE … … 68 72 zdt_6 = zdt / 6._wp 69 73 70 berg => first_berg ! start from the first berg 71 ! 72 DO WHILE ( ASSOCIATED(berg) ) !== loop over all bergs ==! 74 IF( ASSOCIATED(first_berg) ) THEN 75 berg => first_berg ! start from the first berg 73 76 ! 74 pt => berg%current_point 75 76 ll_bounced = .FALSE. 77 78 79 ! STEP 1 ! 80 ! ====== ! 81 zxi1 = pt%xi ; zuvel1 = pt%uvel !** X1 in (i,j) ; V1 in m/s 82 zyj1 = pt%yj ; zvvel1 = pt%vvel 83 84 85 ! !** A1 = A(X1,V1) 86 CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1, & 87 & zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 ) 88 ! 89 zu1 = zuvel1 / ze1 !** V1 in d(i,j)/dt 90 zv1 = zvvel1 / ze2 91 92 ! STEP 2 ! 93 ! ====== ! 94 ! !** X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 95 ! position using di/dt & djdt ! V2 in m/s 96 zxi2 = zxi1 + zdt_2 * zu1 ; zuvel2 = zuvel1 + zdt_2 * zax1 97 zyj2 = zyj1 + zdt_2 * zv1 ; zvvel2 = zvvel1 + zdt_2 * zay1 98 ! 99 CALL icb_ground( zxi2, zxi1, zu1, & 100 & zyj2, zyj1, zv1, ll_bounced ) 101 102 ! !** A2 = A(X2,V2) 103 CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2, & 104 & zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 ) 105 ! 106 zu2 = zuvel2 / ze1 !** V2 in d(i,j)/dt 107 zv2 = zvvel2 / ze2 108 ! 109 ! STEP 3 ! 110 ! ====== ! 111 ! !** X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2; A3=A(X3) 112 zxi3 = zxi1 + zdt_2 * zu2 ; zuvel3 = zuvel1 + zdt_2 * zax2 113 zyj3 = zyj1 + zdt_2 * zv2 ; zvvel3 = zvvel1 + zdt_2 * zay2 114 ! 115 CALL icb_ground( zxi3, zxi1, zu3, & 116 & zyj3, zyj1, zv3, ll_bounced ) 117 118 ! !** A3 = A(X3,V3) 119 CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3, & 120 & zyj3, ze2, zvvel3, zvvel1, zay3, zdt ) 121 ! 122 zu3 = zuvel3 / ze1 !** V3 in d(i,j)/dt 123 zv3 = zvvel3 / ze2 124 125 ! STEP 4 ! 126 ! ====== ! 127 ! !** X4 = X1+dt*V3 ; V4 = V1+dt*A3 128 zxi4 = zxi1 + zdt * zu3 ; zuvel4 = zuvel1 + zdt * zax3 129 zyj4 = zyj1 + zdt * zv3 ; zvvel4 = zvvel1 + zdt * zay3 130 131 CALL icb_ground( zxi4, zxi1, zu4, & 132 & zyj4, zyj1, zv4, ll_bounced ) 133 134 ! !** A4 = A(X4,V4) 135 CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4, & 136 & zyj4, ze2, zvvel4, zvvel1, zay4, zdt ) 137 138 zu4 = zuvel4 / ze1 !** V4 in d(i,j)/dt 139 zv4 = zvvel4 / ze2 140 141 ! FINAL STEP ! 142 ! ========== ! 143 ! !** Xn = X1+dt*(V1+2*V2+2*V3+V4)/6 144 ! !** Vn = V1+dt*(A1+2*A2+2*A3+A4)/6 145 zxi_n = pt%xi + zdt_6 * ( zu1 + 2.*(zu2 + zu3 ) + zu4 ) 146 zyj_n = pt%yj + zdt_6 * ( zv1 + 2.*(zv2 + zv3 ) + zv4 ) 147 zuvel_n = pt%uvel + zdt_6 * ( zax1 + 2.*(zax2 + zax3) + zax4 ) 148 zvvel_n = pt%vvel + zdt_6 * ( zay1 + 2.*(zay2 + zay3) + zay4 ) 149 150 CALL icb_ground( zxi_n, zxi1, zuvel_n, & 151 & zyj_n, zyj1, zvvel_n, ll_bounced ) 152 153 pt%uvel = zuvel_n !** save in berg structure 154 pt%vvel = zvvel_n 155 pt%xi = zxi_n 156 pt%yj = zyj_n 157 158 ! update actual position 159 pt%lon = icb_utl_bilin_x(glamt, pt%xi, pt%yj ) 160 pt%lat = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' ) 161 162 berg => berg%next ! switch to the next berg 163 ! 164 END DO !== end loop over all bergs ==! 77 DO WHILE ( ASSOCIATED(berg) ) !== loop over all bergs ==! 78 ! 79 pt => berg%current_point 80 ! 81 ! STEP 1 ! 82 ! ====== ! 83 zxi1 = pt%xi ; zuvel1 = pt%uvel !** X1 in (i,j) ; V1 in m/s 84 zyj1 = pt%yj ; zvvel1 = pt%vvel 85 ! 86 ! !** A1 = A(X1,V1) 87 CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1, & 88 & zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 ) 89 ! 90 zu1 = zuvel1 / ze1 !** V1 in d(i,j)/dt 91 zv1 = zvvel1 / ze2 92 ! 93 ! STEP 2 ! 94 ! ====== ! 95 ! !** X2 = X1+dt/2*V1 ; V2 = V1+dt/2*A1 96 ! position using di/dt & djdt ! V2 in m/s 97 zxi2 = zxi1 + zdt_2 * zu1 ; zuvel2 = zuvel1 + zdt_2 * zax1 98 zyj2 = zyj1 + zdt_2 * zv1 ; zvvel2 = zvvel1 + zdt_2 * zay1 99 ! 100 ! check ground iceberg 101 CALL icb_ground( zxi2, zxi1, zu1, & 102 & zyj2, zyj1, zv1, ll_bounced ) 103 ! 104 ! SAVE intermediaite properties 105 pt%xRK1 = (/ zxi1, zuvel1, zu1, zax1 /); pt%yRK1 = (/ zyj1, zvvel1, zv1, zay1 /); 106 pt%xRK2 = (/ zxi2, zuvel2, 0.0, 0.0 /); pt%yRK2 = (/ zyj2, zvvel2, 0.0, 0.0 /); 107 pt%xi = zxi2 ; pt%yj = zyj2 ; 108 ! 109 berg => berg%next ! switch to the next berg 110 ! 111 END DO 112 END IF 113 ! 114 ! EXCHANGE BERG (in case first step fall out of the inner domain) 115 IF( lk_mpp ) THEN ; CALL icb_lbc_mpp() ! Send bergs to other PEs 116 ELSE ; CALL icb_lbc() ! Deal with any cyclic boundaries in non-mpp case 117 ENDIF 118 ! 119 IF( ASSOCIATED(first_berg) ) THEN 120 berg => first_berg ! start from the first berg 121 DO WHILE ( ASSOCIATED(berg) ) !== loop over all bergs ==! 122 ! 123 pt => berg%current_point 124 ! 125 ! LOAD intermediaite properties 126 zxi1 = pt%xRK1(1) ; zuvel1 = pt%xRK1(2) ; zu1 = pt%xRK1(3) ; zax1 = pt%xRK1(4) 127 zxi2 = pt%xRK2(1) ; zuvel2 = pt%xRK2(2) ; zu2 = pt%xRK2(3) ; zax2 = pt%xRK2(4) 128 zyj1 = pt%yRK1(1) ; zvvel1 = pt%yRK1(2) ; zv1 = pt%yRK1(3) ; zay1 = pt%yRK1(4) 129 zyj2 = pt%yRK2(1) ; zvvel2 = pt%yRK2(2) ; zv2 = pt%yRK2(3) ; zay2 = pt%yRK2(4) 130 ! 131 ! !** A2 = A(X2,V2) 132 CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2, & 133 & zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 ) 134 ! 135 zu2 = zuvel2 / ze1 !** V2 in d(i,j)/dt 136 zv2 = zvvel2 / ze2 137 ! 138 ! STEP 3 ! 139 ! ====== ! 140 ! !** X3 = X1+dt/2*V2 ; V3 = V1+dt/2*A2; A3=A(X3) 141 zxi3 = zxi1 + zdt_2 * zu2 ; zuvel3 = zuvel1 + zdt_2 * zax2 142 zyj3 = zyj1 + zdt_2 * zv2 ; zvvel3 = zvvel1 + zdt_2 * zay2 143 ! 144 CALL icb_ground( zxi3, zxi1, zu3, & 145 & zyj3, zyj1, zv3, ll_bounced ) 146 ! 147 ! SAVE intermediaite properties 148 pt%xRK2(3:4) = (/ zu2, zax2 /); pt%yRK2(3:4) = (/ zv2, zay2 /); 149 pt%xRK3 = (/ zxi3, zuvel3, 0.0, 0.0 /); pt%yRK3 = (/ zyj3, zvvel3, 0.0, 0.0 /); 150 pt%xi = zxi3 ; pt%yj = zyj3 ; 151 ! 152 berg => berg%next ! switch to the next berg 153 ! 154 END DO 155 END IF 156 ! 157 ! EXCHANGE BERG (in case first step fall out of the inner domain) 158 IF( lk_mpp ) THEN ; CALL icb_lbc_mpp() ! Send bergs to other PEs 159 ELSE ; CALL icb_lbc() ! Deal with any cyclic boundaries in non-mpp case 160 ENDIF 161 ! 162 IF( ASSOCIATED(first_berg) ) THEN 163 berg => first_berg ! start from the first berg 164 DO WHILE ( ASSOCIATED(berg) ) !== loop over all bergs ==! 165 ! 166 pt => berg%current_point 167 ! 168 ! LOAD intermediaite properties 169 zxi1 = pt%xRK1(1) ; zuvel1 = pt%xRK1(2) ; zu1 = pt%xRK1(3) ; zax1 = pt%xRK1(4) 170 zxi3 = pt%xRK3(1) ; zuvel3 = pt%xRK3(2) ; zu3 = pt%xRK3(3) ; zax3 = pt%xRK3(4) 171 zyj1 = pt%yRK1(1) ; zvvel1 = pt%yRK1(2) ; zv1 = pt%yRK1(3) ; zay1 = pt%yRK1(4) 172 zyj3 = pt%yRK3(1) ; zvvel3 = pt%yRK3(2) ; zv3 = pt%yRK3(3) ; zay3 = pt%yRK3(4) 173 ! 174 ! !** A3 = A(X3,V3) 175 CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3, & 176 & zyj3, ze2, zvvel3, zvvel1, zay3, zdt ) 177 ! 178 zu3 = zuvel3 / ze1 !** V3 in d(i,j)/dt 179 zv3 = zvvel3 / ze2 180 ! 181 ! STEP 4 ! 182 ! ====== ! 183 ! !** X4 = X1+dt*V3 ; V4 = V1+dt*A3 184 zxi4 = zxi1 + zdt * zu3 ; zuvel4 = zuvel1 + zdt * zax3 185 zyj4 = zyj1 + zdt * zv3 ; zvvel4 = zvvel1 + zdt * zay3 186 ! 187 CALL icb_ground( zxi4, zxi1, zu4, & 188 & zyj4, zyj1, zv4, ll_bounced ) 189 ! 190 ! SAVE intermediaite properties 191 pt%xRK3(3:4) = (/ zu3, zax3 /); pt%yRK3(3:4) = (/ zv3, zay3 /); 192 pt%xRK4 = (/ zxi4, zuvel4, 0.0, 0.0 /); pt%yRK4 = (/ zyj4, zvvel4, 0.0, 0.0 /); 193 pt%xi = zxi4 ; pt%yj = zyj4 ; 194 ! 195 berg => berg%next ! switch to the next berg 196 ! 197 END DO 198 END IF 199 ! 200 ! EXCHANGE BERG (in case first step fall out of the inner domain) 201 IF( lk_mpp ) THEN ; CALL icb_lbc_mpp() ! Send bergs to other PEs 202 ELSE ; CALL icb_lbc() ! Deal with any cyclic boundaries in non-mpp case 203 ENDIF 204 ! 205 IF( ASSOCIATED(first_berg) ) THEN 206 berg => first_berg ! start from the first berg 207 DO WHILE ( ASSOCIATED(berg) ) !== loop over all bergs ==! 208 ! 209 pt => berg%current_point 210 ! 211 ! LOAD intermediaite properties 212 zxi1 = pt%xRK1(1) ; zuvel1 = pt%xRK1(2) ; zu1 = pt%xRK1(3) ; zax1 = pt%xRK1(4) 213 zxi2 = pt%xRK2(1) ; zuvel2 = pt%xRK2(2) ; zu2 = pt%xRK2(3) ; zax2 = pt%xRK2(4) 214 zxi3 = pt%xRK3(1) ; zuvel3 = pt%xRK3(2) ; zu3 = pt%xRK3(3) ; zax3 = pt%xRK3(4) 215 zxi4 = pt%xRK4(1) ; zuvel4 = pt%xRK4(2) ; zu4 = pt%xRK4(3) ; zax4 = pt%xRK4(4) 216 zyj1 = pt%yRK1(1) ; zvvel1 = pt%yRK1(2) ; zv1 = pt%yRK1(3) ; zay1 = pt%yRK1(4) 217 zyj2 = pt%yRK2(1) ; zvvel2 = pt%yRK2(2) ; zv2 = pt%yRK2(3) ; zay2 = pt%yRK2(4) 218 zyj3 = pt%yRK3(1) ; zvvel3 = pt%yRK3(2) ; zv3 = pt%yRK3(3) ; zay3 = pt%yRK3(4) 219 zyj4 = pt%yRK4(1) ; zvvel4 = pt%yRK4(2) ; zv4 = pt%yRK4(3) ; zay4 = pt%yRK4(4) 220 ! 221 ! !** A4 = A(X4,V4) 222 CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4, & 223 & zyj4, ze2, zvvel4, zvvel1, zay4, zdt ) 224 ! 225 zu4 = zuvel4 / ze1 !** V4 in d(i,j)/dt 226 zv4 = zvvel4 / ze2 227 ! 228 ! FINAL STEP ! 229 ! ========== ! 230 ! !** Xn = X1+dt*(V1+2*V2+2*V3+V4)/6 231 ! !** Vn = V1+dt*(A1+2*A2+2*A3+A4)/6 232 zxi_n = zxi1 + zdt_6 * ( zu1 + 2._wp*(zu2 + zu3 ) + zu4 ) 233 zyj_n = zyj1 + zdt_6 * ( zv1 + 2._wp*(zv2 + zv3 ) + zv4 ) 234 zuvel_n = zuvel1 + zdt_6 * ( zax1 + 2._wp*(zax2 + zax3) + zax4 ) 235 zvvel_n = zvvel1 + zdt_6 * ( zay1 + 2._wp*(zay2 + zay3) + zay4 ) 236 ! 237 CALL icb_ground( zxi_n, zxi1, zuvel_n, & 238 & zyj_n, zyj1, zvvel_n, ll_bounced ) 239 ! 240 ! !** save in berg structure 241 pt%xi = zxi_n ; pt%yj = zyj_n ; 242 pt%uvel = zuvel_n ; pt%vvel = zvvel_n 243 ! 244 berg => berg%next ! switch to the next berg 245 ! 246 END DO 247 END IF 248 ! 249 ! EXCHANGE BERG (in case first step fall out of the inner domain) 250 IF( lk_mpp ) THEN ; CALL icb_lbc_mpp() ! Send bergs to other PEs 251 ELSE ; CALL icb_lbc() ! Deal with any cyclic boundaries in non-mpp case 252 ENDIF 253 ! 254 IF( ASSOCIATED(first_berg) ) THEN 255 berg => first_berg ! start from the first berg 256 DO WHILE ( ASSOCIATED(berg) ) !== loop over all bergs ==! 257 ! 258 pt => berg%current_point 259 ! 260 ! update actual position 261 pt%lon = icb_utl_bilin_x(glamt, pt%xi, pt%yj ) 262 pt%lat = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' ) 263 ! 264 berg => berg%next ! switch to the next berg 265 ! 266 END DO !== end loop over all bergs ==! 267 END IF 165 268 ! 166 269 END SUBROUTINE icb_dyn … … 182 285 LOGICAL , INTENT( out) :: ld_bounced ! bounced indicator 183 286 ! 184 INTEGER :: ii, ii 0185 INTEGER :: ij, ij 0287 INTEGER :: ii, iig, iig0 288 INTEGER :: ij, ijg, ijg0 186 289 INTEGER :: ibounce_method 187 290 !!---------------------------------------------------------------------- 188 291 ! 292 ! (PM) is this variable used somewhere ? 189 293 ld_bounced = .FALSE. 190 294 ! 191 ii0 = INT( pi0+0.5 ) ; ij0 = INT( pj0+0.5 ) ! initial gridpoint position (T-cell) 192 ii = INT( pi +0.5 ) ; ij = INT( pj +0.5 ) ! current - - 193 ! 194 IF( ii == ii0 .AND. ij == ij0 ) RETURN ! berg remains in the same cell 295 ! work with global position in case pi and pi0 are on different domains 296 iig0 = INT( pi0 + 0.5_wp ) ; ijg0 = INT( pj0 + 0.5_wp ) ! initial gridpoint position (T-cell) 297 iig = INT( pi + 0.5_wp ) ; ijg = INT( pj + 0.5_wp ) ! current - - 298 ! 299 IF( iig == iig0 .AND. ijg == ijg0 ) RETURN ! berg remains in the same cell 195 300 ! 196 301 ! map into current processor 197 ii0 = mi1( ii0 ) 198 ij0 = mj1( ij0 ) 199 ii = mi1( ii ) 200 ij = mj1( ij ) 302 ii = mi1( iig ) 303 ij = mj1( ijg ) 201 304 ! 202 305 IF( tmask(ii,ij,1) /= 0._wp ) RETURN ! berg reach a new t-cell, but an ocean one … … 221 324 pv = 0._wp 222 325 CASE ( 2 ) 223 IF( ii 0 /= ii) THEN326 IF( iig0 /= iig ) THEN 224 327 pi = pi0 ! return back to the initial position 225 328 pu = 0._wp ! zeroing of velocity in the direction of the grounding 226 329 ENDIF 227 IF( ij 0 /= ij) THEN330 IF( ijg0 /= ijg ) THEN 228 331 pj = pj0 ! return back to the initial position 229 332 pv = 0._wp ! zeroing of velocity in the direction of the grounding -
NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icblbc.F90
r10570 r10679 56 56 57 57 INTEGER, PARAMETER, PRIVATE :: jp_delta_buf = 25 ! Size by which to increment buffers 58 INTEGER, PARAMETER, PRIVATE :: jp_buffer_width = 15+nkounts ! items to store for each berg58 INTEGER, PARAMETER, PRIVATE :: jp_buffer_width = 47+nkounts ! items to store for each berg 59 59 60 60 #endif … … 90 90 DO WHILE( ASSOCIATED(this) ) 91 91 pt => this%current_point 92 iine = INT( pt%xi + 0.5 )92 iine = INT( pt%xi + 0.5_wp ) 93 93 IF( iine > mig(nicbei) ) THEN 94 94 pt%xi = ricb_right + MOD(pt%xi, 1._wp ) - 1._wp … … 125 125 DO WHILE( ASSOCIATED(this) ) 126 126 pt => this%current_point 127 ijne = INT( pt%yj + 0.5 )127 ijne = INT( pt%yj + 0.5_wp ) 128 128 IF( ijne .GT. mjg(nicbej) ) THEN 129 129 ! 130 iine = INT( pt%xi + 0.5 )130 iine = INT( pt%xi + 0.5_wp ) 131 131 ipts = nicbfldpts (mi1(iine)) 132 132 ! … … 232 232 DO WHILE (ASSOCIATED(this)) 233 233 pt => this%current_point 234 iine = INT( pt%xi + 0.5 )234 iine = INT( pt%xi + 0.5_wp ) 235 235 IF( ipe_E >= 0 .AND. iine > mig(nicbei) ) THEN 236 236 tmpberg => this … … 242 242 ENDIF 243 243 ! deal with periodic case 244 tmpberg%current_point%xi = ricb_right + MOD(tmpberg%current_point%xi, 1._wp ) - 1._wp 244 tmpberg%current_point%xi = ricb_right + MOD(tmpberg%current_point%xi, 1._wp ) - 1._wp 245 tmpberg%current_point%xRK1(1) = ricb_right + MOD(tmpberg%current_point%xRK1(1), 1._wp ) - 1._wp 246 tmpberg%current_point%xRK2(1) = ricb_right + MOD(tmpberg%current_point%xRK2(1), 1._wp ) - 1._wp 247 tmpberg%current_point%xRK3(1) = ricb_right + MOD(tmpberg%current_point%xRK3(1), 1._wp ) - 1._wp 248 tmpberg%current_point%xRK4(1) = ricb_right + MOD(tmpberg%current_point%xRK4(1), 1._wp ) - 1._wp 245 249 ! now pack it into buffer and delete from list 246 250 CALL icb_pack_into_buffer( tmpberg, obuffer_e, ibergs_to_send_e) … … 255 259 ENDIF 256 260 ! deal with periodic case 257 tmpberg%current_point%xi = ricb_left + MOD(tmpberg%current_point%xi, 1._wp ) 261 tmpberg%current_point%xi = ricb_left + MOD(tmpberg%current_point%xi, 1._wp ) 262 tmpberg%current_point%xRK1(1) = ricb_left + MOD(tmpberg%current_point%xRK1(1), 1._wp ) 263 tmpberg%current_point%xRK2(1) = ricb_left + MOD(tmpberg%current_point%xRK2(1), 1._wp ) 264 tmpberg%current_point%xRK3(1) = ricb_left + MOD(tmpberg%current_point%xRK3(1), 1._wp ) 265 tmpberg%current_point%xRK4(1) = ricb_left + MOD(tmpberg%current_point%xRK4(1), 1._wp ) 258 266 ! now pack it into buffer and delete from list 259 267 CALL icb_pack_into_buffer( tmpberg, obuffer_w, ibergs_to_send_w) … … 370 378 DO WHILE (ASSOCIATED(this)) 371 379 pt => this%current_point 372 ijne = INT( pt%yj + 0.5 )380 ijne = INT( pt%yj + 0.5_wp ) 373 381 IF( ipe_N >= 0 .AND. ijne .GT. mjg(nicbej) ) THEN 374 382 tmpberg => this … … 537 545 DO WHILE (ASSOCIATED(this)) 538 546 pt => this%current_point 539 iine = INT( pt%xi + 0.5 )540 ijne = INT( pt%yj + 0.5 )547 iine = INT( pt%xi + 0.5_wp ) 548 ijne = INT( pt%yj + 0.5_wp ) 541 549 IF( iine .LT. mig(nicbdi) .OR. & 542 550 iine .GT. mig(nicbei) .OR. & … … 544 552 ijne .GT. mjg(nicbej)) THEN 545 553 i = i + 1 546 WRITE(numicb,*) 'berg lost in halo: ', this%number(:),iine,ijne 547 WRITE(numicb,*) ' ', nimpp, njmpp 548 WRITE(numicb,*) ' ', nicbdi, nicbei, nicbdj, nicbej 554 WRITE(numicb,*) 'berg lost in halo: ', this%number(:),iine,ijne, pt%xi, pt%yj 555 WRITE(numicb,*) 'nimpp, njmpp ', nimpp, njmpp 556 WRITE(numicb,*) 'nicb di, ei, dj, ej ', nicbdi , nicbei , nicbdj , nicbej 557 WRITE(numicb,*) 'mijg di, ei, dj, ej ', mig(nicbdi), mig(nicbei), mjg(nicbdj), mjg(nicbej) 549 558 CALL flush( numicb ) 550 559 ENDIF … … 611 620 DO WHILE (ASSOCIATED(this)) 612 621 pt => this%current_point 613 iine = INT( pt%xi + 0.5 )614 ijne = INT( pt%yj + 0.5 )622 iine = INT( pt%xi + 0.5_wp ) 623 ijne = INT( pt%yj + 0.5_wp ) 615 624 iproc = nicbflddest(mi1(iine)) 616 625 IF( ijne .GT. mjg(nicbej) ) THEN … … 691 700 DO WHILE (ASSOCIATED(this)) 692 701 pt => this%current_point 693 iine = INT( pt%xi + 0.5 )694 ijne = INT( pt%yj + 0.5 )702 iine = INT( pt%xi + 0.5_wp ) 703 ijne = INT( pt%yj + 0.5_wp ) 695 704 ipts = nicbfldpts (mi1(iine)) 696 705 iproc = nicbflddest(mi1(iine)) … … 791 800 IF( .NOT. ASSOCIATED(pbuff) ) CALL icb_increase_buffer( pbuff, jp_delta_buf ) 792 801 IF( kb .GT. pbuff%size ) CALL icb_increase_buffer( pbuff, jp_delta_buf ) 802 IF( kb .GT. pbuff%size ) PRINT *, 'SHOULD NOT SEE THIS' 793 803 794 804 !! pack points into buffer … … 809 819 pbuff%data(14,kb) = berg%current_point%heat_density 810 820 811 pbuff%data(15,kb) = berg%mass_scaling 821 pbuff%data(15:18,kb) = berg%current_point%xRK1 822 pbuff%data(19:22,kb) = berg%current_point%xRK2 823 pbuff%data(23:26,kb) = berg%current_point%xRK3 824 pbuff%data(27:30,kb) = berg%current_point%xRK4 825 826 pbuff%data(31:34,kb) = berg%current_point%yRK1 827 pbuff%data(35:38,kb) = berg%current_point%yRK2 828 pbuff%data(39:42,kb) = berg%current_point%yRK3 829 pbuff%data(43:46,kb) = berg%current_point%yRK4 830 831 pbuff%data(47,kb) = berg%mass_scaling 812 832 DO k=1,nkounts 813 pbuff%data( 15+k,kb) = REAL( berg%number(k), wp )833 pbuff%data(47+k,kb) = REAL( berg%number(k), wp ) 814 834 END DO 815 835 ! … … 844 864 pt%heat_density = pbuff%data(14,kb) 845 865 846 currentberg%mass_scaling = pbuff%data(15,kb) 866 pt%xRK1 = pbuff%data(15:18,kb) 867 pt%xRK2 = pbuff%data(19:22,kb) 868 pt%xRK3 = pbuff%data(23:26,kb) 869 pt%xRK4 = pbuff%data(27:30,kb) 870 871 pt%yRK1 = pbuff%data(31:34,kb) 872 pt%yRK2 = pbuff%data(35:38,kb) 873 pt%yRK3 = pbuff%data(39:42,kb) 874 pt%yRK4 = pbuff%data(43:46,kb) 875 876 currentberg%mass_scaling = pbuff%data(47,kb) 847 877 DO ik = 1, nkounts 848 currentberg%number(ik) = INT( pbuff%data( 15+ik,kb) )878 currentberg%number(ik) = INT( pbuff%data(47+ik,kb) ) 849 879 END DO 850 880 ! -
NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbstp.F90
r10570 r10679 29 29 USE icbclv ! iceberg: calving routines 30 30 USE icbthm ! iceberg: thermodynamics routines 31 USE icblbc ! iceberg: lateral boundary routines (including mpp)32 31 USE icbtrj ! iceberg: trajectory I/O routines 33 32 USE icbdia ! iceberg: budget … … 91 90 9100 FORMAT('kt= ',i8, ' day= ',i8,' secs=',i8) 92 91 ! 93 ! !* copy nemo forcing arrays into iceberg versions with extra halo94 CALL icb_utl_copy() ! only necessary for variables not on T points95 !96 !97 92 ! !== process icebergs ==! 98 93 ! ! … … 104 99 ! !== For each berg, evolve ==! 105 100 ! 106 IF( ASSOCIATED(first_berg) ) CALL icb_dyn( kt ) ! ice berg dynamics 107 108 IF( lk_mpp ) THEN ; CALL icb_lbc_mpp() ! Send bergs to other PEs 109 ELSE ; CALL icb_lbc() ! Deal with any cyclic boundaries in non-mpp case 110 ENDIF 101 CALL icb_dyn( kt ) ! ice berg dynamics 111 102 112 103 IF( ASSOCIATED(first_berg) ) CALL icb_thm( kt ) ! Ice berg thermodynamics (melting) + rolling -
NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/ICB/icbutl.F90
r10570 r10679 31 31 PRIVATE 32 32 33 PUBLIC icb_utl_copy ! routine called in icbstp module34 33 PUBLIC icb_utl_interp ! routine called in icbdyn, icbthm modules 35 34 PUBLIC icb_utl_bilin ! routine called in icbini, icbdyn modules … … 54 53 CONTAINS 55 54 56 SUBROUTINE icb_utl_copy()57 !!----------------------------------------------------------------------58 !! *** ROUTINE icb_utl_copy ***59 !!60 !! ** Purpose : iceberg initialization.61 !!62 !! ** Method : - blah blah63 !!----------------------------------------------------------------------64 #if defined key_si365 REAL(wp), DIMENSION(jpi,jpj) :: zssh_lead_m ! ocean surface (ssh_m) if ice is not embedded66 ! ! ocean surface in leads if ice is embedded67 #endif68 ! copy nemo forcing arrays into iceberg versions with extra halo69 ! only necessary for variables not on T points70 ! and ssh which is used to calculate gradients71 72 uo_e(:,:) = 0._wp ; uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1)73 vo_e(:,:) = 0._wp ; vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1)74 ff_e(:,:) = 0._wp ; ff_e(1:jpi,1:jpj) = ff_f (:,:)75 tt_e(:,:) = 0._wp ; tt_e(1:jpi,1:jpj) = sst_m(:,:)76 fr_e(:,:) = 0._wp ; fr_e(1:jpi,1:jpj) = fr_i (:,:)77 ua_e(:,:) = 0._wp ; ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk78 va_e(:,:) = 0._wp ; va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk79 !80 CALL lbc_lnk_icb( 'icbutl', uo_e, 'U', -1._wp, 1, 1 )81 CALL lbc_lnk_icb( 'icbutl', vo_e, 'V', -1._wp, 1, 1 )82 CALL lbc_lnk_icb( 'icbutl', ff_e, 'F', +1._wp, 1, 1 )83 CALL lbc_lnk_icb( 'icbutl', ua_e, 'U', -1._wp, 1, 1 )84 CALL lbc_lnk_icb( 'icbutl', va_e, 'V', -1._wp, 1, 1 )85 CALL lbc_lnk_icb( 'icbutl', fr_e, 'T', +1._wp, 1, 1 )86 CALL lbc_lnk_icb( 'icbutl', tt_e, 'T', +1._wp, 1, 1 )87 #if defined key_si388 hicth(:,:) = 0._wp ; hicth(1:jpi,1:jpj) = hm_i (:,:)89 ui_e(:,:) = 0._wp ; ui_e(1:jpi, 1:jpj) = u_ice(:,:)90 vi_e(:,:) = 0._wp ; vi_e(1:jpi, 1:jpj) = v_ice(:,:)91 !92 ! compute ssh slope using ssh_lead if embedded93 zssh_lead_m(:,:) = ice_var_sshdyn(ssh_m, snwice_mass, snwice_mass_b)94 ssh_e(:,:) = 0._wp ; ssh_e(1:jpi, 1:jpj) = zssh_lead_m(:,:) * tmask(:,:,1)95 !96 CALL lbc_lnk_icb( 'icbutl', hicth, 'T', +1._wp, 1, 1 )97 CALL lbc_lnk_icb( 'icbutl', ui_e , 'U', -1._wp, 1, 1 )98 CALL lbc_lnk_icb( 'icbutl', vi_e , 'V', -1._wp, 1, 1 )99 #else100 ssh_e(:,:) = 0._wp ; ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1)101 #endif102 103 !! special for ssh which is used to calculate slope104 !! so fudge some numbers all the way around the boundary105 ssh_e(0 , :) = ssh_e(1 , :)106 ssh_e(jpi+1, :) = ssh_e(jpi, :)107 ssh_e(: , 0) = ssh_e(: , 1)108 ssh_e(: ,jpj+1) = ssh_e(: ,jpj)109 ssh_e(0,0) = ssh_e(1,1)110 ssh_e(jpi+1,0) = ssh_e(jpi,1)111 ssh_e(0,jpj+1) = ssh_e(1,jpj)112 ssh_e(jpi+1,jpj+1) = ssh_e(jpi,jpj)113 CALL lbc_lnk_icb( 'icbutl', ssh_e, 'T', +1._wp, 1, 1 )114 !115 END SUBROUTINE icb_utl_copy116 117 118 55 SUBROUTINE icb_utl_interp( pi, pe1, puo, pui, pua, pssh_i, & 119 56 & pj, pe2, pvo, pvi, pva, pssh_j, & … … 140 77 ! 141 78 REAL(wp) :: zcd, zmod ! local scalars 79 REAL(wp), DIMENSION(jpi,jpj) :: zssh 142 80 !!---------------------------------------------------------------------- 143 81 … … 145 83 pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 146 84 ! 147 puo = icb_utl_bilin _h( uo_e, pi, pj, 'U' ) ! ocean velocities148 pvo = icb_utl_bilin _h( vo_e, pi, pj, 'V' )149 psst = icb_utl_bilin _h( tt_e, pi, pj, 'T' ) ! SST150 pcn = icb_utl_bilin _h( fr_e, pi, pj, 'T' ) ! ice concentration151 pff = icb_utl_bilin _h( ff_e, pi, pj, 'F' ) ! Coriolis parameter152 ! 153 pua = icb_utl_bilin _h( ua_e, pi, pj, 'U' ) ! 10m wind154 pva = icb_utl_bilin _h( va_e, pi, pj, 'V' ) ! here (ua,va) are stress => rough conversion from stress to speed85 puo = icb_utl_bilin( ssu_m(:,:) * umask(:,:,1), pi, pj, 'U' ) ! ocean velocities 86 pvo = icb_utl_bilin( ssv_m(:,:) * vmask(:,:,1), pi, pj, 'V' ) 87 psst = icb_utl_bilin( sst_m(:,:), pi, pj, 'T' ) ! SST 88 pcn = icb_utl_bilin( fr_i(:,:) , pi, pj, 'T' ) ! ice concentration 89 pff = icb_utl_bilin( ff_f(:,:) , pi, pj, 'F' ) ! Coriolis parameter 90 ! 91 pua = icb_utl_bilin( utau(:,:) * umask(:,:,1), pi, pj, 'U' ) ! 10m wind 92 pva = icb_utl_bilin( vtau(:,:) * vmask(:,:,1), pi, pj, 'V' ) ! here (ua,va) are stress => rough conversion from stress to speed 155 93 zcd = 1.22_wp * 1.5e-3_wp ! air density * drag coefficient 156 94 zmod = 1._wp / MAX( 1.e-20, SQRT( zcd * SQRT( pua*pua + pva*pva) ) ) … … 159 97 160 98 #if defined key_si3 161 pui = icb_utl_bilin_h( ui_e , pi, pj, 'U' ) ! sea-ice velocities 162 pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V' ) 163 phi = icb_utl_bilin_h( hicth, pi, pj, 'T' ) ! ice thickness 99 pui = icb_utl_bilin( u_ice , pi, pj, 'U' ) ! sea-ice velocities 100 pvi = icb_utl_bilin( v_ice , pi, pj, 'V' ) 101 phi = icb_utl_bilin( hm_i , pi, pj, 'T' ) ! ice thickness 102 ! Estimate SSH gradient in i- and j-direction (centred evaluation) 103 ! compute ssh slope using ssh_lead if embedded 104 zssh(:,:) = ice_var_sshdyn(ssh_m, snwice_mass, snwice_mass_b) 164 105 #else 165 106 pui = 0._wp 166 107 pvi = 0._wp 167 108 phi = 0._wp 109 zssh(:,:) = ssh_m(:,:) 168 110 #endif 111 zssh(:,:) = zssh(:,:) * tmask(:,:,1) 169 112 170 113 ! Estimate SSH gradient in i- and j-direction (centred evaluation) 171 pssh_i = ( icb_utl_bilin _h( ssh_e, pi+0.1_wp, pj, 'T' ) - &172 & icb_utl_bilin _h( ssh_e, pi-0.1_wp, pj, 'T' ) ) / ( 0.2_wp * pe1 )173 pssh_j = ( icb_utl_bilin _h( ssh_e, pi, pj+0.1_wp, 'T' ) - &174 & icb_utl_bilin _h( ssh_e, pi, pj-0.1_wp, 'T' ) ) / ( 0.2_wp * pe2 )114 pssh_i = ( icb_utl_bilin( zssh(:,:), pi+0.1_wp, pj, 'T' ) - & 115 & icb_utl_bilin( zssh(:,:), pi-0.1_wp, pj, 'T' ) ) / ( 0.2_wp * pe1 ) 116 pssh_j = ( icb_utl_bilin( zssh(:,:), pi, pj+0.1_wp, 'T' ) - & 117 & icb_utl_bilin( zssh(:,:), pi, pj-0.1_wp, 'T' ) ) / ( 0.2_wp * pe2 ) 175 118 ! 176 119 END SUBROUTINE icb_utl_interp 177 178 179 REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type )180 !!----------------------------------------------------------------------181 !! *** FUNCTION icb_utl_bilin ***182 !!183 !! ** Purpose : bilinear interpolation at berg location depending on the grid-point type184 !! this version deals with extra halo points185 !!186 !! !!gm CAUTION an optional argument should be added to handle187 !! the slip/no-slip conditions ==>>> to be done later188 !!189 !!----------------------------------------------------------------------190 REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) :: pfld ! field to be interpolated191 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential192 CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points193 !194 INTEGER :: ii, ij ! local integer195 REAL(wp) :: zi, zj ! local real196 !!----------------------------------------------------------------------197 !198 SELECT CASE ( cd_type )199 CASE ( 'T' )200 ! note that here there is no +0.5 added201 ! since we're looking for four T points containing quadrant we're in of202 ! current T cell203 ii = MAX(1, INT( pi ))204 ij = MAX(1, INT( pj )) ! T-point205 zi = pi - REAL(ii,wp)206 zj = pj - REAL(ij,wp)207 CASE ( 'U' )208 ii = MAX(1, INT( pi-0.5 ))209 ij = MAX(1, INT( pj )) ! U-point210 zi = pi - 0.5 - REAL(ii,wp)211 zj = pj - REAL(ij,wp)212 CASE ( 'V' )213 ii = MAX(1, INT( pi ))214 ij = MAX(1, INT( pj-0.5 )) ! V-point215 zi = pi - REAL(ii,wp)216 zj = pj - 0.5 - REAL(ij,wp)217 CASE ( 'F' )218 ii = MAX(1, INT( pi-0.5 ))219 ij = MAX(1, INT( pj-0.5 )) ! F-point220 zi = pi - 0.5 - REAL(ii,wp)221 zj = pj - 0.5 - REAL(ij,wp)222 END SELECT223 !224 ! find position in this processor. Prevent near edge problems (see #1389)225 !226 IF ( ii < mig( 1 ) ) THEN ; ii = 1227 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi228 ELSE ; ii = mi1(ii)229 ENDIF230 IF ( ij < mjg( 1 ) ) THEN ; ij = 1231 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj232 ELSE ; ij = mj1(ij)233 ENDIF234 !235 IF( ii == jpi ) ii = ii-1236 IF( ij == jpj ) ij = ij-1237 !238 icb_utl_bilin_h = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) &239 & + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) * zj240 !241 END FUNCTION icb_utl_bilin_h242 243 120 244 121 REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type ) … … 270 147 zj = pj - REAL(ij,wp) 271 148 CASE ( 'U' ) 272 ii = MAX(1, INT( pi-0.5 ))273 ij = MAX(1, INT( pj )) ! U-point274 zi = pi - 0.5 - REAL(ii,wp)149 ii = MAX(1, INT( pi-0.5_wp )) 150 ij = MAX(1, INT( pj )) ! U-point 151 zi = pi - 0.5_wp - REAL(ii,wp) 275 152 zj = pj - REAL(ij,wp) 276 153 CASE ( 'V' ) 277 ii = MAX(1, INT( pi ))278 ij = MAX(1, INT( pj-0.5 )) ! V-point154 ii = MAX(1, INT( pi )) 155 ij = MAX(1, INT( pj-0.5_wp )) ! V-point 279 156 zi = pi - REAL(ii,wp) 280 zj = pj - 0.5 - REAL(ij,wp)157 zj = pj - 0.5_wp - REAL(ij,wp) 281 158 CASE ( 'F' ) 282 ii = MAX(1, INT( pi-0.5 ))283 ij = MAX(1, INT( pj-0.5 )) ! F-point284 zi = pi - 0.5 - REAL(ii,wp)285 zj = pj - 0.5 - REAL(ij,wp)159 ii = MAX(1, INT( pi-0.5_wp )) 160 ij = MAX(1, INT( pj-0.5_wp )) ! F-point 161 zi = pi - 0.5_wp - REAL(ii,wp) 162 zj = pj - 0.5_wp - REAL(ij,wp) 286 163 END SELECT 287 164 ! 288 ! find position in this processor. Prevent near edge problems (see #1389) 289 IF ( ii < mig( 1 ) ) THEN ; ii = 1 290 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 291 ELSE ; ii = mi1(ii) 292 ENDIF 293 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 294 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 295 ELSE ; ij = mj1(ij) 296 ENDIF 297 ! 298 IF( ii == jpi ) ii = ii-1 299 IF( ij == jpj ) ij = ij-1 165 ! find position in this processor. 166 ii = mi1(ii) ; ij = mj1(ij) 300 167 ! 301 168 icb_utl_bilin = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) & … … 333 200 zj = pj - REAL(ij,wp) 334 201 ! 335 ! find position in this processor. Prevent near edge problems (see #1389) 336 IF ( ii < mig( 1 ) ) THEN ; ii = 1 337 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 338 ELSE ; ii = mi1(ii) 339 ENDIF 340 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 341 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 342 ELSE ; ij = mj1(ij) 343 ENDIF 344 ! 345 IF( ii == jpi ) ii = ii-1 346 IF( ij == jpj ) ij = ij-1 202 ! find position in this processor. 203 ii = mi1(ii) ; ij = mj1(ij) 347 204 ! 348 205 z4(1) = pfld(ii ,ij ) … … 392 249 zi = pi - REAL(ii,wp) !!gm use here mig, mjg arrays 393 250 zj = pj - REAL(ij,wp) 394 395 ! find position in this processor. Prevent near edge problems (see #1389) 396 IF ( ii < mig( 1 ) ) THEN ; ii = 1 397 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 398 ELSE ; ii = mi1(ii) 399 ENDIF 400 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 401 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 402 ELSE ; ij = mj1(ij) 403 ENDIF 404 ! 405 IF( ii == jpi ) ii = ii-1 406 IF( ij == jpj ) ij = ij-1 251 ! 252 ! find position in this processor. 253 ii = mi1(ii) ; ij = mj1(ij) 407 254 ! 408 255 IF( 0.0_wp <= zi .AND. zi < 0.5_wp ) THEN -
NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/LBC/lbclnk.F90
r10425 r10679 41 41 END INTERFACE 42 42 ! 43 INTERFACE lbc_lnk_icb44 MODULE PROCEDURE mpp_lnk_2d_icb45 END INTERFACE46 43 47 44 PUBLIC lbc_lnk ! ocean/ice lateral boundary conditions 48 45 PUBLIC lbc_lnk_multi ! modified ocean/ice lateral boundary conditions 49 46 PUBLIC lbc_bdy_lnk ! ocean lateral BDY boundary conditions 50 PUBLIC lbc_lnk_icb ! iceberg lateral boundary conditions51 47 52 48 !!---------------------------------------------------------------------- … … 93 89 END INTERFACE 94 90 ! 95 INTERFACE lbc_lnk_icb96 MODULE PROCEDURE lbc_lnk_2d_icb97 END INTERFACE98 91 99 92 PUBLIC lbc_lnk ! ocean/ice lateral boundary conditions 100 93 PUBLIC lbc_lnk_multi ! modified ocean/ice lateral boundary conditions 101 94 PUBLIC lbc_bdy_lnk ! ocean lateral BDY boundary conditions 102 PUBLIC lbc_lnk_icb ! iceberg lateral boundary conditions103 95 104 96 !!---------------------------------------------------------------------- … … 214 206 215 207 216 !!gm This routine should be removed with an optional halos size added in argument of generic routines217 218 SUBROUTINE lbc_lnk_2d_icb( cdname, pt2d, cd_type, psgn, ki, kj )219 !!----------------------------------------------------------------------220 CHARACTER(len=*) , INTENT(in ) :: cdname ! name of the calling subroutine221 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt2d ! 2D array on which the lbc is applied222 CHARACTER(len=1) , INTENT(in ) :: cd_type ! nature of pt3d grid-points223 REAL(wp) , INTENT(in ) :: psgn ! sign used across north fold224 INTEGER , INTENT(in ) :: ki, kj ! sizes of extra halo (not needed in non-mpp)225 !!----------------------------------------------------------------------226 CALL lbc_lnk_2d( cdname, pt2d, cd_type, psgn )227 END SUBROUTINE lbc_lnk_2d_icb228 !!gm end229 208 230 209 #endif -
NEMO/branches/2019/fix_ticket2238_solution2/src/OCE/LBC/lib_mpp.F90
r10538 r10679 41 41 !! mynode : indentify the processor unit 42 42 !! mpp_lnk : interface (defined in lbclnk) for message passing of 2d or 3d arrays (mpp_lnk_2d, mpp_lnk_3d) 43 !! mpp_lnk_icb : interface for message passing of 2d arrays with extra halo for icebergs (mpp_lnk_2d_icb)44 43 !! mpprecv : 45 44 !! mppsend : … … 54 53 !! mppstop : 55 54 !! mpp_ini_north : initialisation of north fold 56 !! mpp_lbc_north_icb : alternative to mpp_nfd for extra outer halo with icebergs57 55 !!---------------------------------------------------------------------- 58 56 USE dom_oce ! ocean space and time domain … … 80 78 PUBLIC mynode, mppstop, mppsync, mpp_comm_free 81 79 PUBLIC mpp_ini_north 82 PUBLIC mpp_lnk_2d_icb83 PUBLIC mpp_lbc_north_icb84 80 PUBLIC mpp_min, mpp_max, mpp_sum, mpp_minloc, mpp_maxloc 85 81 PUBLIC mpp_delay_max, mpp_delay_sum, mpp_delay_rcv … … 1185 1181 ! 1186 1182 END SUBROUTINE DDPDD_MPI 1187 1188 1189 SUBROUTINE mpp_lbc_north_icb( pt2d, cd_type, psgn, kextj)1190 !!---------------------------------------------------------------------1191 !! *** routine mpp_lbc_north_icb ***1192 !!1193 !! ** Purpose : Ensure proper north fold horizontal bondary condition1194 !! in mpp configuration in case of jpn1 > 1 and for 2d1195 !! array with outer extra halo1196 !!1197 !! ** Method : North fold condition and mpp with more than one proc1198 !! in i-direction require a specific treatment. We gather1199 !! the 4+kextj northern lines of the global domain on 11200 !! processor and apply lbc north-fold on this sub array.1201 !! Then we scatter the north fold array back to the processors.1202 !! This routine accounts for an extra halo with icebergs1203 !! and assumes ghost rows and columns have been suppressed.1204 !!1205 !!----------------------------------------------------------------------1206 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt2d ! 2D array with extra halo1207 CHARACTER(len=1) , INTENT(in ) :: cd_type ! nature of pt3d grid-points1208 ! ! = T , U , V , F or W -points1209 REAL(wp) , INTENT(in ) :: psgn ! = -1. the sign change across the1210 !! ! north fold, = 1. otherwise1211 INTEGER , INTENT(in ) :: kextj ! Extra halo width at north fold1212 !1213 INTEGER :: ji, jj, jr1214 INTEGER :: ierr, itaille, ildi, ilei, iilb1215 INTEGER :: ipj, ij, iproc1216 !1217 REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: ztab_e, znorthloc_e1218 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: znorthgloio_e1219 !!----------------------------------------------------------------------1220 !1221 ipj=41222 ALLOCATE( ztab_e(jpiglo, 1-kextj:ipj+kextj) , &1223 & znorthloc_e(jpimax, 1-kextj:ipj+kextj) , &1224 & znorthgloio_e(jpimax, 1-kextj:ipj+kextj,jpni) )1225 !1226 ztab_e(:,:) = 0._wp1227 znorthloc_e(:,:) = 0._wp1228 !1229 ij = 1 - kextj1230 ! put the last ipj+2*kextj lines of pt2d into znorthloc_e1231 DO jj = jpj - ipj + 1 - kextj , jpj + kextj1232 znorthloc_e(1:jpi,ij)=pt2d(1:jpi,jj)1233 ij = ij + 11234 END DO1235 !1236 itaille = jpimax * ( ipj + 2*kextj )1237 !1238 IF( ln_timing ) CALL tic_tac(.TRUE.)1239 CALL MPI_ALLGATHER( znorthloc_e(1,1-kextj) , itaille, MPI_DOUBLE_PRECISION, &1240 & znorthgloio_e(1,1-kextj,1), itaille, MPI_DOUBLE_PRECISION, &1241 & ncomm_north, ierr )1242 !1243 IF( ln_timing ) CALL tic_tac(.FALSE.)1244 !1245 DO jr = 1, ndim_rank_north ! recover the global north array1246 iproc = nrank_north(jr) + 11247 ildi = nldit (iproc)1248 ilei = nleit (iproc)1249 iilb = nimppt(iproc)1250 DO jj = 1-kextj, ipj+kextj1251 DO ji = ildi, ilei1252 ztab_e(ji+iilb-1,jj) = znorthgloio_e(ji,jj,jr)1253 END DO1254 END DO1255 END DO1256 1257 ! 2. North-Fold boundary conditions1258 ! ----------------------------------1259 CALL lbc_nfd( ztab_e(:,1-kextj:ipj+kextj), cd_type, psgn, kextj )1260 1261 ij = 1 - kextj1262 !! Scatter back to pt2d1263 DO jj = jpj - ipj + 1 - kextj , jpj + kextj1264 DO ji= 1, jpi1265 pt2d(ji,jj) = ztab_e(ji+nimpp-1,ij)1266 END DO1267 ij = ij +11268 END DO1269 !1270 DEALLOCATE( ztab_e, znorthloc_e, znorthgloio_e )1271 !1272 END SUBROUTINE mpp_lbc_north_icb1273 1274 1275 SUBROUTINE mpp_lnk_2d_icb( cdname, pt2d, cd_type, psgn, kexti, kextj )1276 !!----------------------------------------------------------------------1277 !! *** routine mpp_lnk_2d_icb ***1278 !!1279 !! ** Purpose : Message passing management for 2d array (with extra halo for icebergs)1280 !! This routine receives a (1-kexti:jpi+kexti,1-kexti:jpj+kextj)1281 !! array (usually (0:jpi+1, 0:jpj+1)) from lbc_lnk_icb calls.1282 !!1283 !! ** Method : Use mppsend and mpprecv function for passing mask1284 !! between processors following neighboring subdomains.1285 !! domain parameters1286 !! jpi : first dimension of the local subdomain1287 !! jpj : second dimension of the local subdomain1288 !! kexti : number of columns for extra outer halo1289 !! kextj : number of rows for extra outer halo1290 !! nbondi : mark for "east-west local boundary"1291 !! nbondj : mark for "north-south local boundary"1292 !! noea : number for local neighboring processors1293 !! nowe : number for local neighboring processors1294 !! noso : number for local neighboring processors1295 !! nono : number for local neighboring processors1296 !!----------------------------------------------------------------------1297 CHARACTER(len=*) , INTENT(in ) :: cdname ! name of the calling subroutine1298 REAL(wp), DIMENSION(1-kexti:jpi+kexti,1-kextj:jpj+kextj), INTENT(inout) :: pt2d ! 2D array with extra halo1299 CHARACTER(len=1) , INTENT(in ) :: cd_type ! nature of ptab array grid-points1300 REAL(wp) , INTENT(in ) :: psgn ! sign used across the north fold1301 INTEGER , INTENT(in ) :: kexti ! extra i-halo width1302 INTEGER , INTENT(in ) :: kextj ! extra j-halo width1303 !1304 INTEGER :: jl ! dummy loop indices1305 INTEGER :: imigr, iihom, ijhom ! local integers1306 INTEGER :: ipreci, iprecj ! - -1307 INTEGER :: ml_req1, ml_req2, ml_err ! for key_mpi_isend1308 INTEGER, DIMENSION(MPI_STATUS_SIZE) :: ml_stat ! for key_mpi_isend1309 !!1310 REAL(wp), DIMENSION(1-kexti:jpi+kexti,nn_hls+kextj,2) :: r2dns, r2dsn1311 REAL(wp), DIMENSION(1-kextj:jpj+kextj,nn_hls+kexti,2) :: r2dwe, r2dew1312 !!----------------------------------------------------------------------1313 1314 ipreci = nn_hls + kexti ! take into account outer extra 2D overlap area1315 iprecj = nn_hls + kextj1316 1317 IF( narea == 1 .AND. numcom == -1 ) CALL mpp_report( cdname, 1, 1, 1, ld_lbc = .TRUE. )1318 1319 ! 1. standard boundary treatment1320 ! ------------------------------1321 ! Order matters Here !!!!1322 !1323 ! ! East-West boundaries1324 ! !* Cyclic east-west1325 IF( l_Iperio ) THEN1326 pt2d(1-kexti: 1 ,:) = pt2d(jpim1-kexti: jpim1 ,:) ! east1327 pt2d( jpi :jpi+kexti,:) = pt2d( 2 :2+kexti,:) ! west1328 !1329 ELSE !* closed1330 IF( .NOT. cd_type == 'F' ) pt2d( 1-kexti :nn_hls ,:) = 0._wp ! east except at F-point1331 pt2d(jpi-nn_hls+1:jpi+kexti,:) = 0._wp ! west1332 ENDIF1333 ! ! North-South boundaries1334 IF( l_Jperio ) THEN !* cyclic (only with no mpp j-split)1335 pt2d(:,1-kextj: 1 ) = pt2d(:,jpjm1-kextj: jpjm1) ! north1336 pt2d(:, jpj :jpj+kextj) = pt2d(:, 2 :2+kextj) ! south1337 ELSE !* closed1338 IF( .NOT. cd_type == 'F' ) pt2d(:, 1-kextj :nn_hls ) = 0._wp ! north except at F-point1339 pt2d(:,jpj-nn_hls+1:jpj+kextj) = 0._wp ! south1340 ENDIF1341 !1342 1343 ! north fold treatment1344 ! -----------------------1345 IF( npolj /= 0 ) THEN1346 !1347 SELECT CASE ( jpni )1348 CASE ( 1 ) ; CALL lbc_nfd ( pt2d(1:jpi,1:jpj+kextj), cd_type, psgn, kextj )1349 CASE DEFAULT ; CALL mpp_lbc_north_icb( pt2d(1:jpi,1:jpj+kextj), cd_type, psgn, kextj )1350 END SELECT1351 !1352 ENDIF1353 1354 ! 2. East and west directions exchange1355 ! ------------------------------------1356 ! we play with the neigbours AND the row number because of the periodicity1357 !1358 SELECT CASE ( nbondi ) ! Read Dirichlet lateral conditions1359 CASE ( -1, 0, 1 ) ! all exept 2 (i.e. close case)1360 iihom = jpi-nreci-kexti1361 DO jl = 1, ipreci1362 r2dew(:,jl,1) = pt2d(nn_hls+jl,:)1363 r2dwe(:,jl,1) = pt2d(iihom +jl,:)1364 END DO1365 END SELECT1366 !1367 ! ! Migrations1368 imigr = ipreci * ( jpj + 2*kextj )1369 !1370 IF( ln_timing ) CALL tic_tac(.TRUE.)1371 !1372 SELECT CASE ( nbondi )1373 CASE ( -1 )1374 CALL mppsend( 2, r2dwe(1-kextj,1,1), imigr, noea, ml_req1 )1375 CALL mpprecv( 1, r2dew(1-kextj,1,2), imigr, noea )1376 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)1377 CASE ( 0 )1378 CALL mppsend( 1, r2dew(1-kextj,1,1), imigr, nowe, ml_req1 )1379 CALL mppsend( 2, r2dwe(1-kextj,1,1), imigr, noea, ml_req2 )1380 CALL mpprecv( 1, r2dew(1-kextj,1,2), imigr, noea )1381 CALL mpprecv( 2, r2dwe(1-kextj,1,2), imigr, nowe )1382 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)1383 IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)1384 CASE ( 1 )1385 CALL mppsend( 1, r2dew(1-kextj,1,1), imigr, nowe, ml_req1 )1386 CALL mpprecv( 2, r2dwe(1-kextj,1,2), imigr, nowe )1387 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)1388 END SELECT1389 !1390 IF( ln_timing ) CALL tic_tac(.FALSE.)1391 !1392 ! ! Write Dirichlet lateral conditions1393 iihom = jpi - nn_hls1394 !1395 SELECT CASE ( nbondi )1396 CASE ( -1 )1397 DO jl = 1, ipreci1398 pt2d(iihom+jl,:) = r2dew(:,jl,2)1399 END DO1400 CASE ( 0 )1401 DO jl = 1, ipreci1402 pt2d(jl-kexti,:) = r2dwe(:,jl,2)1403 pt2d(iihom+jl,:) = r2dew(:,jl,2)1404 END DO1405 CASE ( 1 )1406 DO jl = 1, ipreci1407 pt2d(jl-kexti,:) = r2dwe(:,jl,2)1408 END DO1409 END SELECT1410 1411 1412 ! 3. North and south directions1413 ! -----------------------------1414 ! always closed : we play only with the neigbours1415 !1416 IF( nbondj /= 2 ) THEN ! Read Dirichlet lateral conditions1417 ijhom = jpj-nrecj-kextj1418 DO jl = 1, iprecj1419 r2dsn(:,jl,1) = pt2d(:,ijhom +jl)1420 r2dns(:,jl,1) = pt2d(:,nn_hls+jl)1421 END DO1422 ENDIF1423 !1424 ! ! Migrations1425 imigr = iprecj * ( jpi + 2*kexti )1426 !1427 IF( ln_timing ) CALL tic_tac(.TRUE.)1428 !1429 SELECT CASE ( nbondj )1430 CASE ( -1 )1431 CALL mppsend( 4, r2dsn(1-kexti,1,1), imigr, nono, ml_req1 )1432 CALL mpprecv( 3, r2dns(1-kexti,1,2), imigr, nono )1433 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)1434 CASE ( 0 )1435 CALL mppsend( 3, r2dns(1-kexti,1,1), imigr, noso, ml_req1 )1436 CALL mppsend( 4, r2dsn(1-kexti,1,1), imigr, nono, ml_req2 )1437 CALL mpprecv( 3, r2dns(1-kexti,1,2), imigr, nono )1438 CALL mpprecv( 4, r2dsn(1-kexti,1,2), imigr, noso )1439 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)1440 IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)1441 CASE ( 1 )1442 CALL mppsend( 3, r2dns(1-kexti,1,1), imigr, noso, ml_req1 )1443 CALL mpprecv( 4, r2dsn(1-kexti,1,2), imigr, noso )1444 IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)1445 END SELECT1446 !1447 IF( ln_timing ) CALL tic_tac(.FALSE.)1448 !1449 ! ! Write Dirichlet lateral conditions1450 ijhom = jpj - nn_hls1451 !1452 SELECT CASE ( nbondj )1453 CASE ( -1 )1454 DO jl = 1, iprecj1455 pt2d(:,ijhom+jl) = r2dns(:,jl,2)1456 END DO1457 CASE ( 0 )1458 DO jl = 1, iprecj1459 pt2d(:,jl-kextj) = r2dsn(:,jl,2)1460 pt2d(:,ijhom+jl) = r2dns(:,jl,2)1461 END DO1462 CASE ( 1 )1463 DO jl = 1, iprecj1464 pt2d(:,jl-kextj) = r2dsn(:,jl,2)1465 END DO1466 END SELECT1467 !1468 END SUBROUTINE mpp_lnk_2d_icb1469 1470 1183 1471 1184 SUBROUTINE mpp_report( cdname, kpk, kpl, kpf, ld_lbc, ld_glb, ld_dlg )
Note: See TracChangeset
for help on using the changeset viewer.