Changeset 10726 for NEMO/releases/release-4.0/src/OCE/ICB/icbutl.F90
- Timestamp:
- 2019-02-27T16:06:35+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/releases/release-4.0/src/OCE/ICB/icbutl.F90
r10570 r10726 70 70 ! and ssh which is used to calculate gradients 71 71 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 sbcblk72 uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 73 vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 74 ff_e(1:jpi,1:jpj) = ff_f (:,:) 75 tt_e(1:jpi,1:jpj) = sst_m(:,:) 76 fr_e(1:jpi,1:jpj) = fr_i (:,:) 77 ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk 78 va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk 79 79 ! 80 80 CALL lbc_lnk_icb( 'icbutl', uo_e, 'U', -1._wp, 1, 1 ) … … 86 86 CALL lbc_lnk_icb( 'icbutl', tt_e, 'T', +1._wp, 1, 1 ) 87 87 #if defined key_si3 88 hi cth(:,:) = 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(:,:)88 hi_e(1:jpi, 1:jpj) = hm_i (:,:) 89 ui_e(1:jpi, 1:jpj) = u_ice(:,:) 90 vi_e(1:jpi, 1:jpj) = v_ice(:,:) 91 91 ! 92 92 ! compute ssh slope using ssh_lead if embedded 93 93 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', hi cth, 'T', +1._wp, 1, 1 )94 ssh_e(1:jpi, 1:jpj) = zssh_lead_m(:,:) * tmask(:,:,1) 95 ! 96 CALL lbc_lnk_icb( 'icbutl', hi_e , 'T', +1._wp, 1, 1 ) 97 97 CALL lbc_lnk_icb( 'icbutl', ui_e , 'U', -1._wp, 1, 1 ) 98 98 CALL lbc_lnk_icb( 'icbutl', vi_e , 'V', -1._wp, 1, 1 ) 99 99 #else 100 ssh_e( :,:) = 0._wp ; ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1)100 ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1) 101 101 #endif 102 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 102 CALL lbc_lnk_icb( 'icbutl', ssh_e, 'T', +1._wp, 1, 1 ) 114 103 ! … … 131 120 !! is half the off shore value, wile the normal-to-the-coast value is zero. 132 121 !! This is OK as a starting point. 122 !! !!pm HARD CODED: - rho_air now computed in sbcblk (what are the effect ?) 123 !! - drag coefficient (should it be namelist parameter ?) 133 124 !! 134 125 !!---------------------------------------------------------------------- … … 142 133 !!---------------------------------------------------------------------- 143 134 144 pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj ) ! scale factors135 pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj ) ! scale factors 145 136 pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 146 137 ! 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 speed155 zcd = 1.22_wp * 1.5e-3_wp ! air density * drag coefficient138 puo = icb_utl_bilin_h( uo_e, pi, pj, 'U', .false. ) ! ocean velocities 139 pvo = icb_utl_bilin_h( vo_e, pi, pj, 'V', .false. ) 140 psst = icb_utl_bilin_h( tt_e, pi, pj, 'T', .true. ) ! SST 141 pcn = icb_utl_bilin_h( fr_e, pi, pj, 'T', .true. ) ! ice concentration 142 pff = icb_utl_bilin_h( ff_e, pi, pj, 'F', .false. ) ! Coriolis parameter 143 ! 144 pua = icb_utl_bilin_h( ua_e, pi, pj, 'U', .true. ) ! 10m wind 145 pva = icb_utl_bilin_h( va_e, pi, pj, 'V', .true. ) ! here (ua,va) are stress => rough conversion from stress to speed 146 zcd = 1.22_wp * 1.5e-3_wp ! air density * drag coefficient 156 147 zmod = 1._wp / MAX( 1.e-20, SQRT( zcd * SQRT( pua*pua + pva*pva) ) ) 157 148 pua = pua * zmod ! note: stress module=0 necessarly implies ua=va=0 … … 159 150 160 151 #if defined key_si3 161 pui = icb_utl_bilin_h( ui_e , pi, pj, 'U' )! sea-ice velocities162 pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V' )163 phi = icb_utl_bilin_h( hi cth, pi, pj, 'T' )! ice thickness152 pui = icb_utl_bilin_h( ui_e , pi, pj, 'U', .false. ) ! sea-ice velocities 153 pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V', .false. ) 154 phi = icb_utl_bilin_h( hi_e , pi, pj, 'T', .true. ) ! ice thickness 164 155 #else 165 156 pui = 0._wp … … 169 160 170 161 ! 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 )162 pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T', .true. ) - & 163 & icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T', .true. ) ) / ( 0.2_wp * pe1 ) 164 pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T', .true. ) - & 165 & icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T', .true. ) ) / ( 0.2_wp * pe2 ) 175 166 ! 176 167 END SUBROUTINE icb_utl_interp 177 168 178 169 179 REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type )170 REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type, plmask ) 180 171 !!---------------------------------------------------------------------- 181 172 !! *** FUNCTION icb_utl_bilin *** … … 191 182 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 192 183 CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points 184 LOGICAL , INTENT(in) :: plmask ! special treatment of mask point 193 185 ! 194 186 INTEGER :: ii, ij ! local integer 195 187 REAL(wp) :: zi, zj ! local real 188 REAL(wp) :: zw1, zw2, zw3, zw4 189 REAL(wp), DIMENSION(4) :: zmask 196 190 !!---------------------------------------------------------------------- 197 191 ! … … 201 195 ! since we're looking for four T points containing quadrant we're in of 202 196 ! current T cell 203 ii = MAX( 1, INT( pi ))204 ij = MAX( 1, INT( pj )) ! T-point197 ii = MAX(0, INT( pi )) 198 ij = MAX(0, INT( pj )) ! T-point 205 199 zi = pi - REAL(ii,wp) 206 200 zj = pj - REAL(ij,wp) 207 201 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)202 ii = MAX(0, INT( pi-0.5_wp )) 203 ij = MAX(0, INT( pj )) ! U-point 204 zi = pi - 0.5_wp - REAL(ii,wp) 211 205 zj = pj - REAL(ij,wp) 212 206 CASE ( 'V' ) 213 ii = MAX( 1, INT( pi ))214 ij = MAX( 1, INT( pj-0.5)) ! V-point207 ii = MAX(0, INT( pi )) 208 ij = MAX(0, INT( pj-0.5_wp )) ! V-point 215 209 zi = pi - REAL(ii,wp) 216 zj = pj - 0.5 - REAL(ij,wp)210 zj = pj - 0.5_wp - REAL(ij,wp) 217 211 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)212 ii = MAX(0, INT( pi-0.5_wp )) 213 ij = MAX(0, INT( pj-0.5_wp )) ! F-point 214 zi = pi - 0.5_wp - REAL(ii,wp) 215 zj = pj - 0.5_wp - REAL(ij,wp) 222 216 END SELECT 223 217 ! 224 218 ! find position in this processor. Prevent near edge problems (see #1389) 225 ! 226 IF ( ii < mig( 1 ) ) THEN ; ii = 1 227 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 228 ELSE ; ii = mi1(ii) 229 ENDIF 230 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 231 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 232 ELSE ; ij = mj1(ij) 233 ENDIF 234 ! 235 IF( ii == jpi ) ii = ii-1 236 IF( ij == jpj ) ij = ij-1 237 ! 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 ) * zj 219 ! (PM) will be useless if extra halo is used in NEMO 220 ! 221 IF ( ii <= mig(1)-1 ) THEN ; ii = 0 222 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 223 ELSE ; ii = mi1(ii) 224 ENDIF 225 IF ( ij <= mjg(1)-1 ) THEN ; ij = 0 226 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 227 ELSE ; ij = mj1(ij) 228 ENDIF 229 ! 230 ! define mask array 231 IF (plmask) THEN 232 ! land value is not used in the interpolation 233 SELECT CASE ( cd_type ) 234 CASE ( 'T' ) 235 zmask = (/tmask_e(ii,ij), tmask_e(ii+1,ij), tmask_e(ii,ij+1), tmask_e(ii+1,ij+1)/) 236 CASE ( 'U' ) 237 zmask = (/umask_e(ii,ij), umask_e(ii+1,ij), umask_e(ii,ij+1), umask_e(ii+1,ij+1)/) 238 CASE ( 'V' ) 239 zmask = (/vmask_e(ii,ij), vmask_e(ii+1,ij), vmask_e(ii,ij+1), vmask_e(ii+1,ij+1)/) 240 CASE ( 'F' ) 241 ! F case only used for coriolis, ff_f is not mask so zmask = 1 242 zmask = 1. 243 END SELECT 244 ELSE 245 ! land value is used during interpolation 246 zmask = 1. 247 END iF 248 ! 249 ! compute weight 250 zw1 = zmask(1) * (1._wp-zi) * (1._wp-zj) 251 zw2 = zmask(2) * zi * (1._wp-zj) 252 zw3 = zmask(3) * (1._wp-zi) * zj 253 zw4 = zmask(4) * zi * zj 254 ! 255 ! compute interpolated value 256 icb_utl_bilin_h = ( pfld(ii,ij)*zw1 + pfld(ii+1,ij)*zw2 + pfld(ii,ij+1)*zw3 + pfld(ii+1,ij+1)*zw4 ) / MAX(1.e-20, zw1+zw2+zw3+zw4) 240 257 ! 241 258 END FUNCTION icb_utl_bilin_h … … 372 389 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 373 390 ! 374 INTEGER :: ii, ij, icase ! local integer391 INTEGER :: ii, ij, icase, ierr ! local integer 375 392 ! 376 393 ! weights corresponding to corner points of a T cell quadrant … … 394 411 395 412 ! 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 413 ! 414 ierr = 0 415 IF ( ii < mig( 1 ) ) THEN ; ii = 1 ; ierr = ierr + 1 416 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi ; ierr = ierr + 1 398 417 ELSE ; ii = mi1(ii) 399 418 ENDIF 400 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 401 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 419 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 ; ierr = ierr + 1 420 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj ; ierr = ierr + 1 402 421 ELSE ; ij = mj1(ij) 403 422 ENDIF 404 423 ! 405 IF( ii == jpi ) ii = ii-1 406 IF( ij == jpj ) ij = ij-1 424 IF( ii == jpi ) THEN ; ii = ii-1 ; ierr = ierr + 1 ; END IF 425 IF( ij == jpj ) THEN ; ij = ij-1 ; ierr = ierr + 1 ; END IF 426 ! 427 IF ( ierr > 0 ) CALL ctl_stop('STOP','icb_utl_bilin_e: an icebergs coordinates is out of valid range (out of bound error)') 407 428 ! 408 429 IF( 0.0_wp <= zi .AND. zi < 0.5_wp ) THEN … … 436 457 ENDIF 437 458 ! 438 icb_utl_bilin_e = ( ze01 * (1. -zi) + ze11 * zi ) *zj &439 & + ( ze00 * (1. -zi) + ze10 * zi ) * (1.-zj)459 icb_utl_bilin_e = ( ze01 * (1._wp-zi) + ze11 * zi ) * zj & 460 & + ( ze00 * (1._wp-zi) + ze10 * zi ) * (1._wp-zj) 440 461 ! 441 462 END FUNCTION icb_utl_bilin_e
Note: See TracChangeset
for help on using the changeset viewer.