Changeset 14038 for NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ICB/icbutl.F90
- Timestamp:
- 2020-12-03T12:25:19+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/ICB/icbutl.F90
r13281 r14038 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 21 USE oce, ONLY: ts, uu, vv 18 22 USE dom_oce ! ocean domain 19 23 USE in_out_manager ! IO parameters … … 31 35 PRIVATE 32 36 37 INTERFACE icb_utl_bilin_h 38 MODULE PROCEDURE icb_utl_bilin_2d_h, icb_utl_bilin_3d_h 39 END INTERFACE 40 33 41 PUBLIC icb_utl_copy ! routine called in icbstp module 42 PUBLIC icb_utl_getkb ! routine called in icbdyn and icbthm modules 43 PUBLIC test_icb_utl_getkb ! routine called in icbdyn and icbthm modules 44 PUBLIC icb_utl_zavg ! routine called in icbdyn and icbthm modules 34 45 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 46 PUBLIC icb_utl_bilin_h ! routine called in icbdyn module 37 47 PUBLIC icb_utl_add ! routine called in icbini.F90, icbclv, icblbc and icbrst modules 38 48 PUBLIC icb_utl_delete ! routine called in icblbc, icbthm modules … … 54 64 CONTAINS 55 65 56 SUBROUTINE icb_utl_copy( )66 SUBROUTINE icb_utl_copy( Kmm ) 57 67 !!---------------------------------------------------------------------- 58 68 !! *** ROUTINE icb_utl_copy *** … … 62 72 !! ** Method : - blah blah 63 73 !!---------------------------------------------------------------------- 74 REAL(wp), DIMENSION(0:jpi+1,0:jpj+1) :: ztmp 64 75 #if defined key_si3 65 76 REAL(wp), DIMENSION(jpi,jpj) :: zssh_lead_m ! ocean surface (ssh_m) if ice is not embedded 66 77 ! ! ocean surface in leads if ice is embedded 67 78 #endif 79 INTEGER :: jk ! vertical loop index 80 INTEGER :: Kmm ! ocean time levelindex 81 ! 68 82 ! copy nemo forcing arrays into iceberg versions with extra halo 69 83 ! only necessary for variables not on T points 70 84 ! and ssh which is used to calculate gradients 71 72 uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 73 vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 85 ! 86 ! surface forcing 87 ! 88 ssu_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 89 ssv_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 90 sst_e(1:jpi,1:jpj) = sst_m(:,:) 91 sss_e(1:jpi,1:jpj) = sss_m(:,:) 92 fr_e (1:jpi,1:jpj) = fr_i (:,:) 93 ua_e (1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk 94 va_e (1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk 74 95 ff_e(1:jpi,1:jpj) = ff_f (:,:) 75 tt_e(1:jpi,1:jpj) = sst_m(:,:) 76 ss_e(1:jpi,1:jpj) = sss_m(:,:) 77 fr_e(1:jpi,1:jpj) = fr_i (:,:) 78 ua_e(1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk 79 va_e(1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk 80 ! 81 CALL lbc_lnk_icb( 'icbutl', uo_e, 'U', -1._wp, 1, 1 ) 82 CALL lbc_lnk_icb( 'icbutl', vo_e, 'V', -1._wp, 1, 1 ) 83 CALL lbc_lnk_icb( 'icbutl', ff_e, 'F', +1._wp, 1, 1 ) 84 CALL lbc_lnk_icb( 'icbutl', ua_e, 'U', -1._wp, 1, 1 ) 85 CALL lbc_lnk_icb( 'icbutl', va_e, 'V', -1._wp, 1, 1 ) 86 CALL lbc_lnk_icb( 'icbutl', fr_e, 'T', +1._wp, 1, 1 ) 87 CALL lbc_lnk_icb( 'icbutl', tt_e, 'T', +1._wp, 1, 1 ) 88 CALL lbc_lnk_icb( 'icbutl', ss_e, 'T', +1._wp, 1, 1 ) 96 ! 97 CALL lbc_lnk_icb( 'icbutl', ssu_e, 'U', -1._wp, 1, 1 ) 98 CALL lbc_lnk_icb( 'icbutl', ssv_e, 'V', -1._wp, 1, 1 ) 99 CALL lbc_lnk_icb( 'icbutl', ua_e , 'U', -1._wp, 1, 1 ) 100 CALL lbc_lnk_icb( 'icbutl', va_e , 'V', -1._wp, 1, 1 ) 89 101 #if defined key_si3 90 102 hi_e(1:jpi, 1:jpj) = hm_i (:,:) … … 96 108 ssh_e(1:jpi, 1:jpj) = zssh_lead_m(:,:) * tmask(:,:,1) 97 109 ! 98 CALL lbc_lnk_icb( 'icbutl', hi_e , 'T', +1._wp, 1, 1 )99 110 CALL lbc_lnk_icb( 'icbutl', ui_e , 'U', -1._wp, 1, 1 ) 100 111 CALL lbc_lnk_icb( 'icbutl', vi_e , 'V', -1._wp, 1, 1 ) 101 112 #else 102 ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1) 113 ssh_e(1:jpi, 1:jpj) = ssh_m(:,:) * tmask(:,:,1) 103 114 #endif 104 CALL lbc_lnk_icb( 'icbutl', ssh_e, 'T', +1._wp, 1, 1 ) 115 ! 116 ! (PM) could be improve with a 3d lbclnk gathering both variables 117 ! should be done once extra haloe generalised 118 IF ( ln_M2016 ) THEN 119 DO jk = 1,jpk 120 ! uoce 121 ztmp(1:jpi,1:jpj) = uu(:,:,jk,Kmm) 122 CALL lbc_lnk_icb( 'icbutl', ztmp, 'U', -1._wp, 1, 1 ) 123 uoce_e(:,:,jk) = ztmp(:,:) 124 ! 125 ! voce 126 ztmp(1:jpi,1:jpj) = vv(:,:,jk,Kmm) 127 CALL lbc_lnk_icb( 'icbutl', ztmp, 'V', -1._wp, 1, 1 ) 128 voce_e(:,:,jk) = ztmp(:,:) 129 END DO 130 toce_e(1:jpi,1:jpj,1:jpk) = ts(:,:,:,1,Kmm) 131 e3t_e (1:jpi,1:jpj,1:jpk) = e3t(:,:,:,Kmm) 132 END IF 105 133 ! 106 134 END SUBROUTINE icb_utl_copy 107 135 108 136 109 SUBROUTINE icb_utl_interp( pi, pe1, puo, pui, pua, pssh_i, & 110 & pj, pe2, pvo, pvi, pva, pssh_j, & 111 & psst, pcn, phi, pff, psss ) 137 SUBROUTINE icb_utl_interp( pi, pj, pe1 , pssu, pui, pua, pssh_i, & 138 & pe2 , pssv, pvi, pva, pssh_j, & 139 & psst, psss, pcn, phi, pff , & 140 & plon, plat, ptoce, puoce, pvoce, pe3t ) 112 141 !!---------------------------------------------------------------------- 113 142 !! *** ROUTINE icb_utl_interp *** … … 127 156 !!---------------------------------------------------------------------- 128 157 REAL(wp), INTENT(in ) :: pi , pj ! position in (i,j) referential 129 REAL(wp), INTENT( out) :: pe1, pe2 ! i- and j scale factors 130 REAL(wp), INTENT( out) :: puo, pvo, pui, pvi, pua, pva ! ocean, ice and wind speeds 131 REAL(wp), INTENT( out) :: pssh_i, pssh_j ! ssh i- & j-gradients 132 REAL(wp), INTENT( out) :: psst, pcn, phi, pff, psss ! SST, ice concentration, ice thickness, Coriolis, SSS 133 ! 158 REAL(wp), INTENT( out), OPTIONAL :: pe1, pe2 ! i- and j scale factors 159 REAL(wp), INTENT( out), OPTIONAL :: pssu, pssv, pui, pvi, pua, pva ! ocean, ice and wind speeds 160 REAL(wp), INTENT( out), OPTIONAL :: pssh_i, pssh_j ! ssh i- & j-gradients 161 REAL(wp), INTENT( out), OPTIONAL :: psst, psss, pcn, phi, pff ! SST, SSS, ice concentration, ice thickness, Coriolis 162 REAL(wp), INTENT( out), OPTIONAL :: plat, plon ! position 163 REAL(wp), DIMENSION(jpk), INTENT( out), OPTIONAL :: ptoce, puoce, pvoce, pe3t ! 3D variables 164 ! 165 REAL(wp), DIMENSION(4) :: zwT , zwU , zwV , zwF ! interpolation weight 166 REAL(wp), DIMENSION(4) :: zmskF, zmskU, zmskV, zmskT ! mask 167 REAL(wp), DIMENSION(4) :: zwTp, zmskTp, zwTm, zmskTm 168 REAL(wp), DIMENSION(4,jpk) :: zw1d 169 INTEGER :: iiT, iiU, iiV, iiF, ijT, ijU, ijV, ijF ! bottom left corner 170 INTEGER :: iiTp, iiTm, ijTp, ijTm 134 171 REAL(wp) :: zcd, zmod ! local scalars 135 172 !!---------------------------------------------------------------------- 136 137 pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj ) ! scale factors 138 pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 139 ! 140 puo = icb_utl_bilin_h( uo_e, pi, pj, 'U', .false. ) ! ocean velocities 141 pvo = icb_utl_bilin_h( vo_e, pi, pj, 'V', .false. ) 142 psst = icb_utl_bilin_h( tt_e, pi, pj, 'T', .true. ) ! SST 143 psss = icb_utl_bilin_h( ss_e, pi, pj, 'T', .true. ) ! SSS 144 pcn = icb_utl_bilin_h( fr_e, pi, pj, 'T', .true. ) ! ice concentration 145 pff = icb_utl_bilin_h( ff_e, pi, pj, 'F', .false. ) ! Coriolis parameter 146 ! 147 pua = icb_utl_bilin_h( ua_e, pi, pj, 'U', .true. ) ! 10m wind 148 pva = icb_utl_bilin_h( va_e, pi, pj, 'V', .true. ) ! here (ua,va) are stress => rough conversion from stress to speed 149 zcd = 1.22_wp * 1.5e-3_wp ! air density * drag coefficient 150 zmod = 1._wp / MAX( 1.e-20, SQRT( zcd * SQRT( pua*pua + pva*pva) ) ) 151 pua = pua * zmod ! note: stress module=0 necessarly implies ua=va=0 152 pva = pva * zmod 153 173 ! 174 ! get position, weight and mask 175 CALL icb_utl_pos( pi, pj, 'T', iiT, ijT, zwT, zmskT ) 176 CALL icb_utl_pos( pi, pj, 'U', iiU, ijU, zwU, zmskU ) 177 CALL icb_utl_pos( pi, pj, 'V', iiV, ijV, zwV, zmskV ) 178 CALL icb_utl_pos( pi, pj, 'F', iiF, ijF, zwF, zmskF ) 179 ! 180 ! metrics and coordinates 181 IF ( PRESENT(pe1 ) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj ) ! scale factors 182 IF ( PRESENT(pe2 ) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 183 IF ( PRESENT(plon) ) plon= icb_utl_bilin_h( rlon_e, iiT, ijT, zwT, .true. ) 184 IF ( PRESENT(plat) ) plat= icb_utl_bilin_h( rlat_e, iiT, ijT, zwT, .false. ) 185 ! 186 IF ( PRESENT(pssu) ) pssu = icb_utl_bilin_h( ssu_e, iiU, ijU, zwU , .false. ) ! ocean velocities 187 IF ( PRESENT(pssv) ) pssv = icb_utl_bilin_h( ssv_e, iiV, ijV, zwV , .false. ) ! 188 IF ( PRESENT(psst) ) psst = icb_utl_bilin_h( sst_e, iiT, ijT, zwT * zmskT, .false. ) ! sst 189 IF ( PRESENT(psss) ) psss = icb_utl_bilin_h( sss_e, iiT, ijT, zwT * zmskT, .false. ) ! sss 190 IF ( PRESENT(pcn ) ) pcn = icb_utl_bilin_h( fr_e , iiT, ijT, zwT * zmskT, .false. ) ! ice concentration 191 IF ( PRESENT(pff ) ) pff = icb_utl_bilin_h( ff_e , iiF, ijF, zwF , .false. ) ! Coriolis parameter 192 ! 193 IF ( PRESENT(pua) .AND. PRESENT(pva) ) THEN 194 pua = icb_utl_bilin_h( ua_e, iiU, ijU, zwU * zmskU, .false. ) ! 10m wind 195 pva = icb_utl_bilin_h( va_e, iiV, ijV, zwV * zmskV, .false. ) ! here (ua,va) are stress => rough conversion from stress to speed 196 zcd = 1.22_wp * 1.5e-3_wp ! air density * drag coefficient 197 zmod = 1._wp / MAX( 1.e-20, SQRT( zcd * SQRT( pua*pua + pva*pva) ) ) 198 pua = pua * zmod ! note: stress module=0 necessarly implies ua=va=0 199 pva = pva * zmod 200 END IF 201 ! 154 202 #if defined key_si3 155 pui = icb_utl_bilin_h( ui_e , pi, pj, 'U', .false. )! sea-ice velocities156 pvi = icb_utl_bilin_h( vi_e , pi, pj, 'V', .false. )157 phi = icb_utl_bilin_h( hi_e , pi, pj, 'T', .true. )! ice thickness203 IF ( PRESENT(pui) ) pui = icb_utl_bilin_h( ui_e , iiU, ijU, zwU , .false. ) ! sea-ice velocities 204 IF ( PRESENT(pvi) ) pvi = icb_utl_bilin_h( vi_e , iiV, ijV, zwV , .false. ) 205 IF ( PRESENT(phi) ) phi = icb_utl_bilin_h( hi_e , iiT, ijT, zwT * zmskT, .false. ) ! ice thickness 158 206 #else 159 pui = 0._wp160 pvi = 0._wp161 phi = 0._wp207 IF ( PRESENT(pui) ) pui = 0._wp 208 IF ( PRESENT(pvi) ) pvi = 0._wp 209 IF ( PRESENT(phi) ) phi = 0._wp 162 210 #endif 163 211 ! 164 212 ! Estimate SSH gradient in i- and j-direction (centred evaluation) 165 pssh_i = ( icb_utl_bilin_h( ssh_e, pi+0.1_wp, pj, 'T', .true. ) - & 166 & icb_utl_bilin_h( ssh_e, pi-0.1_wp, pj, 'T', .true. ) ) / ( 0.2_wp * pe1 ) 167 pssh_j = ( icb_utl_bilin_h( ssh_e, pi, pj+0.1_wp, 'T', .true. ) - & 168 & icb_utl_bilin_h( ssh_e, pi, pj-0.1_wp, 'T', .true. ) ) / ( 0.2_wp * pe2 ) 213 IF ( PRESENT(pssh_i) .AND. PRESENT(pssh_j) ) THEN 214 CALL icb_utl_pos( pi+0.1, pj , 'T', iiTp, ijTp, zwTp, zmskTp ) 215 CALL icb_utl_pos( pi-0.1, pj , 'T', iiTm, ijTm, zwTm, zmskTm ) 216 ! 217 IF ( .NOT. PRESENT(pe1) ) pe1 = icb_utl_bilin_e( e1t, e1u, e1v, e1f, pi, pj ) 218 pssh_i = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) - & 219 & icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. ) ) / ( 0.2_wp * pe1 ) 220 ! 221 CALL icb_utl_pos( pi , pj+0.1, 'T', iiTp, ijTp, zwTp, zmskTp ) 222 CALL icb_utl_pos( pi , pj-0.1, 'T', iiTm, ijTm, zwTm, zmskTm ) 223 ! 224 IF ( .NOT. PRESENT(pe2) ) pe2 = icb_utl_bilin_e( e2t, e2u, e2v, e2f, pi, pj ) 225 pssh_j = ( icb_utl_bilin_h( ssh_e, iiTp, ijTp, zwTp*zmskTp, .false. ) - & 226 & icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. ) ) / ( 0.2_wp * pe2 ) 227 END IF 228 ! 229 ! 3d interpolation 230 IF ( PRESENT(puoce) .AND. PRESENT(pvoce) ) THEN 231 ! no need to mask as 0 is a valid data for land 232 zw1d(1,:) = zwU(1) ; zw1d(2,:) = zwU(2) ; zw1d(3,:) = zwU(3) ; zw1d(4,:) = zwU(4) ; 233 puoce(:) = icb_utl_bilin_h( uoce_e , iiU, ijU, zw1d ) 234 235 zw1d(1,:) = zwV(1) ; zw1d(2,:) = zwV(2) ; zw1d(3,:) = zwV(3) ; zw1d(4,:) = zwV(4) ; 236 pvoce(:) = icb_utl_bilin_h( voce_e , iiV, ijV, zw1d ) 237 END IF 238 239 IF ( PRESENT(ptoce) ) THEN 240 ! for temperature we need to mask the weight properly 241 ! no need of extra halo as it is a T point variable 242 zw1d(1,:) = tmask(iiT ,ijT ,:) * zwT(1) * zmskT(1) 243 zw1d(2,:) = tmask(iiT+1,ijT ,:) * zwT(2) * zmskT(2) 244 zw1d(3,:) = tmask(iiT ,ijT+1,:) * zwT(3) * zmskT(3) 245 zw1d(4,:) = tmask(iiT+1,ijT+1,:) * zwT(4) * zmskT(4) 246 ptoce(:) = icb_utl_bilin_h( toce_e , iiT, ijT, zw1d ) 247 END IF 248 ! 249 IF ( PRESENT(pe3t) ) pe3t(:) = e3t_e(iiT,ijT,:) ! as in Nacho tarball need to be fix once we are able to reproduce Nacho results 169 250 ! 170 251 END SUBROUTINE icb_utl_interp 171 252 172 173 REAL(wp) FUNCTION icb_utl_bilin_h( pfld, pi, pj, cd_type, plmask ) 253 SUBROUTINE icb_utl_pos( pi, pj, cd_type, kii, kij, pw, pmsk ) 174 254 !!---------------------------------------------------------------------- 175 255 !! *** FUNCTION icb_utl_bilin *** … … 182 262 !! 183 263 !!---------------------------------------------------------------------- 184 REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) :: pfld ! field to be interpolated 185 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 186 CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points 187 LOGICAL , INTENT(in) :: plmask ! special treatment of mask point 188 ! 189 INTEGER :: ii, ij ! local integer 190 REAL(wp) :: zi, zj ! local real 191 REAL(wp) :: zw1, zw2, zw3, zw4 192 REAL(wp), DIMENSION(4) :: zmask 264 REAL(wp) , INTENT(IN) :: pi, pj ! targeted coordinates in (i,j) referential 265 CHARACTER(len=1) , INTENT(IN) :: cd_type ! point type 266 REAL(wp), DIMENSION(4), INTENT(OUT) :: pw, pmsk ! weight and mask 267 INTEGER , INTENT(OUT) :: kii, kij ! bottom left corner position in local domain 268 ! 269 REAL(wp) :: zwi, zwj ! distance to bottom left corner 270 INTEGER :: ierr 271 ! 193 272 !!---------------------------------------------------------------------- 194 273 ! … … 198 277 ! since we're looking for four T points containing quadrant we're in of 199 278 ! current T cell 200 ii = MAX(0, INT( pi))201 ij = MAX(0, INT( pj)) ! T-point202 z i = pi - REAL(ii,wp)203 z j = pj - REAL(ij,wp)279 kii = MAX(0, INT( pi )) 280 kij = MAX(0, INT( pj )) ! T-point 281 zwi = pi - REAL(kii,wp) 282 zwj = pj - REAL(kij,wp) 204 283 CASE ( 'U' ) 205 ii = MAX(0, INT( pi-0.5_wp ))206 ij = MAX(0, INT( pj)) ! U-point207 z i = pi - 0.5_wp - REAL(ii,wp)208 z j = pj - REAL(ij,wp)284 kii = MAX(0, INT( pi-0.5_wp )) 285 kij = MAX(0, INT( pj )) ! U-point 286 zwi = pi - 0.5_wp - REAL(kii,wp) 287 zwj = pj - REAL(kij,wp) 209 288 CASE ( 'V' ) 210 ii = MAX(0, INT( pi))211 ij = MAX(0, INT( pj-0.5_wp )) ! V-point212 z i = pi - REAL(ii,wp)213 z j = pj - 0.5_wp - REAL(ij,wp)289 kii = MAX(0, INT( pi )) 290 kij = MAX(0, INT( pj-0.5_wp )) ! V-point 291 zwi = pi - REAL(kii,wp) 292 zwj = pj - 0.5_wp - REAL(kij,wp) 214 293 CASE ( 'F' ) 215 ii = MAX(0, INT( pi-0.5_wp ))216 ij = MAX(0, INT( pj-0.5_wp )) ! F-point217 z i = pi - 0.5_wp - REAL(ii,wp)218 z j = pj - 0.5_wp - REAL(ij,wp)294 kii = MAX(0, INT( pi-0.5_wp )) 295 kij = MAX(0, INT( pj-0.5_wp )) ! F-point 296 zwi = pi - 0.5_wp - REAL(kii,wp) 297 zwj = pj - 0.5_wp - REAL(kij,wp) 219 298 END SELECT 299 ! 300 ! compute weight 301 pw(1) = (1._wp-zwi) * (1._wp-zwj) 302 pw(2) = zwi * (1._wp-zwj) 303 pw(3) = (1._wp-zwi) * zwj 304 pw(4) = zwi * zwj 305 ! 306 ! find position in this processor. Prevent near edge problems (see #1389) 307 ! 308 IF (TRIM(cd_type) == 'T' ) THEN 309 ierr = 0 310 IF ( kii < mig( 1 ) ) THEN ; ierr = ierr + 1 311 ELSEIF( kii >= mig(jpi) ) THEN ; ierr = ierr + 1 312 ENDIF 313 ! 314 IF ( kij < mjg( 1 ) ) THEN ; ierr = ierr + 1 315 ELSEIF( kij >= mjg(jpj) ) THEN ; ierr = ierr + 1 316 ENDIF 317 ! 318 IF ( ierr > 0 ) THEN 319 WRITE(numout,*) 'bottom left corner T point out of bound' 320 WRITE(numout,*) pi, kii, mig( 1 ), mig(jpi) 321 WRITE(numout,*) pj, kij, mjg( 1 ), mjg(jpj) 322 WRITE(numout,*) pmsk 323 CALL ctl_stop('STOP','icb_utl_bilin_h: an icebergs coordinates is out of valid range (out of bound error)') 324 END IF 325 END IF 220 326 ! 221 327 ! find position in this processor. Prevent near edge problems (see #1389) 222 328 ! (PM) will be useless if extra halo is used in NEMO 223 329 ! 224 IF ( ii <= mig(1)-1 ) THEN ;ii = 0225 ELSEIF( ii > mig(jpi) ) THEN ;ii = jpi226 ELSE ; ii = mi1(ii)330 IF ( kii <= mig(1)-1 ) THEN ; kii = 0 331 ELSEIF( kii > mig(jpi) ) THEN ; kii = jpi 332 ELSE ; kii = mi1(kii) 227 333 ENDIF 228 IF ( ij <= mjg(1)-1 ) THEN ;ij = 0229 ELSEIF( ij > mjg(jpj) ) THEN ;ij = jpj230 ELSE ; ij = mj1(ij)334 IF ( kij <= mjg(1)-1 ) THEN ; kij = 0 335 ELSEIF( kij > mjg(jpj) ) THEN ; kij = jpj 336 ELSE ; kij = mj1(kij) 231 337 ENDIF 232 338 ! 233 339 ! define mask array 234 IF (plmask) THEN 235 ! land value is not used in the interpolation 236 SELECT CASE ( cd_type ) 237 CASE ( 'T' ) 238 zmask = (/tmask_e(ii,ij), tmask_e(ii+1,ij), tmask_e(ii,ij+1), tmask_e(ii+1,ij+1)/) 239 CASE ( 'U' ) 240 zmask = (/umask_e(ii,ij), umask_e(ii+1,ij), umask_e(ii,ij+1), umask_e(ii+1,ij+1)/) 241 CASE ( 'V' ) 242 zmask = (/vmask_e(ii,ij), vmask_e(ii+1,ij), vmask_e(ii,ij+1), vmask_e(ii+1,ij+1)/) 243 CASE ( 'F' ) 244 ! F case only used for coriolis, ff_f is not mask so zmask = 1 245 zmask = 1. 246 END SELECT 247 ELSE 248 ! land value is used during interpolation 249 zmask = 1. 250 END iF 251 ! 252 ! compute weight 253 zw1 = zmask(1) * (1._wp-zi) * (1._wp-zj) 254 zw2 = zmask(2) * zi * (1._wp-zj) 255 zw3 = zmask(3) * (1._wp-zi) * zj 256 zw4 = zmask(4) * zi * zj 257 ! 258 ! compute interpolated value 259 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) 260 ! 261 END FUNCTION icb_utl_bilin_h 262 263 264 REAL(wp) FUNCTION icb_utl_bilin( pfld, pi, pj, cd_type ) 340 ! land value is not used in the interpolation 341 SELECT CASE ( cd_type ) 342 CASE ( 'T' ) 343 pmsk = (/tmask_e(kii,kij), tmask_e(kii+1,kij), tmask_e(kii,kij+1), tmask_e(kii+1,kij+1)/) 344 CASE ( 'U' ) 345 pmsk = (/umask_e(kii,kij), umask_e(kii+1,kij), umask_e(kii,kij+1), umask_e(kii+1,kij+1)/) 346 CASE ( 'V' ) 347 pmsk = (/vmask_e(kii,kij), vmask_e(kii+1,kij), vmask_e(kii,kij+1), vmask_e(kii+1,kij+1)/) 348 CASE ( 'F' ) 349 ! F case only used for coriolis, ff_f is not mask so zmask = 1 350 pmsk = 1. 351 END SELECT 352 END SUBROUTINE icb_utl_pos 353 354 REAL(wp) FUNCTION icb_utl_bilin_2d_h( pfld, pii, pij, pw, pllon ) 265 355 !!---------------------------------------------------------------------- 266 356 !! *** FUNCTION icb_utl_bilin *** 267 357 !! 268 358 !! ** Purpose : bilinear interpolation at berg location depending on the grid-point type 359 !! this version deals with extra halo points 269 360 !! 270 361 !! !!gm CAUTION an optional argument should be added to handle … … 272 363 !! 273 364 !!---------------------------------------------------------------------- 274 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pfld ! field to be interpolated 275 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 276 CHARACTER(len=1) , INTENT(in) :: cd_type ! type of pfld array grid-points: = T , U , V or F points 277 ! 278 INTEGER :: ii, ij ! local integer 279 REAL(wp) :: zi, zj ! local real 280 !!---------------------------------------------------------------------- 281 ! 282 SELECT CASE ( cd_type ) 283 CASE ( 'T' ) 284 ! note that here there is no +0.5 added 285 ! since we're looking for four T points containing quadrant we're in of 286 ! current T cell 287 ii = MAX(1, INT( pi )) 288 ij = MAX(1, INT( pj )) ! T-point 289 zi = pi - REAL(ii,wp) 290 zj = pj - REAL(ij,wp) 291 CASE ( 'U' ) 292 ii = MAX(1, INT( pi-0.5 )) 293 ij = MAX(1, INT( pj )) ! U-point 294 zi = pi - 0.5 - REAL(ii,wp) 295 zj = pj - REAL(ij,wp) 296 CASE ( 'V' ) 297 ii = MAX(1, INT( pi )) 298 ij = MAX(1, INT( pj-0.5 )) ! V-point 299 zi = pi - REAL(ii,wp) 300 zj = pj - 0.5 - REAL(ij,wp) 301 CASE ( 'F' ) 302 ii = MAX(1, INT( pi-0.5 )) 303 ij = MAX(1, INT( pj-0.5 )) ! F-point 304 zi = pi - 0.5 - REAL(ii,wp) 305 zj = pj - 0.5 - REAL(ij,wp) 306 END SELECT 307 ! 308 ! find position in this processor. Prevent near edge problems (see #1389) 309 IF ( ii < mig( 1 ) ) THEN ; ii = 1 310 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 311 ELSE ; ii = mi1(ii) 365 REAL(wp), DIMENSION(0:jpi+1,0:jpj+1), INTENT(in) :: pfld ! field to be interpolated 366 REAL(wp), DIMENSION(4) , INTENT(in) :: pw ! weight 367 LOGICAL , INTENT(in) :: pllon ! input data is a longitude 368 INTEGER , INTENT(in) :: pii, pij ! bottom left corner 369 ! 370 REAL(wp), DIMENSION(4) :: zdat ! input data 371 !!---------------------------------------------------------------------- 372 ! 373 ! data 374 zdat(1) = pfld(pii ,pij ) 375 zdat(2) = pfld(pii+1,pij ) 376 zdat(3) = pfld(pii ,pij+1) 377 zdat(4) = pfld(pii+1,pij+1) 378 ! 379 IF( pllon .AND. MAXVAL(zdat) - MINVAL(zdat) > 90._wp ) THEN 380 WHERE( zdat < 0._wp ) zdat = zdat + 360._wp 312 381 ENDIF 313 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 314 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 315 ELSE ; ij = mj1(ij) 316 ENDIF 317 ! 318 IF( ii == jpi ) ii = ii-1 319 IF( ij == jpj ) ij = ij-1 320 ! 321 icb_utl_bilin = ( pfld(ii,ij ) * (1.-zi) + pfld(ii+1,ij ) * zi ) * (1.-zj) & 322 & + ( pfld(ii,ij+1) * (1.-zi) + pfld(ii+1,ij+1) * zi ) * zj 323 ! 324 END FUNCTION icb_utl_bilin 325 326 327 REAL(wp) FUNCTION icb_utl_bilin_x( pfld, pi, pj ) 328 !!---------------------------------------------------------------------- 329 !! *** FUNCTION icb_utl_bilin_x *** 382 ! 383 ! compute interpolated value 384 icb_utl_bilin_2d_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)) 385 ! 386 IF( pllon .AND. icb_utl_bilin_2d_h > 180._wp ) icb_utl_bilin_2d_h = icb_utl_bilin_2d_h - 360._wp 387 ! 388 END FUNCTION icb_utl_bilin_2d_h 389 390 FUNCTION icb_utl_bilin_3d_h( pfld, pii, pij, pw ) 391 !!---------------------------------------------------------------------- 392 !! *** FUNCTION icb_utl_bilin *** 330 393 !! 331 394 !! ** Purpose : bilinear interpolation at berg location depending on the grid-point type 332 !! Special case for interpolating longitude395 !! this version deals with extra halo points 333 396 !! 334 397 !! !!gm CAUTION an optional argument should be added to handle … … 336 399 !! 337 400 !!---------------------------------------------------------------------- 338 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pfld ! field to be interpolated 339 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 340 ! 341 INTEGER :: ii, ij ! local integer 342 REAL(wp) :: zi, zj ! local real 343 REAL(wp) :: zret ! local real 344 REAL(wp), DIMENSION(4) :: z4 345 !!---------------------------------------------------------------------- 346 ! 347 ! note that here there is no +0.5 added 348 ! since we're looking for four T points containing quadrant we're in of 349 ! current T cell 350 ii = MAX(1, INT( pi )) 351 ij = MAX(1, INT( pj )) ! T-point 352 zi = pi - REAL(ii,wp) 353 zj = pj - REAL(ij,wp) 354 ! 355 ! find position in this processor. Prevent near edge problems (see #1389) 356 IF ( ii < mig( 1 ) ) THEN ; ii = 1 357 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi 358 ELSE ; ii = mi1(ii) 359 ENDIF 360 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 361 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj 362 ELSE ; ij = mj1(ij) 363 ENDIF 364 ! 365 IF( ii == jpi ) ii = ii-1 366 IF( ij == jpj ) ij = ij-1 367 ! 368 z4(1) = pfld(ii ,ij ) 369 z4(2) = pfld(ii+1,ij ) 370 z4(3) = pfld(ii ,ij+1) 371 z4(4) = pfld(ii+1,ij+1) 372 IF( MAXVAL(z4) - MINVAL(z4) > 90._wp ) THEN 373 WHERE( z4 < 0._wp ) z4 = z4 + 360._wp 374 ENDIF 375 ! 376 zret = (z4(1) * (1.-zi) + z4(2) * zi) * (1.-zj) + (z4(3) * (1.-zi) + z4(4) * zi) * zj 377 IF( zret > 180._wp ) zret = zret - 360._wp 378 icb_utl_bilin_x = zret 379 ! 380 END FUNCTION icb_utl_bilin_x 381 401 REAL(wp), DIMENSION(0:jpi+1,0:jpj+1, jpk), INTENT(in) :: pfld ! field to be interpolated 402 REAL(wp), DIMENSION(4,jpk) , INTENT(in) :: pw ! weight 403 INTEGER , INTENT(in) :: pii, pij ! bottom left corner 404 REAL(wp), DIMENSION(jpk) :: icb_utl_bilin_3d_h 405 ! 406 REAL(wp), DIMENSION(4,jpk) :: zdat ! input data 407 INTEGER :: jk 408 !!---------------------------------------------------------------------- 409 ! 410 ! data 411 zdat(1,:) = pfld(pii ,pij ,:) 412 zdat(2,:) = pfld(pii+1,pij ,:) 413 zdat(3,:) = pfld(pii ,pij+1,:) 414 zdat(4,:) = pfld(pii+1,pij+1,:) 415 ! 416 ! compute interpolated value 417 DO jk=1,jpk 418 icb_utl_bilin_3d_h(jk) = ( zdat(1,jk)*pw(1,jk) + zdat(2,jk)*pw(2,jk) + zdat(3,jk)*pw(3,jk) + zdat(4,jk)*pw(4,jk) ) & 419 & / MAX(1.e-20, pw(1,jk)+pw(2,jk)+pw(3,jk)+pw(4,jk)) 420 END DO 421 ! 422 END FUNCTION icb_utl_bilin_3d_h 382 423 383 424 REAL(wp) FUNCTION icb_utl_bilin_e( pet, peu, pev, pef, pi, pj ) … … 390 431 !!---------------------------------------------------------------------- 391 432 REAL(wp), DIMENSION(:,:), INTENT(in) :: pet, peu, pev, pef ! horizontal scale factor to be interpolated at t-,u-,v- & f-pts 392 REAL(wp) , INTENT(in) :: pi, pj ! targeted coordinates in (i,j) referential 393 ! 394 INTEGER :: ii, ij, icase, ierr ! local integer 433 REAL(wp) , INTENT(IN) :: pi , pj ! iceberg position 395 434 ! 396 435 ! weights corresponding to corner points of a T cell quadrant 397 436 REAL(wp) :: zi, zj ! local real 437 INTEGER :: ii, ij ! bottom left corner coordinate in local domain 398 438 ! 399 439 ! values at corner points of a T cell quadrant … … 402 442 !!---------------------------------------------------------------------- 403 443 ! 444 ! cannot used iiT because need ii/ij reltaive to global indices not local one 404 445 ii = MAX(1, INT( pi )) ; ij = MAX(1, INT( pj )) ! left bottom T-point (i,j) indices 405 446 ! 406 447 ! fractional box spacing 407 448 ! 0 <= zi < 0.5, 0 <= zj < 0.5 --> NW quadrant of current T cell … … 413 454 zj = pj - REAL(ij,wp) 414 455 415 ! find position in this processor. Prevent near edge problems (see #1389) 416 ! 417 ierr = 0 418 IF ( ii < mig( 1 ) ) THEN ; ii = 1 ; ierr = ierr + 1 419 ELSEIF( ii > mig(jpi) ) THEN ; ii = jpi ; ierr = ierr + 1 420 ELSE ; ii = mi1(ii) 421 ENDIF 422 IF ( ij < mjg( 1 ) ) THEN ; ij = 1 ; ierr = ierr + 1 423 ELSEIF( ij > mjg(jpj) ) THEN ; ij = jpj ; ierr = ierr + 1 424 ELSE ; ij = mj1(ij) 425 ENDIF 426 ! 427 IF( ii == jpi ) THEN ; ii = ii-1 ; ierr = ierr + 1 ; END IF 428 IF( ij == jpj ) THEN ; ij = ij-1 ; ierr = ierr + 1 ; END IF 429 ! 430 IF ( ierr > 0 ) CALL ctl_stop('STOP','icb_utl_bilin_e: an icebergs coordinates is out of valid range (out of bound error)') 456 ! conversion to local domain (no need to do a sanity check already done in icbpos) 457 ii = mi1(ii) 458 ij = mj1(ij) 431 459 ! 432 460 IF( 0.0_wp <= zi .AND. zi < 0.5_wp ) THEN … … 465 493 END FUNCTION icb_utl_bilin_e 466 494 495 SUBROUTINE icb_utl_getkb( kb, pe3, pD ) 496 !!---------------------------------------------------------------------- 497 !! *** ROUTINE icb_utl_getkb *** 498 !! 499 !! ** Purpose : compute the latest level affected by icb 500 !! 501 !!---------------------------------------------------------------------- 502 INTEGER, INTENT(out):: kb 503 REAL(wp), DIMENSION(:), INTENT(in) :: pe3 504 REAL(wp), INTENT(in) :: pD 505 !! 506 INTEGER :: jk 507 REAL(wp) :: zdepw 508 !!---------------------------------------------------------------------- 509 !! 510 zdepw = pe3(1) ; kb = 2 511 DO WHILE ( zdepw < pD) 512 zdepw = zdepw + pe3(kb) 513 kb = kb + 1 514 END DO 515 kb = MIN(kb - 1,jpk) 516 END SUBROUTINE 517 518 SUBROUTINE icb_utl_zavg(pzavg, pdat, pe3, pD, kb ) 519 !!---------------------------------------------------------------------- 520 !! *** ROUTINE icb_utl_getkb *** 521 !! 522 !! ** Purpose : compute the vertical average of ocean properties affected by icb 523 !! 524 !!---------------------------------------------------------------------- 525 INTEGER, INTENT(in ) :: kb ! deepest level affected by icb 526 REAL(wp), DIMENSION(:), INTENT(in ) :: pe3, pdat ! vertical profile 527 REAL(wp), INTENT(in ) :: pD ! draft 528 REAL(wp), INTENT(out) :: pzavg ! z average 529 !!---------------------------------------------------------------------- 530 INTEGER :: jk 531 REAL(wp) :: zdep 532 !!---------------------------------------------------------------------- 533 pzavg = 0.0 ; zdep = 0.0 534 DO jk = 1,kb-1 535 pzavg = pzavg + pe3(jk)*pdat(jk) 536 zdep = zdep + pe3(jk) 537 END DO 538 ! if kb is limited by mbkt => bottom value is used between bathy and icb tail 539 ! if kb not limited by mbkt => ocean value over mask is used (ie 0.0 for u, v) 540 pzavg = ( pzavg + (pD - zdep)*pdat(kb)) / pD 541 END SUBROUTINE 467 542 468 543 SUBROUTINE icb_utl_add( bergvals, ptvals ) … … 653 728 WRITE(numicb, 9200) kt, berg%number(1), & 654 729 pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel, & 655 pt% uo, pt%vo, pt%ua, pt%va, pt%ui, pt%vi730 pt%ssu, pt%ssv, pt%ua, pt%va, pt%ui, pt%vi 656 731 CALL flush( numicb ) 657 732 9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4)) … … 679 754 WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea 680 755 WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' ) & 681 & 'timestep', 'number', 'xi,yj','lon,lat','u,v',' uo,vo','ua,va','ui,vi'756 & 'timestep', 'number', 'xi,yj','lon,lat','u,v','ssu,ssv','ua,va','ui,vi' 682 757 ENDIF 683 758 DO WHILE( ASSOCIATED(this) ) … … 823 898 END FUNCTION icb_utl_heat 824 899 900 SUBROUTINE test_icb_utl_getkb 901 !!---------------------------------------------------------------------- 902 !! *** FUNCTION test_icb_utl_getkb *** 903 !! 904 !! ** Purpose : Test routine icb_utl_getkb, icb_utl_zavg 905 !! ** Methode : Call each subroutine with specific input data 906 !! What should be the output is easy to determined and check 907 !! if NEMO return the correct answer. 908 !! ** Comments : not called, if needed a CALL test_icb_utl_getkb need to be added in icb_step 909 !!---------------------------------------------------------------------- 910 INTEGER :: ikb 911 REAL(wp) :: zD, zout 912 REAL(wp), DIMENSION(jpk) :: ze3, zin 913 WRITE(numout,*) 'Test icb_utl_getkb : ' 914 zD = 0.0 ; ze3= 20.0 915 WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1) 916 CALL icb_utl_getkb(ikb, ze3, zD) 917 WRITE(numout,*) 'OUTPUT : kb = ',ikb 918 919 zD = 8000000.0 ; ze3= 20.0 920 WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1) 921 CALL icb_utl_getkb(ikb, ze3, zD) 922 WRITE(numout,*) 'OUTPUT : kb = ',ikb 923 924 zD = 80.0 ; ze3= 20.0 925 WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1) 926 CALL icb_utl_getkb(ikb, ze3, zD) 927 WRITE(numout,*) 'OUTPUT : kb = ',ikb 928 929 zD = 85.0 ; ze3= 20.0 930 WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1) 931 CALL icb_utl_getkb(ikb, ze3, zD) 932 WRITE(numout,*) 'OUTPUT : kb = ',ikb 933 934 zD = 75.0 ; ze3= 20.0 935 WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1) 936 CALL icb_utl_getkb(ikb, ze3, zD) 937 WRITE(numout,*) 'OUTPUT : kb = ',ikb 938 939 WRITE(numout,*) '==================================' 940 WRITE(numout,*) 'Test icb_utl_zavg' 941 zD = 0.0 ; ze3= 20.0 ; zin=1.0 942 CALL icb_utl_getkb(ikb, ze3, zD) 943 CALL icb_utl_zavg(zout, zin, ze3, zD, ikb) 944 WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb 945 WRITE(numout,*) 'OUTPUT : zout = ',zout 946 947 zD = 50.0 ; ze3= 20.0 ; zin=1.0; zin(3:jpk) = 0.0 948 CALL icb_utl_getkb(ikb, ze3, zD) 949 CALL icb_utl_zavg(zout, zin, ze3, zD, ikb) 950 WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb 951 WRITE(numout,*) 'OUTPUT : zout = ',zout 952 CALL FLUSH(numout) 953 954 zD = 80.0 ; ze3= 20.0 ; zin=1.0; zin(3:jpk) = 0.0 955 CALL icb_utl_getkb(ikb, ze3, zD) 956 CALL icb_utl_zavg(zout, zin, ze3, zD, ikb) 957 WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb 958 WRITE(numout,*) 'OUTPUT : zout = ',zout 959 960 zD = 80 ; ze3= 20.0 ; zin=1.0 ; zin(3:jpk) = 0.0 961 CALL icb_utl_getkb(ikb, ze3, zD) 962 ikb = 2 963 CALL icb_utl_zavg(zout, zin, ze3, zD, ikb) 964 WRITE(numout,*) 'INPUT : zD = ',zD,' ze3 = ',ze3(1),' zin = ', zin, ' ikb = ',ikb 965 WRITE(numout,*) 'OUTPUT : zout = ',zout 966 967 CALL FLUSH(numout) 968 969 END SUBROUTINE test_icb_utl_getkb 970 825 971 !!====================================================================== 826 972 END MODULE icbutl
Note: See TracChangeset
for help on using the changeset viewer.
