- Timestamp:
- 2018-01-06T15:18:23+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ICB/icbutl.F90
r9019 r9190 9 9 !! - ! 2011-04 (Alderson) Split into separate modules 10 10 !!---------------------------------------------------------------------- 11 11 12 !!---------------------------------------------------------------------- 12 13 !! icb_utl_interp : … … 48 49 !! $Id$ 49 50 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 50 !!------------------------------------------------------------------------- 51 51 !!---------------------------------------------------------------------- 52 52 CONTAINS 53 53 … … 65 65 ! and ssh which is used to calculate gradients 66 66 67 uo_e(:,:) = 0._wp ; uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1)68 vo_e(:,:) = 0._wp ; vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1)69 ff_e(:,:) = 0._wp ; ff_e(1:jpi,1:jpj) = ff_f (:,:)70 tt_e(:,:) = 0._wp ; tt_e(1:jpi,1:jpj) = sst_m(:,:)71 fr_e(:,:) = 0._wp ; fr_e(1:jpi,1:jpj) = fr_i (:,:)72 ua_e(:,:) = 0._wp ; ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk73 va_e(:,:) = 0._wp ; va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk74 67 uo_e(:,:) = 0._wp ; uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 68 vo_e(:,:) = 0._wp ; vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 69 ff_e(:,:) = 0._wp ; ff_e(1:jpi,1:jpj) = ff_f (:,:) 70 tt_e(:,:) = 0._wp ; tt_e(1:jpi,1:jpj) = sst_m(:,:) 71 fr_e(:,:) = 0._wp ; fr_e(1:jpi,1:jpj) = fr_i (:,:) 72 ua_e(:,:) = 0._wp ; ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk 73 va_e(:,:) = 0._wp ; va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk 74 ! 75 75 CALL lbc_lnk_icb( uo_e, 'U', -1._wp, 1, 1 ) 76 76 CALL lbc_lnk_icb( vo_e, 'V', -1._wp, 1, 1 ) … … 84 84 ui_e(:,:) = 0._wp ; ui_e(1:jpi, 1:jpj) = u_ice(:,:) 85 85 vi_e(:,:) = 0._wp ; vi_e(1:jpi, 1:jpj) = v_ice(:,:) 86 CALL lbc_lnk_icb(hicth, 'T', +1._wp, 1, 1 ) 87 CALL lbc_lnk_icb( ui_e, 'U', -1._wp, 1, 1 ) 88 CALL lbc_lnk_icb( vi_e, 'V', -1._wp, 1, 1 ) 86 ! 87 CALL lbc_lnk_icb( hicth, 'T', +1._wp, 1, 1 ) 88 CALL lbc_lnk_icb( ui_e , 'U', -1._wp, 1, 1 ) 89 CALL lbc_lnk_icb( vi_e , 'V', -1._wp, 1, 1 ) 89 90 #endif 90 91 … … 149 150 150 151 #if defined key_lim3 151 pui = icb_utl_bilin_h( ui_e , pi, pj, 'U' ) ! sea-ice velocities152 pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V' )153 phi = icb_utl_bilin_h( hicth, pi, pj, 'T' ) ! ice thickness152 pui = icb_utl_bilin_h( ui_e , pi, pj, 'U' ) ! sea-ice velocities 153 pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V' ) 154 phi = icb_utl_bilin_h( hicth, pi, pj, 'T' ) ! ice thickness 154 155 #else 155 156 pui = 0._wp … … 160 161 ! Estimate SSH gradient in i- and j-direction (centred evaluation) 161 162 pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T' ) - & 162 &icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' ) ) / ( 0.2_wp * pe1 )163 & icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T' ) ) / ( 0.2_wp * pe1 ) 163 164 pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T' ) - & 164 &icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' ) ) / ( 0.2_wp * pe2 )165 & icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T' ) ) / ( 0.2_wp * pe2 ) 165 166 ! 166 167 END SUBROUTINE icb_utl_interp … … 181 182 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 182 183 CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points 184 ! 185 INTEGER :: ii, ij ! local integer 186 REAL(wp) :: zi, zj ! local real 187 !!---------------------------------------------------------------------- 188 ! 189 SELECT CASE ( cd_type ) 190 CASE ( 'T' ) 191 ! note that here there is no +0.5 added 192 ! since we're looking for four T points containing quadrant we're in of 193 ! current T cell 194 ii = MAX(1, INT( pi )) 195 ij = MAX(1, INT( pj )) ! T-point 196 zi = pi - REAL(ii,wp) 197 zj = pj - REAL(ij,wp) 198 CASE ( 'U' ) 199 ii = MAX(1, INT( pi-0.5 )) 200 ij = MAX(1, INT( pj )) ! U-point 201 zi = pi - 0.5 - REAL(ii,wp) 202 zj = pj - REAL(ij,wp) 203 CASE ( 'V' ) 204 ii = MAX(1, INT( pi )) 205 ij = MAX(1, INT( pj-0.5 )) ! V-point 206 zi = pi - REAL(ii,wp) 207 zj = pj - 0.5 - REAL(ij,wp) 208 CASE ( 'F' ) 209 ii = MAX(1, INT( pi-0.5 )) 210 ij = MAX(1, INT( pj-0.5 )) ! F-point 211 zi = pi - 0.5 - REAL(ii,wp) 212 zj = pj - 0.5 - REAL(ij,wp) 213 END SELECT 214 ! 215 ! find position in this processor. Prevent near edge problems (see #1389) 216 ! 217 IF ( ii < mig( 1 ) ) THEN ; ii = 1 218 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 219 ELSE ; ii = mi1(ii) 220 ENDIF 221 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 222 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 223 ELSE ; ij = mj1(ij) 224 ENDIF 225 ! 226 IF( ii == jpi ) ii = ii-1 227 IF( ij == jpj ) ij = ij-1 228 ! 229 icb_utl_bilin_h = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) & 230 & + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) * zj 231 ! 232 END FUNCTION icb_utl_bilin_h 233 234 235 REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type ) 236 !!---------------------------------------------------------------------- 237 !! *** FUNCTION icb_utl_bilin *** 238 !! 239 !! ** Purpose : bilinear interpolation at berg location depending on the grid-point type 240 !! 241 !! !!gm CAUTION an optional argument should be added to handle 242 !! the slip/no-slip conditions ==>>> to be done later 243 !! 244 !!---------------------------------------------------------------------- 245 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pfld ! field to be interpolated 246 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 247 CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points 183 248 ! 184 249 INTEGER :: ii, ij ! local integer … … 213 278 ! 214 279 ! find position in this processor. Prevent near edge problems (see #1389) 215 216 if (ii.lt.mig(1)) then 217 ii = 1 218 else if (ii.gt.mig(jpi)) then 219 ii = jpi 220 else 221 ii = mi1( ii ) 222 end if 223 224 if (ij.lt.mjg(1)) then 225 ij = 1 226 else if (ij.gt.mjg(jpj)) then 227 ij = jpj 228 else 229 ij = mj1( ij ) 230 end if 231 232 if (ij.eq.jpj) ij=ij-1 233 if (ii.eq.jpi) ii=ii-1 234 235 ! 236 icb_utl_bilin_h = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) & 237 & + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) * zj 238 ! 239 END FUNCTION icb_utl_bilin_h 240 241 242 REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type ) 243 !!---------------------------------------------------------------------- 244 !! *** FUNCTION icb_utl_bilin *** 245 !! 246 !! ** Purpose : bilinear interpolation at berg location depending on the grid-point type 247 !! 248 !! !!gm CAUTION an optional argument should be added to handle 249 !! the slip/no-slip conditions ==>>> to be done later 250 !! 251 !!---------------------------------------------------------------------- 252 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pfld ! field to be interpolated 253 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 254 CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points 255 ! 256 INTEGER :: ii, ij ! local integer 257 REAL(wp) :: zi, zj ! local real 258 !!---------------------------------------------------------------------- 259 ! 260 SELECT CASE ( cd_type ) 261 CASE ( 'T' ) 262 ! note that here there is no +0.5 added 263 ! since we're looking for four T points containing quadrant we're in of 264 ! current T cell 265 ii = MAX(1, INT( pi )) 266 ij = MAX(1, INT( pj )) ! T-point 267 zi = pi - REAL(ii,wp) 268 zj = pj - REAL(ij,wp) 269 CASE ( 'U' ) 270 ii = MAX(1, INT( pi-0.5 )) 271 ij = MAX(1, INT( pj )) ! U-point 272 zi = pi - 0.5 - REAL(ii,wp) 273 zj = pj - REAL(ij,wp) 274 CASE ( 'V' ) 275 ii = MAX(1, INT( pi )) 276 ij = MAX(1, INT( pj-0.5 )) ! V-point 277 zi = pi - REAL(ii,wp) 278 zj = pj - 0.5 - REAL(ij,wp) 279 CASE ( 'F' ) 280 ii = MAX(1, INT( pi-0.5 )) 281 ij = MAX(1, INT( pj-0.5 )) ! F-point 282 zi = pi - 0.5 - REAL(ii,wp) 283 zj = pj - 0.5 - REAL(ij,wp) 284 END SELECT 285 ! 286 ! find position in this processor. Prevent near edge problems (see #1389) 287 288 if (ii.lt.mig(1)) then 289 ii = 1 290 else if (ii.gt.mig(jpi)) then 291 ii = jpi 292 else 293 ii = mi1( ii ) 294 end if 295 296 if (ij.lt.mjg(1)) then 297 ij = 1 298 else if (ij.gt.mjg(jpj)) then 299 ij = jpj 300 else 301 ij = mj1( ij ) 302 end if 303 304 if (ij.eq.jpj) ij=ij-1 305 if (ii.eq.jpi) ii=ii-1 306 280 IF ( ii < mig( 1 ) ) THEN ; ii = 1 281 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 282 ELSE ; ii = mi1(ii) 283 ENDIF 284 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 285 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 286 ELSE ; ij = mj1(ij) 287 ENDIF 288 ! 289 IF( ii == jpi ) ii = ii-1 290 IF( ij == jpj ) ij = ij-1 291 ! 307 292 icb_utl_bilin = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) & 308 293 & + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) * zj … … 340 325 ! 341 326 ! find position in this processor. Prevent near edge problems (see #1389) 342 343 if (ii.lt.mig(1)) then 344 ii = 1 345 else if (ii.gt.mig(jpi)) then 346 ii = jpi 347 else 348 ii = mi1( ii ) 349 end if 350 351 if (ij.lt.mjg(1)) then 352 ij = 1 353 else if (ij.gt.mjg(jpj)) then 354 ij = jpj 355 else 356 ij = mj1( ij ) 357 end if 358 359 if (ij.eq.jpj) ij=ij-1 360 if (ii.eq.jpi) ii=ii-1 361 327 IF ( ii < mig( 1 ) ) THEN ; ii = 1 328 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 329 ELSE ; ii = mi1(ii) 330 ENDIF 331 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 332 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 333 ELSE ; ij = mj1(ij) 334 ENDIF 335 ! 336 IF( ii == jpi ) ii = ii-1 337 IF( ij == jpj ) ij = ij-1 338 ! 362 339 z4(1) = pfld(ii ,ij ) 363 340 z4(2) = pfld(ii+1,ij ) … … 408 385 409 386 ! find position in this processor. Prevent near edge problems (see #1389) 410 411 if (ii.lt.mig(1)) then 412 ii = 1 413 else if (ii.gt.mig(jpi)) then 414 ii = jpi 415 else 416 ii = mi1( ii ) 417 end if 418 419 if (ij.lt.mjg(1)) then 420 ij = 1 421 else if (ij.gt.mjg(jpj)) then 422 ij = jpj 423 else 424 ij = mj1( ij ) 425 end if 426 427 if (ij.eq.jpj) ij=ij-1 428 if (ii.eq.jpi) ii=ii-1 429 387 IF ( ii < mig( 1 ) ) THEN ; ii = 1 388 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 389 ELSE ; ii = mi1(ii) 390 ENDIF 391 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 392 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 393 ELSE ; ij = mj1(ij) 394 ENDIF 395 ! 396 IF( ii == jpi ) ii = ii-1 397 IF( ij == jpj ) ij = ij-1 398 ! 430 399 IF( 0.0_wp <= zi .AND. zi < 0.5_wp ) THEN 431 400 IF( 0.0_wp <= zj .AND. zj < 0.5_wp ) THEN ! NE quadrant
Note: See TracChangeset
for help on using the changeset viewer.
