Changeset 13265 for NEMO/branches/2020
- Timestamp:
- 2020-07-08T13:34:28+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB/icb_oce.F90
r12472 r13265 86 86 ! particularly for MPP when iceberg can lie inside T grid but outside U, V, or f grid 87 87 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: uo_e, vo_e 88 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ff_e,tt_e, fr_e88 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: tt_e, fr_e 89 89 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ua_e, va_e 90 90 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ssh_e 91 91 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: tmask_e, umask_e, vmask_e 92 REAl(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: rlon_e, rlat_e, ff_e 93 ! prepration to Nacho work 94 ! tt3d_e 95 ! uo3d_e 96 ! vo3d_e 97 ! 92 98 #if defined key_si3 || defined key_cice 93 99 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hi_e, ui_e, vi_e … … 181 187 & hi_e(0:jpi+1,0:jpj+1) , & 182 188 #endif 183 & f f_e(0:jpi+1,0:jpj+1) , fr_e(0:jpi+1,0:jpj+1) ,&189 & fr_e(0:jpi+1,0:jpj+1) , & 184 190 & tt_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) , & 185 191 & first_width(nclasses) , first_length(nclasses) , & … … 188 194 icb_alloc = icb_alloc + ill 189 195 196 ! prepration to Nacho work 197 ! tt3d_e 198 ! uo3d_e 199 ! vo3d_e 200 ! 190 201 ALLOCATE( tmask_e(0:jpi+1,0:jpj+1), umask_e(0:jpi+1,0:jpj+1), vmask_e(0:jpi+1,0:jpj+1), & 191 & STAT=ill)202 & rlon_e(0:jpi+1,0:jpj+1) , rlat_e(0:jpi+1,0:jpj+1) , ff_e(0:jpi+1,0:jpj+1) , STAT=ill) 192 203 icb_alloc = icb_alloc + ill 193 204 -
NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB/icbdyn.F90
r10570 r13265 14 14 USE dom_oce ! NEMO ocean domain 15 15 USE phycst ! NEMO physical constants 16 USE in_out_manager ! IO parameters 16 17 ! 17 18 USE icb_oce ! define iceberg arrays … … 155 156 pt%xi = zxi_n 156 157 pt%yj = zyj_n 157 158 ! update actual position159 pt%lon = icb_utl_bilin_x(glamt, pt%xi, pt%yj )160 pt%lat = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' )161 158 162 159 berg => berg%next ! switch to the next berg … … 270 267 ! Interpolate gridded fields to berg 271 268 nknberg = berg%number(1) 272 CALL icb_utl_interp( pxi, pe1, zuo, zui, zua, zssh_x, & 273 & pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff ) 269 CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2, & ! scale factor 270 & puo=zuo, pui=zui, pua=zua, & ! oce/ice/atm velocities 271 & pvo=zvo, pvi=zvi, pva=zva, & ! oce/ice/atm velocities 272 & pssh_i=zssh_x, pssh_j=zssh_y, & ! ssh gradient 273 & phi=zhi, pff=zff ) ! ice thickness and coriolis 274 275 IF (lwp) THEN 276 WRITE(numout,*) 'icbdyn ', pxi, pyj 277 WRITE(numout,*) pe1, zuo, zui, zua, zssh_x 278 WRITE(numout,*) pe2, zvo, zvi, zva, zssh_y 279 WRITE(numout,*) zhi, zff 280 WRITE(numout,*) '' 281 END IF 282 274 283 275 284 zM = berg%current_point%mass -
NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB/icbini.F90
r12489 r13265 82 82 ff_e(:,:) = 0._wp ; tt_e(:,:) = 0._wp ; 83 83 fr_e(:,:) = 0._wp ; 84 ! prepration to Nacho work 85 ! e3t_e 86 ! tt3d_e 87 ! uo3d_e 88 ! vo3d_e 89 ! 84 90 #if defined key_si3 85 91 hi_e(:,:) = 0._wp ; … … 242 248 CALL lbc_lnk_icb( 'icbini', umask_e, 'T', +1._wp, 1, 1 ) 243 249 CALL lbc_lnk_icb( 'icbini', vmask_e, 'T', +1._wp, 1, 1 ) 244 ! 250 251 ! definition of extended lat/lon array needed by icb_bilin_h 252 rlon_e(:,:) = 0._wp ; rlon_e(1:jpi,1:jpj) = glamt(:,:) 253 rlat_e(:,:) = 0._wp ; rlat_e(1:jpi,1:jpj) = gphit(:,:) 254 CALL lbc_lnk_icb( 'icbini', rlon_e, 'T', +1._wp, 1, 1 ) 255 CALL lbc_lnk_icb( 'icbini', rlat_e, 'T', +1._wp, 1, 1 ) 256 ! 257 ! definnitionn of extennded ff_f array needed by icb_utl_interp 258 ff_e(:,:) = 0._wp ; ff_e(1:jpi,1:jpj) = ff_f(:,:) 259 CALL lbc_lnk_icb( 'icbini', ff_e, 'F', +1._wp, 1, 1 ) 260 245 261 ! assign each new iceberg with a unique number constructed from the processor number 246 262 ! and incremented by the total number of processors … … 338 354 localpt%xi = REAL( mig(ji), wp ) 339 355 localpt%yj = REAL( mjg(jj), wp ) 340 localpt%lon = icb_utl_bilin(glamt, localpt%xi, localpt%yj, 'T' ) 341 localpt%lat = icb_utl_bilin(gphit, localpt%xi, localpt%yj, 'T' ) 356 CALL icb_utl_interp( localpt%xi, localpt%yj, plat=localpt%lat, plon=localpt%lon ) 342 357 localpt%mass = rn_initial_mass (iberg) 343 358 localpt%thickness = rn_initial_thickness(iberg) -
NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB/icbthm.F90
r13226 r13265 83 83 pt => this%current_point 84 84 nknberg = this%number(1) 85 CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x, & 86 & pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y, & 87 & pt%sst, pt%cn, pt%hi, zff ) 85 IF ( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) THEN 86 CALL icb_utl_interp( pt%xi, pt%yj, pe1=pt%e1, pe2=pt%e2, & 87 & puo=pt%uo, pui=pt%ui, pua=pt%ua, pssh_i=pt%ssh_x, & 88 & pvo=pt%vo, pvi=pt%vi, pva=pt%va, pssh_j=pt%ssh_y, & 89 & psst=pt%sst, pcn=pt%cn, phi=pt%hi, & 90 & plat=pt%lat, plon=pt%lon ) 91 ELSE 92 CALL icb_utl_interp( pt%xi, pt%yj, & ! position 93 & puo=pt%uo, pua=pt%ua, & ! oce/atm velocities 94 & pvo=pt%vo, pva=pt%va, & ! oce/atm velocities 95 & psst=pt%sst, pcn=pt%cn ) ! sst and ice concentration 96 ! preparation Nacho add t3d and uo, vo, to basal 97 END IF 88 98 ! 89 99 zSST = pt%sst -
NEMO/branches/2020/tickets_2494_2375/src/OCE/ICB/icbutl.F90
r10702 r13265 8 8 !! - ! Removal of mapping from another grid 9 9 !! - ! 2011-04 (Alderson) Split into separate modules 10 !! 4.2 ! 2020-07 (P. Mathiot) simplification of interpolation routine 11 !! ! and add Nacho Merino work 10 12 !!---------------------------------------------------------------------- 11 13 12 14 !!---------------------------------------------------------------------- 13 15 !! icb_utl_interp : 14 !! icb_utl_bilin : 15 !! icb_utl_bilin_e : 16 !! icb_utl_pos : compute bottom left corner indice, weight and mask 17 !! icb_utl_bilin_h : interpolation field to icb position 18 !! icb_utl_bilin_e : interpolation of scale factor to icb position 16 19 !!---------------------------------------------------------------------- 17 20 USE par_oce ! ocean parameters … … 33 36 PUBLIC icb_utl_copy ! routine called in icbstp module 34 37 PUBLIC icb_utl_interp ! routine called in icbdyn, icbthm modules 35 PUBLIC icb_utl_bilin ! routine called in icbini, icbdyn modules 36 PUBLIC icb_utl_bilin_x ! routine called in icbdyn module 38 PUBLIC icb_utl_bilin_h ! routine called in icbdyn module 37 39 PUBLIC icb_utl_add ! routine called in icbini.F90, icbclv, icblbc and icbrst modules 38 40 PUBLIC icb_utl_delete ! routine called in icblbc, icbthm modules … … 72 74 uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 73 75 vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 74 ff_e(1:jpi,1:jpj) = ff_f (:,:)75 76 tt_e(1:jpi,1:jpj) = sst_m(:,:) 76 77 fr_e(1:jpi,1:jpj) = fr_i (:,:) … … 80 81 CALL lbc_lnk_icb( 'icbutl', uo_e, 'U', -1._wp, 1, 1 ) 81 82 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 83 CALL lbc_lnk_icb( 'icbutl', ua_e, 'U', -1._wp, 1, 1 ) 84 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 85 #if defined key_si3 88 86 hi_e(1:jpi, 1:jpj) = hm_i (:,:) … … 94 92 ssh_e(1:jpi, 1:jpj) = zssh_lead_m(:,:) * tmask(:,:,1) 95 93 ! 96 CALL lbc_lnk_icb( 'icbutl', hi_e , 'T', +1._wp, 1, 1 )97 94 CALL lbc_lnk_icb( 'icbutl', ui_e , 'U', -1._wp, 1, 1 ) 98 95 CALL lbc_lnk_icb( 'icbutl', vi_e , 'V', -1._wp, 1, 1 ) 99 96 #else 100 ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1) 97 ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1) 101 98 #endif 102 CALL lbc_lnk_icb( 'icbutl', ssh_e, 'T', +1._wp, 1, 1 )103 99 ! 104 100 END SUBROUTINE icb_utl_copy 105 101 106 102 107 SUBROUTINE icb_utl_interp( pi, p e1, puo, pui, pua, pssh_i, &108 & pj,pe2, pvo, pvi, pva, pssh_j, &109 & psst, pcn, phi, pff)103 SUBROUTINE icb_utl_interp( pi, pj, pe1, puo, pui, pua, pssh_i, & 104 & pe2, pvo, pvi, pva, pssh_j, & 105 & psst, pcn, phi, pff, plon, plat) 110 106 !!---------------------------------------------------------------------- 111 107 !! *** ROUTINE icb_utl_interp *** … … 125 121 !!---------------------------------------------------------------------- 126 122 REAL(wp), INTENT(in ) :: pi , pj ! position in (i,j) referential 127 REAL(wp), INTENT( out) :: pe1, pe2 ! i- and j scale factors 128 REAL(wp), INTENT( out) :: puo, pvo, pui, pvi, pua, pva ! ocean, ice and wind speeds 129 REAL(wp), INTENT( out) :: pssh_i, pssh_j ! ssh i- & j-gradients 130 REAL(wp), INTENT( out) :: psst, pcn, phi, pff ! SST, ice concentration, ice thickness, Coriolis 131 ! 123 REAL(wp), INTENT( out), OPTIONAL :: pe1, pe2 ! i- and j scale factors 124 REAL(wp), INTENT( out), OPTIONAL :: puo, pvo, pui, pvi, pua, pva ! ocean, ice and wind speeds 125 REAL(wp), INTENT( out), OPTIONAL :: pssh_i, pssh_j ! ssh i- & j-gradients 126 REAL(wp), INTENT( out), OPTIONAL :: psst, pcn, phi, pff ! SST, ice concentration, ice thickness, Coriolis 127 REAL(wp), INTENT( out), OPTIONAL :: plat, plon ! position 128 ! 129 REAL(wp), DIMENSION(4) :: zwT , zwU , zwV , zwF ! interpolation weight 130 REAL(wp), DIMENSION(4) :: zmskF, zmskU, zmskV, zmskT ! mask 131 REAL(wp), DIMENSION(4) :: zwTp, zmskTp, zwTm, zmskTm 132 INTEGER :: iiT, iiU, iiV, iiF, ijT, ijU, ijV, ijF ! bottom left corner 133 INTEGER :: iiTp, iiTm, ijTp, ijTm 132 134 REAL(wp) :: zcd, zmod ! local scalars 133 135 !!---------------------------------------------------------------------- 134 135 pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj ) ! scale factors 136 pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 137 ! 138 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 147 zmod = 1._wp / MAX( 1.e-20, SQRT( zcd * SQRT( pua*pua + pva*pva) ) ) 148 pua = pua * zmod ! note: stress module=0 necessarly implies ua=va=0 149 pva = pva * zmod 150 136 ! 137 ! get position, weight and mask 138 CALL icb_utl_pos( pi, pj, 'T', iiT, ijT, zwT, zmskT ) 139 CALL icb_utl_pos( pi, pj, 'U', iiU, ijU, zwU, zmskU ) 140 CALL icb_utl_pos( pi, pj, 'V', iiV, ijV, zwV, zmskV ) 141 CALL icb_utl_pos( pi, pj, 'F', iiF, ijF, zwF, zmskF ) 142 ! 143 ! metrics and coordinates 144 IF ( PRESENT(pe1 ) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj ) ! scale factors 145 IF ( PRESENT(pe2 ) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 146 IF ( PRESENT(plon) ) plon= icb_utl_bilin_h( rlon_e, iiT, ijT, zwT, .true. ) 147 IF ( PRESENT(plat) ) plat= icb_utl_bilin_h( rlat_e, iiT, ijT, zwT, .false. ) 148 ! 149 IF ( PRESENT(puo ) ) puo = icb_utl_bilin_h( uo_e, iiU, ijU, zwU , .false. ) ! ocean velocities 150 IF ( PRESENT(pvo ) ) pvo = icb_utl_bilin_h( vo_e, iiV, ijV, zwV , .false. ) ! 151 IF ( PRESENT(psst) ) psst = icb_utl_bilin_h( tt_e, iiT, ijT, zwT * zmskT, .false. ) ! sst 152 IF ( PRESENT(pcn ) ) pcn = icb_utl_bilin_h( fr_e, iiT, ijT, zwT * zmskT, .false. ) ! ice concentration 153 IF ( PRESENT(pff ) ) pff = icb_utl_bilin_h( ff_e, iiF, ijF, zwF , .false. ) ! Coriolis parameter 154 ! 155 IF ( PRESENT(pua) .AND. PRESENT(pva) ) THEN 156 pua = icb_utl_bilin_h( ua_e, iiU, ijU, zwU * zmskU, .false. ) ! 10m wind 157 pva = icb_utl_bilin_h( va_e, iiV, ijV, zwV * zmskV, .false. ) ! here (ua,va) are stress => rough conversion from stress to speed 158 zcd = 1.22_wp * 1.5e-3_wp ! air density * drag coefficient 159 zmod = 1._wp / MAX( 1.e-20, SQRT( zcd * SQRT( pua*pua + pva*pva) ) ) 160 pua = pua * zmod ! note: stress module=0 necessarly implies ua=va=0 161 pva = pva * zmod 162 END IF 163 ! 151 164 #if defined key_si3 152 pui = icb_utl_bilin_h( ui_e , pi, pj, 'U', .false. )! sea-ice velocities153 pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V', .false. )154 phi = icb_utl_bilin_h( hi_e , pi, pj, 'T', .true. )! ice thickness165 IF ( PRESENT(pui) ) pui = icb_utl_bilin_h( ui_e , iiU, ijU, zwU , .false. ) ! sea-ice velocities 166 IF ( PRESENT(pvi) ) pvi = icb_utl_bilin_h( vi_e , iiV, ijV, zwV , .false. ) 167 IF ( PRESENT(phi) ) phi = icb_utl_bilin_h( hi_e , iiT, ijT, zwT * zmskT, .false. ) ! ice thickness 155 168 #else 156 pui = 0._wp157 pvi = 0._wp158 phi = 0._wp169 IF ( PRESENT(pui) ) pui = 0._wp 170 IF ( PRESENT(pvi) ) pvi = 0._wp 171 IF ( PRESENT(phi) ) phi = 0._wp 159 172 #endif 160 173 ! 161 174 ! Estimate SSH gradient in i- and j-direction (centred evaluation) 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 IF ( PRESENT(pssh_i) .AND. PRESENT(pssh_j) ) THEN 176 CALL icb_utl_pos( pi+0.1, pj , 'T', iiTp, ijTp, zwTp, zmskTp ) 177 CALL icb_utl_pos( pi-0.1, pj , 'T', iiTm, ijTm, zwTm, zmskTm ) 178 ! 179 IF ( .NOT. PRESENT(pe1) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj ) 180 pssh_i = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) - & 181 & icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. ) ) / ( 0.2_wp * pe1 ) 182 ! 183 CALL icb_utl_pos( pi , pj+0.1, 'T', iiTp, ijTp, zwTp, zmskTp ) 184 CALL icb_utl_pos( pi , pj-0.1, 'T', iiTm, ijTm, zwTm, zmskTm ) 185 ! 186 IF ( .NOT. PRESENT(pe2) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 187 pssh_j = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) - & 188 & icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. ) ) / ( 0.2_wp * pe2 ) 189 END IF 166 190 ! 167 191 END SUBROUTINE icb_utl_interp 168 192 169 170 REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type, plmask ) 193 SUBROUTINE icb_utl_pos( pi, pj, cd_type, kii, kij, pw, pmsk ) 171 194 !!---------------------------------------------------------------------- 172 195 !! *** FUNCTION icb_utl_bilin *** … … 179 202 !! 180 203 !!---------------------------------------------------------------------- 181 REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) :: pfld ! field to be interpolated 182 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 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 185 ! 186 INTEGER :: ii, ij ! local integer 187 REAL(wp) :: zi, zj ! local real 188 REAL(wp) :: zw1, zw2, zw3, zw4 189 REAL(wp), DIMENSION(4) :: zmask 204 REAL(wp) , INTENT(IN) :: pi, pj ! targeted coordinates in (i,j) referential 205 CHARACTER(len=1) , INTENT(IN) :: cd_type ! point type 206 REAL(wp), DIMENSION(4), INTENT(OUT) :: pw, pmsk ! weight and mask 207 INTEGER , INTENT(OUT) :: kii, kij ! bottom left corner position in local domain 208 ! 209 REAL(wp) :: zwi, zwj ! distance to bottom left corner 210 INTEGER :: ierr 211 ! 190 212 !!---------------------------------------------------------------------- 191 213 ! … … 195 217 ! since we're looking for four T points containing quadrant we're in of 196 218 ! current T cell 197 ii = MAX(0, INT( pi ))198 ij = MAX(0, INT( pj )) ! T-point199 z i = pi - REAL(ii,wp)200 z j = pj - REAL(ij,wp)219 kii = MAX(0, INT( pi )) 220 kij = MAX(0, INT( pj )) ! T-point 221 zwi = pi - REAL(kii,wp) 222 zwj = pj - REAL(kij,wp) 201 223 CASE ( 'U' ) 202 ii = MAX(0, INT( pi-0.5_wp ))203 ij = MAX(0, INT( pj )) ! U-point204 z i = pi - 0.5_wp - REAL(ii,wp)205 z j = pj - REAL(ij,wp)224 kii = MAX(0, INT( pi-0.5_wp )) 225 kij = MAX(0, INT( pj )) ! U-point 226 zwi = pi - 0.5_wp - REAL(kii,wp) 227 zwj = pj - REAL(kij,wp) 206 228 CASE ( 'V' ) 207 ii = MAX(0, INT( pi ))208 ij = MAX(0, INT( pj-0.5_wp )) ! V-point209 z i = pi - REAL(ii,wp)210 z j = pj - 0.5_wp - REAL(ij,wp)229 kii = MAX(0, INT( pi )) 230 kij = MAX(0, INT( pj-0.5_wp )) ! V-point 231 zwi = pi - REAL(kii,wp) 232 zwj = pj - 0.5_wp - REAL(kij,wp) 211 233 CASE ( 'F' ) 212 ii = MAX(0, INT( pi-0.5_wp ))213 ij = MAX(0, INT( pj-0.5_wp )) ! F-point214 z i = pi - 0.5_wp - REAL(ii,wp)215 z j = pj - 0.5_wp - REAL(ij,wp)234 kii = MAX(0, INT( pi-0.5_wp )) 235 kij = MAX(0, INT( pj-0.5_wp )) ! F-point 236 zwi = pi - 0.5_wp - REAL(kii,wp) 237 zwj = pj - 0.5_wp - REAL(kij,wp) 216 238 END SELECT 239 ! 240 ! compute weight 241 pw(1) = (1._wp-zwi) * (1._wp-zwj) 242 pw(2) = zwi * (1._wp-zwj) 243 pw(3) = (1._wp-zwi) * zwj 244 pw(4) = zwi * zwj 245 ! 246 ! find position in this processor. Prevent near edge problems (see #1389) 247 ! 248 IF (TRIM(cd_type) == 'T' ) THEN 249 ierr = 0 250 IF ( kii < mig( 1 ) ) THEN ; ierr = ierr + 1 251 ELSEIF( kii >= mig(jpi) ) THEN ; ierr = ierr + 1 252 ENDIF 253 ! 254 IF ( kij < mjg( 1 ) ) THEN ; ierr = ierr + 1 255 ELSEIF( kij >= mjg(jpj) ) THEN ; ierr = ierr + 1 256 ENDIF 257 ! 258 IF ( ierr > 0 ) THEN 259 WRITE(numout,*) 'bottom left corner T point out of bound' 260 WRITE(numout,*) kii, mig( 1 ), mig(jpi) 261 WRITE(numout,*) kij, mjg( 1 ), mjg(jpj) 262 CALL ctl_stop('STOP','icb_utl_bilin_h: an icebergs coordinates is out of valid range (out of bound error)') 263 END IF 264 END IF 217 265 ! 218 266 ! find position in this processor. Prevent near edge problems (see #1389) 219 267 ! (PM) will be useless if extra halo is used in NEMO 220 268 ! 221 IF ( ii <= mig(1)-1 ) THEN ;ii = 0222 ELSEIF( ii > mig(jpi) ) THEN ;ii = jpi223 ELSE ; ii = mi1(ii)224 ENDIF 225 IF ( ij <= mjg(1)-1 ) THEN ;ij = 0226 ELSEIF( ij > mjg(jpj) ) THEN ;ij = jpj227 ELSE ; ij = mj1(ij)269 IF ( kii <= mig(1)-1 ) THEN ; kii = 0 270 ELSEIF( kii > mig(jpi) ) THEN ; kii = jpi 271 ELSE ; kii = mi1(kii) 272 ENDIF 273 IF ( kij <= mjg(1)-1 ) THEN ; kij = 0 274 ELSEIF( kij > mjg(jpj) ) THEN ; kij = jpj 275 ELSE ; kij = mj1(kij) 228 276 ENDIF 229 277 ! 230 278 ! 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) 257 ! 258 END FUNCTION icb_utl_bilin_h 259 260 261 REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type ) 279 ! land value is not used in the interpolation 280 SELECT CASE ( cd_type ) 281 CASE ( 'T' ) 282 pmsk = (/tmask_e(kii,kij), tmask_e(kii+1,kij), tmask_e(kii,kij+1), tmask_e(kii+1,kij+1)/) 283 CASE ( 'U' ) 284 pmsk = (/umask_e(kii,kij), umask_e(kii+1,kij), umask_e(kii,kij+1), umask_e(kii+1,kij+1)/) 285 CASE ( 'V' ) 286 pmsk = (/vmask_e(kii,kij), vmask_e(kii+1,kij), vmask_e(kii,kij+1), vmask_e(kii+1,kij+1)/) 287 CASE ( 'F' ) 288 ! F case only used for coriolis, ff_f is not mask so zmask = 1 289 pmsk = 1. 290 END SELECT 291 END SUBROUTINE icb_utl_pos 292 293 REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pii, pij, pw, pllon ) 262 294 !!---------------------------------------------------------------------- 263 295 !! *** FUNCTION icb_utl_bilin *** 264 296 !! 265 297 !! ** Purpose : bilinear interpolation at berg location depending on the grid-point type 298 !! this version deals with extra halo points 266 299 !! 267 300 !! !!gm CAUTION an optional argument should be added to handle … … 269 302 !! 270 303 !!---------------------------------------------------------------------- 271 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pfld ! field to be interpolated 272 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 273 CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points 274 ! 275 INTEGER :: ii, ij ! local integer 276 REAL(wp) :: zi, zj ! local real 277 !!---------------------------------------------------------------------- 278 ! 279 SELECT CASE ( cd_type ) 280 CASE ( 'T' ) 281 ! note that here there is no +0.5 added 282 ! since we're looking for four T points containing quadrant we're in of 283 ! current T cell 284 ii = MAX(1, INT( pi )) 285 ij = MAX(1, INT( pj )) ! T-point 286 zi = pi - REAL(ii,wp) 287 zj = pj - REAL(ij,wp) 288 CASE ( 'U' ) 289 ii = MAX(1, INT( pi-0.5 )) 290 ij = MAX(1, INT( pj )) ! U-point 291 zi = pi - 0.5 - REAL(ii,wp) 292 zj = pj - REAL(ij,wp) 293 CASE ( 'V' ) 294 ii = MAX(1, INT( pi )) 295 ij = MAX(1, INT( pj-0.5 )) ! V-point 296 zi = pi - REAL(ii,wp) 297 zj = pj - 0.5 - REAL(ij,wp) 298 CASE ( 'F' ) 299 ii = MAX(1, INT( pi-0.5 )) 300 ij = MAX(1, INT( pj-0.5 )) ! F-point 301 zi = pi - 0.5 - REAL(ii,wp) 302 zj = pj - 0.5 - REAL(ij,wp) 303 END SELECT 304 ! 305 ! find position in this processor. Prevent near edge problems (see #1389) 306 IF ( ii < mig( 1 ) ) THEN ; ii = 1 307 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 308 ELSE ; ii = mi1(ii) 309 ENDIF 310 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 311 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 312 ELSE ; ij = mj1(ij) 313 ENDIF 314 ! 315 IF( ii == jpi ) ii = ii-1 316 IF( ij == jpj ) ij = ij-1 317 ! 318 icb_utl_bilin = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) & 319 & + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) * zj 320 ! 321 END FUNCTION icb_utl_bilin 322 323 324 REAL(wp) FUNCTION icb_utl_bilin_x( pfld, pi, pj ) 325 !!---------------------------------------------------------------------- 326 !! *** FUNCTION icb_utl_bilin_x *** 327 !! 328 !! ** Purpose : bilinear interpolation at berg location depending on the grid-point type 329 !! Special case for interpolating longitude 330 !! 331 !! !!gm CAUTION an optional argument should be added to handle 332 !! the slip/no-slip conditions ==>>> to be done later 333 !! 334 !!---------------------------------------------------------------------- 335 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pfld ! field to be interpolated 336 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 337 ! 338 INTEGER :: ii, ij ! local integer 339 REAL(wp) :: zi, zj ! local real 340 REAL(wp) :: zret ! local real 341 REAL(wp), DIMENSION(4) :: z4 342 !!---------------------------------------------------------------------- 343 ! 344 ! note that here there is no +0.5 added 345 ! since we're looking for four T points containing quadrant we're in of 346 ! current T cell 347 ii = MAX(1, INT( pi )) 348 ij = MAX(1, INT( pj )) ! T-point 349 zi = pi - REAL(ii,wp) 350 zj = pj - REAL(ij,wp) 351 ! 352 ! find position in this processor. Prevent near edge problems (see #1389) 353 IF ( ii < mig( 1 ) ) THEN ; ii = 1 354 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 355 ELSE ; ii = mi1(ii) 356 ENDIF 357 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 358 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 359 ELSE ; ij = mj1(ij) 360 ENDIF 361 ! 362 IF( ii == jpi ) ii = ii-1 363 IF( ij == jpj ) ij = ij-1 364 ! 365 z4(1) = pfld(ii ,ij ) 366 z4(2) = pfld(ii+1,ij ) 367 z4(3) = pfld(ii ,ij+1) 368 z4(4) = pfld(ii+1,ij+1) 369 IF( MAXVAL(z4) - MINVAL(z4) > 90._wp ) THEN 370 WHERE( z4 < 0._wp ) z4 = z4 + 360._wp 371 ENDIF 372 ! 373 zret = (z4(1) * (1.-zi) + z4(2) * zi) * (1.-zj) + (z4(3) * (1.-zi) + z4(4) * zi) * zj 374 IF( zret > 180._wp ) zret = zret - 360._wp 375 icb_utl_bilin_x = zret 376 ! 377 END FUNCTION icb_utl_bilin_x 378 304 REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) :: pfld ! field to be interpolated 305 REAL(wp), DIMENSION(4) , INTENT(in) :: pw ! weight 306 LOGICAL , INTENT(in) :: pllon ! input data is a longitude 307 INTEGER , INTENT(in) :: pii, pij ! bottom left corner 308 ! 309 REAL(wp), DIMENSION(4) :: zdat ! input data 310 !!---------------------------------------------------------------------- 311 ! 312 ! data 313 zdat(1) = pfld(pii ,pij ) 314 zdat(2) = pfld(pii+1,pij ) 315 zdat(3) = pfld(pii ,pij+1) 316 zdat(4) = pfld(pii+1,pij+1) 317 ! 318 IF( pllon .AND. MAXVAL(zdat) - MINVAL(zdat) > 90._wp ) THEN 319 WHERE( zdat < 0._wp ) zdat = zdat + 360._wp 320 ENDIF 321 ! 322 ! compute interpolated value 323 icb_utl_bilin_h = ( zdat(1)*pw(1) + zdat(2)*pw(2) + zdat(3)*pw(3) + zdat(4)*pw(4) ) / MAX(1.e-20, pw(1)+pw(2)+pw(3)+pw(4)) 324 ! 325 IF( pllon .AND. icb_utl_bilin_h > 180._wp ) icb_utl_bilin_h = icb_utl_bilin_h - 360._wp 326 ! 327 END FUNCTION icb_utl_bilin_h 379 328 380 329 REAL(wp) FUNCTION icb_utl_bilin_e( pet, peu, pev, pef, pi, pj ) … … 387 336 !!---------------------------------------------------------------------- 388 337 REAL(wp), DIMENSION(:,:), INTENT(in) :: pet, peu, pev, pef ! horizontal scale factor to be interpolated at t-,u-,v- & f-pts 389 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 390 ! 391 INTEGER :: ii, ij, icase, ierr ! local integer 338 REAL(wp) , INTENT(IN) :: pi , pj ! iceberg position 392 339 ! 393 340 ! weights corresponding to corner points of a T cell quadrant 394 341 REAL(wp) :: zi, zj ! local real 342 INTEGER :: ii, ij ! bottom left corner coordinate in local domain 395 343 ! 396 344 ! values at corner points of a T cell quadrant … … 399 347 !!---------------------------------------------------------------------- 400 348 ! 349 ! cannot used iiT because need ii/ij reltaive to global indices not local one 401 350 ii = MAX(1, INT( pi )) ; ij = MAX(1, INT( pj )) ! left bottom T-point (i,j) indices 402 351 ! 403 352 ! fractional box spacing 404 353 ! 0 <= zi < 0.5, 0 <= zj < 0.5 --> NW quadrant of current T cell … … 410 359 zj = pj - REAL(ij,wp) 411 360 412 ! find position in this processor. Prevent near edge problems (see #1389) 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 417 ELSE ; ii = mi1(ii) 418 ENDIF 419 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 ; ierr = ierr + 1 420 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj ; ierr = ierr + 1 421 ELSE ; ij = mj1(ij) 422 ENDIF 423 ! 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)') 361 ! conversion to local domain (no need to do a sanity check already done in icbpos) 362 ii = mi1(ii) 363 ij = mj1(ij) 364 365 IF (lwp) THEN 366 WRITE(numout,*) 'bilin_e ', pi, pj 367 WRITE(numout,*) ii, ij, zi, zj 368 WRITE(numout,*) '' 369 END IF 370 428 371 ! 429 372 IF( 0.0_wp <= zi .AND. zi < 0.5_wp ) THEN
Note: See TracChangeset
for help on using the changeset viewer.