Changeset 14031
- Timestamp:
- 2020-12-03T10:39:49+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source
- Files:
-
- 12 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/cfgs/SHARED/namelist_ref
r14015 r14031 629 629 ln_use_calving = .false. ! Use calving data even when nn_test_icebergs > 0 630 630 rn_speed_limit = 0. ! CFL speed limit for a berg 631 631 ! 632 ln_M2016 = .false. ! use Merino et al. (2016) modification (use of 3d ocean data instead of only sea surface data) 633 ln_icb_grd = .false. ! ground icb when icb bottom level hit oce bottom level (need ln_M2016 to be activated) 634 ! 632 635 cn_dir = './' ! root directory for the calving data location 633 636 !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icb_oce.F90
r13286 r14031 57 57 TYPE, PUBLIC :: point !: properties of an individual iceberg (position, mass, size, etc...) 58 58 INTEGER :: year 59 REAL(wp) :: xi , yj ! iceberg coordinates in the (i,j) referential (global)59 REAL(wp) :: xi , yj , zk ! iceberg coordinates in the (i,j) referential (global) and deepest level affected 60 60 REAL(wp) :: e1 , e2 ! horizontal scale factors at the iceberg position 61 61 REAL(wp) :: lon, lat, day ! geographic position 62 62 REAL(wp) :: mass, thickness, width, length, uvel, vvel ! iceberg physical properties 63 REAL(wp) :: uo, vo, ui, vi, ua, va, ssh_x, ssh_y, sst, cn, hi, sss! properties of iceberg environment63 REAL(wp) :: ssu, ssv, ui, vi, ua, va, ssh_x, ssh_y, sst, sss, cn, hi ! properties of iceberg environment 64 64 REAL(wp) :: mass_of_bits, heat_density 65 INTEGER :: kb ! icb bottom level 65 66 END TYPE point 66 67 … … 85 86 ! Extra arrays with bigger halo, needed when interpolating forcing onto iceberg position 86 87 ! 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_e88 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ff_e, tt_e, fr_e, ss_e88 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ssu_e, ssv_e 89 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: sst_e, sss_e, fr_e 89 90 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ua_e, va_e 90 91 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ssh_e 91 92 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: tmask_e, umask_e, vmask_e 93 REAl(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: rlon_e, rlat_e, ff_e 94 REAl(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: uoce_e, voce_e, toce_e, e3t_e 95 ! 92 96 #if defined key_si3 || defined key_cice 93 97 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hi_e, ui_e, vi_e … … 117 121 INTEGER , PUBLIC :: nn_verbose_write !: timesteps between verbose messages 118 122 REAL(wp), PUBLIC :: rn_rho_bergs !: Density of icebergs 123 REAL(wp), PUBLIC :: rho_berg_1_oce !: convertion factor (thickness to draft) (rn_rho_bergs/pp_rho_seawater) 119 124 REAL(wp), PUBLIC :: rn_LoW_ratio !: Initial ratio L/W for newly calved icebergs 120 125 REAL(wp), PUBLIC :: rn_bits_erosion_fraction !: Fraction of erosion melt flux to divert to bergy bits … … 124 129 LOGICAL , PUBLIC :: ln_time_average_weight !: Time average the weight on the ocean !!gm I don't understand that ! 125 130 REAL(wp), PUBLIC :: rn_speed_limit !: CFL speed limit for a berg 131 LOGICAL , PUBLIC :: ln_M2016, ln_icb_grd !: use Nacho's Merino 2016 work 126 132 ! 127 133 ! restart … … 135 141 REAL(wp), DIMENSION(nclasses), PUBLIC :: rn_initial_thickness ! Single instance of an icebergs type initialised in icebergs_init and updated in icebergs_run 136 142 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: src_calving, src_calving_hflx !: accumulate input ice 143 INTEGER , PUBLIC , SAVE :: micbkb !: deepest level affected by icebergs 137 144 INTEGER , PUBLIC , SAVE :: numicb !: iceberg IO 138 145 INTEGER , PUBLIC , SAVE, DIMENSION(nkounts) :: num_bergs !: iceberg counter … … 171 178 ! 172 179 ! expanded arrays for bilinear interpolation 173 ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,&174 & vo_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,&180 ALLOCATE( ssu_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) , & 181 & ssv_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) , & 175 182 #if defined key_si3 || defined key_cice 176 183 & ui_e(0:jpi+1,0:jpj+1) , & … … 178 185 & hi_e(0:jpi+1,0:jpj+1) , & 179 186 #endif 180 & f f_e(0:jpi+1,0:jpj+1) , fr_e(0:jpi+1,0:jpj+1) ,&181 & tt_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,&182 & ss _e(0:jpi+1,0:jpj+1) ,&187 & fr_e(0:jpi+1,0:jpj+1) , & 188 & sst_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) , & 189 & sss_e(0:jpi+1,0:jpj+1) , & 183 190 & first_width(nclasses) , first_length(nclasses) , & 184 191 & src_calving (jpi,jpj) , & … … 186 193 icb_alloc = icb_alloc + ill 187 194 195 IF ( ln_M2016 ) THEN 196 ALLOCATE( uoce_e(0:jpi+1,0:jpj+1,jpk), voce_e(0:jpi+1,0:jpj+1,jpk), & 197 & toce_e(0:jpi+1,0:jpj+1,jpk), e3t_e(0:jpi+1,0:jpj+1,jpk) , STAT=ill ) 198 icb_alloc = icb_alloc + ill 199 END IF 200 ! 188 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), & 189 & 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) 190 203 icb_alloc = icb_alloc + ill 191 204 -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbclv.F90
r13295 r14031 21 21 USE lbclnk ! NEMO boundary exchanges for gridded data 22 22 23 USE icb_oce ! iceberg variables24 23 USE icbdia ! iceberg diagnostics 25 24 USE icbutl ! iceberg utility routines … … 142 141 newpt%mass = rn_initial_mass (jn) 143 142 newpt%thickness = rn_initial_thickness(jn) 143 newpt%kb = 1 ! compute correctly in icbthm if needed 144 144 newpt%width = first_width (jn) 145 145 newpt%length = first_length (jn) … … 172 172 END DO 173 173 ! 174 DO jn = 1, nclasses175 CALL lbc_lnk( 'icbclv', berg_grid%stored_ice(:,:,jn), 'T', 1._wp )176 END DO177 CALL lbc_lnk( 'icbclv', berg_grid%stored_heat, 'T', 1._wp )178 !179 174 IF( nn_verbose_level > 0 .AND. icntmax > 1 ) WRITE(numicb,*) 'icb_clv: icnt=', icnt,' on', narea 180 175 ! -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbdyn.F90
r13281 r14031 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 … … 97 98 zyj2 = zyj1 + zdt_2 * zv1 ; zvvel2 = zvvel1 + zdt_2 * zay1 98 99 ! 99 CALL icb_ground( zxi2, zxi1, zu1, &100 & zyj2, zyj1, zv1, ll_bounced )100 CALL icb_ground( berg, zxi2, zxi1, zu1, & 101 & zyj2, zyj1, zv1, ll_bounced ) 101 102 102 103 ! !** A2 = A(X2,V2) … … 113 114 zyj3 = zyj1 + zdt_2 * zv2 ; zvvel3 = zvvel1 + zdt_2 * zay2 114 115 ! 115 CALL icb_ground( zxi3, zxi1, zu3, &116 & zyj3, zyj1, zv3, ll_bounced )116 CALL icb_ground( berg, zxi3, zxi1, zu3, & 117 & zyj3, zyj1, zv3, ll_bounced ) 117 118 118 119 ! !** A3 = A(X3,V3) … … 129 130 zyj4 = zyj1 + zdt * zv3 ; zvvel4 = zvvel1 + zdt * zay3 130 131 131 CALL icb_ground( zxi4, zxi1, zu4, &132 & zyj4, zyj1, zv4, ll_bounced )132 CALL icb_ground( berg, zxi4, zxi1, zu4, & 133 & zyj4, zyj1, zv4, ll_bounced ) 133 134 134 135 ! !** A4 = A(X4,V4) … … 148 149 zvvel_n = pt%vvel + zdt_6 * ( zay1 + 2.*(zay2 + zay3) + zay4 ) 149 150 150 CALL icb_ground( zxi_n, zxi1, zuvel_n, &151 & zyj_n, zyj1, zvvel_n, ll_bounced )151 CALL icb_ground( berg, zxi_n, zxi1, zuvel_n, & 152 & zyj_n, zyj1, zvvel_n, ll_bounced ) 152 153 153 154 pt%uvel = zuvel_n !** save in berg structure … … 156 157 pt%yj = zyj_n 157 158 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 162 159 berg => berg%next ! switch to the next berg 163 160 ! … … 167 164 168 165 169 SUBROUTINE icb_ground( pi, pi0, pu, &170 & pj, pj0, pv, ld_bounced )166 SUBROUTINE icb_ground( berg, pi, pi0, pu, & 167 & pj, pj0, pv, ld_bounced ) 171 168 !!---------------------------------------------------------------------- 172 169 !! *** ROUTINE icb_ground *** … … 177 174 !! NB two possibilities available one of which is hard-coded here 178 175 !!---------------------------------------------------------------------- 176 TYPE(iceberg ), POINTER, INTENT(in ) :: berg ! berg 177 ! 179 178 REAL(wp), INTENT(inout) :: pi , pj ! current iceberg position 180 179 REAL(wp), INTENT(in ) :: pi0, pj0 ! previous iceberg position … … 184 183 INTEGER :: ii, ii0 185 184 INTEGER :: ij, ij0 185 INTEGER :: ikb 186 186 INTEGER :: ibounce_method 187 ! 188 REAL(wp) :: zD 189 REAL(wp), DIMENSION(jpk) :: ze3t 187 190 !!---------------------------------------------------------------------- 188 191 ! … … 200 203 ij = mj1( ij ) 201 204 ! 202 IF( tmask(ii,ij,1) /= 0._wp ) RETURN ! berg reach a new t-cell, but an ocean one 205 ! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0 206 IF ( ln_M2016 .AND. ln_icb_grd ) THEN 207 ! 208 ! draught (keel depth) 209 zD = rho_berg_1_oce * berg%current_point%thickness 210 ! 211 ! interpol needed data 212 CALL icb_utl_interp( pi, pj, pe3t=ze3t ) 213 ! 214 !compute bottom level 215 CALL icb_utl_getkb( ikb, ze3t, zD ) 216 ! 217 ! berg reach a new t-cell, but an ocean one 218 ! .AND. needed in case berg hit an isf (tmask(ii,ij,1) == 0 and tmask(ii,ij,ikb) /= 0) 219 IF( tmask(ii,ij,ikb) /= 0._wp .AND. tmask(ii,ij,1) /= 0._wp ) RETURN 220 ! 221 ELSE 222 IF( tmask(ii,ij,1) /= 0._wp ) RETURN ! berg reach a new t-cell, but an ocean one 223 END IF 203 224 ! 204 225 ! From here, berg have reach land: treat grounding/bouncing … … 257 278 REAL(wp), PARAMETER :: pp_Cr0 = 0.06_wp ! 258 279 ! 259 INTEGER :: itloop 260 REAL(wp) :: zuo, z ui, zua, zuwave, zssh_x, zsst, zcn, zhi, zsss261 REAL(wp) :: zvo, z vi, zva, zvwave, zssh_y280 INTEGER :: itloop, ikb, jk 281 REAL(wp) :: zuo, zssu, zui, zua, zuwave, zssh_x, zcn, zhi 282 REAL(wp) :: zvo, zssv, zvi, zva, zvwave, zssh_y 262 283 REAL(wp) :: zff, zT, zD, zW, zL, zM, zF 263 284 REAL(wp) :: zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad 264 REAL(wp) :: z_ocn, z_atm, z_ice 285 REAL(wp) :: z_ocn, z_atm, z_ice, zdep 265 286 REAL(wp) :: zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop 266 287 REAL(wp) :: zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi 267 288 REAL(wp) :: zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new 289 REAL(wp), DIMENSION(jpk) :: zuoce, zvoce, ze3t, zdepw 268 290 !!---------------------------------------------------------------------- 269 291 270 292 ! Interpolate gridded fields to berg 271 293 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, zsss ) 294 CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2, & ! scale factor 295 & pssu=zssu, pui=zui, pua=zua, & ! oce/ice/atm velocities 296 & pssv=zssv, pvi=zvi, pva=zva, & ! oce/ice/atm velocities 297 & pssh_i=zssh_x, pssh_j=zssh_y, & ! ssh gradient 298 & phi=zhi, pff=zff) ! ice thickness and coriolis 274 299 275 300 zM = berg%current_point%mass 276 301 zT = berg%current_point%thickness ! total thickness 277 zD = ( rn_rho_bergs / pp_rho_seawater ) * zT! draught (keel depth)302 zD = rho_berg_1_oce * zT ! draught (keel depth) 278 303 zF = zT - zD ! freeboard 279 304 zW = berg%current_point%width … … 282 307 zhi = MIN( zhi , zD ) 283 308 zD_hi = MAX( 0._wp, zD-zhi ) 284 285 286 zuwave = zua - z uo ; zvwave = zva - zvo! Use wind speed rel. to ocean for wave model309 310 ! Wave radiation 311 zuwave = zua - zssu ; zvwave = zva - zssv ! Use wind speed rel. to ocean for wave model 287 312 zwmod = zuwave*zuwave + zvwave*zvwave ! The wave amplitude and length depend on the current; 288 313 ! ! wind speed relative to the ocean. Actually wmod is wmod**2 here. … … 309 334 IF( abs(zui) + abs(zvi) == 0._wp ) z_ice = 0._wp 310 335 336 ! lateral velocities 337 ! default ssu and ssv 338 ! ln_M2016: mean velocity along the profile 339 IF ( ln_M2016 ) THEN 340 ! interpol needed data 341 CALL icb_utl_interp( pxi, pyj, puoce=zuoce, pvoce=zvoce, pe3t=ze3t ) ! 3d velocities 342 343 !compute bottom level 344 CALL icb_utl_getkb( ikb, ze3t, zD ) 345 346 ! compute mean velocity 347 CALL icb_utl_zavg(zuo, zuoce, ze3t, zD, ikb) 348 CALL icb_utl_zavg(zvo, zvoce, ze3t, zD, ikb) 349 ELSE 350 zuo = zssu 351 zvo = zssv 352 END IF 353 311 354 zuveln = puvel ; zvveln = pvvel ! Copy starting uvel, vvel 312 355 ! … … 321 364 ! Explicit accelerations 322 365 !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave & 323 ! -zdrag_ocn*(puvel-z uo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)366 ! -zdrag_ocn*(puvel-zssu) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) 324 367 !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave & 325 ! -zdrag_ocn*(pvvel-z vo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)368 ! -zdrag_ocn*(pvvel-zssv) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) 326 369 zaxe = -grav * zssh_x + zwave_rad * zuwave 327 370 zaye = -grav * zssh_y + zwave_rad * zvwave -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbini.F90
r13295 r14031 73 73 ! 74 74 IF( .NOT. ln_icebergs ) RETURN 75 75 ! 76 76 ! ! allocate gridded fields 77 77 IF( icb_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'icb_alloc : unable to allocate arrays' ) 78 78 ! 79 79 ! ! initialised variable with extra haloes to zero 80 uo_e(:,:) = 0._wp ; vo_e(:,:) = 0._wp ; 81 ua_e(:,:) = 0._wp ; va_e(:,:) = 0._wp ; 82 ff_e(:,:) = 0._wp ; tt_e(:,:) = 0._wp ; 83 fr_e(:,:) = 0._wp ; ss_e(:,:) = 0._wp ; 80 ssu_e(:,:) = 0._wp ; ssv_e(:,:) = 0._wp ; 81 ua_e(:,:) = 0._wp ; va_e(:,:) = 0._wp ; 82 ff_e(:,:) = 0._wp ; sst_e(:,:) = 0._wp ; 83 fr_e(:,:) = 0._wp ; sss_e(:,:) = 0._wp ; 84 ! 85 IF ( ln_M2016 ) THEN 86 toce_e(:,:,:) = 0._wp 87 uoce_e(:,:,:) = 0._wp 88 voce_e(:,:,:) = 0._wp 89 e3t_e(:,:,:) = 0._wp 90 END IF 91 ! 84 92 #if defined key_si3 85 93 hi_e(:,:) = 0._wp ; … … 100 108 first_width (:) = SQRT( rn_initial_mass(:) / ( rn_LoW_ratio * rn_rho_bergs * rn_initial_thickness(:) ) ) 101 109 first_length(:) = rn_LoW_ratio * first_width(:) 110 rho_berg_1_oce = rn_rho_bergs / pp_rho_seawater ! scale factor used for convertion thickness to draft 111 ! 112 ! deepest level affected by icebergs 113 ! can be tuned but the safest is this 114 ! (with z* and z~ the depth of each level change overtime, so the more robust micbkb is jpk) 115 micbkb = jpk 102 116 103 117 berg_grid%calving (:,:) = 0._wp … … 240 254 vmask_e(:,:) = 0._wp ; vmask_e(1:jpi,1:jpj) = vmask(:,:,1) 241 255 CALL lbc_lnk_icb( 'icbini', tmask_e, 'T', +1._wp, 1, 1 ) 242 CALL lbc_lnk_icb( 'icbini', umask_e, 'T', +1._wp, 1, 1 ) 243 CALL lbc_lnk_icb( 'icbini', vmask_e, 'T', +1._wp, 1, 1 ) 244 ! 256 CALL lbc_lnk_icb( 'icbini', umask_e, 'U', +1._wp, 1, 1 ) 257 CALL lbc_lnk_icb( 'icbini', vmask_e, 'V', +1._wp, 1, 1 ) 258 259 ! definition of extended lat/lon array needed by icb_bilin_h 260 rlon_e(:,:) = 0._wp ; rlon_e(1:jpi,1:jpj) = glamt(:,:) 261 rlat_e(:,:) = 0._wp ; rlat_e(1:jpi,1:jpj) = gphit(:,:) 262 CALL lbc_lnk_icb( 'icbini', rlon_e, 'T', +1._wp, 1, 1 ) 263 CALL lbc_lnk_icb( 'icbini', rlat_e, 'T', +1._wp, 1, 1 ) 264 ! 265 ! definnitionn of extennded ff_f array needed by icb_utl_interp 266 ff_e(:,:) = 0._wp ; ff_e(1:jpi,1:jpj) = ff_f(:,:) 267 CALL lbc_lnk_icb( 'icbini', ff_e, 'F', +1._wp, 1, 1 ) 268 245 269 ! assign each new iceberg with a unique number constructed from the processor number 246 270 ! and incremented by the total number of processors … … 338 362 localpt%xi = REAL( mig(ji), wp ) 339 363 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' ) 364 CALL icb_utl_interp( localpt%xi, localpt%yj, plat=localpt%lat, plon=localpt%lon ) 342 365 localpt%mass = rn_initial_mass (iberg) 343 366 localpt%thickness = rn_initial_thickness(iberg) … … 350 373 localpt%uvel = 0._wp 351 374 localpt%vvel = 0._wp 375 localpt%kb = 1 352 376 CALL icb_utl_incr() 353 377 localberg%number(:) = num_bergs(:) … … 383 407 & rn_bits_erosion_fraction , rn_sicn_shift , ln_passive_mode , & 384 408 & ln_time_average_weight , nn_test_icebergs , rn_test_box , & 385 & ln_use_calving , rn_speed_limit , cn_dir, sn_icb , & 386 & cn_icbrst_indir, cn_icbrst_in , cn_icbrst_outdir , cn_icbrst_out 409 & ln_use_calving , rn_speed_limit , cn_dir, sn_icb , ln_M2016 , & 410 & cn_icbrst_indir, cn_icbrst_in , cn_icbrst_outdir , cn_icbrst_out , & 411 & ln_icb_grd 387 412 !!---------------------------------------------------------------------- 388 413 … … 463 488 & 'bits_erosion_fraction = ', rn_bits_erosion_fraction 464 489 490 WRITE(numout,*) ' Use icb module modification from Merino et al. (2016) : ln_M2016 = ', ln_M2016 491 WRITE(numout,*) ' ground icebergs if icb bottom lvl hit the oce bottom level : ln_icb_grd = ', ln_icb_grd 492 465 493 WRITE(numout,*) ' Shift of sea-ice concentration in erosion flux modulation ', & 466 494 & '(0<sicn_shift<1) rn_sicn_shift = ', rn_sicn_shift -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbstp.F90
r11536 r14031 52 52 CONTAINS 53 53 54 SUBROUTINE icb_stp( kt )54 SUBROUTINE icb_stp( kt, Kmm ) 55 55 !!---------------------------------------------------------------------- 56 56 !! *** ROUTINE icb_stp *** … … 61 61 !!---------------------------------------------------------------------- 62 62 INTEGER, INTENT(in) :: kt ! time step index 63 INTEGER, INTENT(in) :: Kmm ! ocean time level index 63 64 ! 64 65 LOGICAL :: ll_sample_traj, ll_budget, ll_verbose ! local logical … … 70 71 ! 71 72 nktberg = kt 73 ! 74 !CALL test_icb_utl_getkb 75 !CALL ctl_stop('end test icb') 72 76 ! 73 77 IF( nn_test_icebergs < 0 .OR. ln_use_calving ) THEN !* read calving data … … 92 96 ! 93 97 ! !* copy nemo forcing arrays into iceberg versions with extra halo 94 CALL icb_utl_copy( ) ! only necessary for variables not on T points98 CALL icb_utl_copy( Kmm ) ! only necessary for variables not on T points 95 99 ! 96 100 ! -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbthm.F90
r13281 r14031 49 49 INTEGER, INTENT(in) :: kt ! timestep number, just passed to icb_utl_print_berg 50 50 ! 51 INTEGER :: ii, ij 52 REAL(wp) :: zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn 51 INTEGER :: ii, ij, jk, ikb 52 REAL(wp) :: zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn, zD, zvb, zub, ztb 53 REAL(wp) :: zMv, zMe, zMb, zmelt, zdvo, zdvob, zdva, zdM, zSs, zdMe, zdMb, zdMv 53 54 REAL(wp) :: zSSS, zfzpt 54 REAL(wp) :: zMv, zMe, zMb, zmelt, zdvo, zdva, zdM, zSs, zdMe, zdMb, zdMv55 55 REAL(wp) :: zMnew, zMnew1, zMnew2, zheat_hcflux, zheat_latent, z1_12 56 56 REAL(wp) :: zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb 57 REAL(wp) :: zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2 57 REAL(wp) :: zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2, zdepw 58 REAL(wp), DIMENSION(jpk) :: ztoce, zuoce, zvoce, ze3t, zzMv 58 59 TYPE(iceberg), POINTER :: this, next 59 60 TYPE(point) , POINTER :: pt … … 85 86 pt => this%current_point 86 87 nknberg = this%number(1) 87 CALL icb_utl_interp( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x, & 88 & pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y, & 89 & pt%sst, pt%cn, pt%hi, zff, pt%sss ) 88 89 CALL icb_utl_interp( pt%xi, pt%yj, & ! position 90 & pssu=pt%ssu, pua=pt%ua, & ! oce/atm velocities 91 & pssv=pt%ssv, pva=pt%va, & ! oce/atm velocities 92 & psst=pt%sst, pcn=pt%cn, & 93 & psss=pt%sss ) 94 95 IF ( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) THEN 96 CALL icb_utl_interp( pt%xi, pt%yj, pe1=pt%e1, pe2=pt%e2, & 97 & pui=pt%ui, pssh_i=pt%ssh_x, & 98 & pvi=pt%vi, pssh_j=pt%ssh_y, & 99 & phi=pt%hi, & 100 & plat=pt%lat, plon=pt%lon ) 101 END IF 90 102 ! 91 103 zSST = pt%sst … … 95 107 zM = pt%mass 96 108 zT = pt%thickness ! total thickness 97 ! D = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth) 98 ! F = zT - D ! freeboard 109 zD = rho_berg_1_oce * zT ! draught (keel depth) 99 110 zW = pt%width 100 111 zL = pt%length … … 108 119 109 120 ! Environment 110 zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 ) 111 zdva = SQRT( (pt%ua -pt%uo)**2 + (pt%va -pt%vo)**2 ) 112 zSs = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva ! Sea state (eqn M.A9) 113 121 ! default sst, ssu and ssv 122 ! ln_M2016: use temp, u and v profile 123 IF ( ln_M2016 ) THEN 124 125 ! load t, u, v and e3 profile at icb position 126 CALL icb_utl_interp( pt%xi, pt%yj, ptoce=ztoce, puoce=zuoce, pvoce=zvoce, pe3t=ze3t ) 127 128 !compute bottom level 129 CALL icb_utl_getkb( pt%kb, ze3t, zD ) 130 131 ikb = MIN(pt%kb,mbkt(ii,ij)) ! limit pt%kb by mbkt 132 ! => bottom temperature used to fill ztoce(mbkt:jpk) 133 ztb = ztoce(ikb) ! basal temperature 134 zub = zuoce(ikb) 135 zvb = zvoce(ikb) 136 ELSE 137 ztb = pt%sst 138 zub = pt%ssu 139 zvb = pt%ssv 140 END IF 141 142 zdvob = SQRT( (pt%uvel-zub)**2 + (pt%vvel-zvb)**2 ) ! relative basal velocity 143 zdva = SQRT( (pt%ua -pt%ssu)**2 + (pt%va -pt%ssv)**2 ) ! relative wind 144 zSs = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva ! Sea state (eqn M.A9) 145 ! 114 146 ! Melt rates in m/s (i.e. division by rday) 115 zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2) , 0._wp ) * z1_rday ! Buoyant convection at sides (eqn M.A10) 147 ! Buoyant convection at sides (eqn M.A10) 148 IF ( ln_M2016 ) THEN 149 ! averaging along all the iceberg draft 150 zzMv(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday 151 CALL icb_utl_zavg(zMv, zzMv, ze3t, zD, ikb ) 152 ELSE 153 zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday 154 END IF 155 ! 156 ! Basal turbulent melting (eqn M.A7 ) 116 157 IF ( zSST > zfzpt ) THEN ! Calculate basal melting only if SST above freezing point 117 zMb = MAX( 0.58_wp*(zdvo **0.8_wp)*(zSST+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday ! Basal turbulent melting (eqn M.A7 )158 zMb = MAX( 0.58_wp*(zdvob**0.8_wp)*(ztb+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday 118 159 ELSE 119 160 zMb = 0._wp ! No basal melting if SST below freezing point 120 161 ENDIF 121 zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3))) , 0._wp ) * z1_rday ! Wave erosion (eqn M.A8 ) 162 ! 163 ! Wave erosion (eqn M.A8 ) 164 zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3))) , 0._wp ) * z1_rday 122 165 123 166 IF( ln_operator_splitting ) THEN ! Operator split update of volume/mass … … 207 250 208 251 ! Rolling 209 zDn = ( rn_rho_bergs / pp_rho_seawater )* zTn ! draught (keel depth)252 zDn = rho_berg_1_oce * zTn ! draught (keel depth) 210 253 IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN 211 254 zT = zTn … … 224 267 225 268 !!gm add a test to avoid over melting ? 269 !!pm I agree, over melting could break conservation (more melt than calving) 226 270 227 271 IF( zMnew <= 0._wp ) THEN ! Delete the berg if completely melted -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbtrj.F90
r13886 r14031 40 40 INTEGER :: numberid, nstepid, nscaling_id 41 41 INTEGER :: nlonid, nlatid, nxid, nyid, nuvelid, nvvelid, nmassid 42 INTEGER :: n uoid, nvoid, nuaid, nvaid, nuiid, nviid42 INTEGER :: nssuid, nssvid, nuaid, nvaid, nuiid, nviid 43 43 INTEGER :: nsshxid, nsshyid, nsstid, ncntid, nthkid 44 44 INTEGER :: nthicknessid, nwidthid, nlengthid … … 111 111 iret = NF90_DEF_VAR( ntrajid, 'uvel' , NF90_DOUBLE, n_dim , nuvelid ) 112 112 iret = NF90_DEF_VAR( ntrajid, 'vvel' , NF90_DOUBLE, n_dim , nvvelid ) 113 iret = NF90_DEF_VAR( ntrajid, ' uto' , NF90_DOUBLE, n_dim , nuoid)114 iret = NF90_DEF_VAR( ntrajid, ' vto' , NF90_DOUBLE, n_dim , nvoid)113 iret = NF90_DEF_VAR( ntrajid, 'ssu' , NF90_DOUBLE, n_dim , nssuid ) 114 iret = NF90_DEF_VAR( ntrajid, 'ssv' , NF90_DOUBLE, n_dim , nssvid ) 115 115 iret = NF90_DEF_VAR( ntrajid, 'uta' , NF90_DOUBLE, n_dim , nuaid ) 116 116 iret = NF90_DEF_VAR( ntrajid, 'vta' , NF90_DOUBLE, n_dim , nvaid ) … … 148 148 iret = NF90_PUT_ATT( ntrajid, nvvelid , 'long_name', 'meridional velocity' ) 149 149 iret = NF90_PUT_ATT( ntrajid, nvvelid , 'units' , 'm/s' ) 150 iret = NF90_PUT_ATT( ntrajid, n uoid, 'long_name', 'ocean u component' )151 iret = NF90_PUT_ATT( ntrajid, n uoid, 'units' , 'm/s' )152 iret = NF90_PUT_ATT( ntrajid, n void, 'long_name', 'ocean v component' )153 iret = NF90_PUT_ATT( ntrajid, n void, 'units' , 'm/s' )150 iret = NF90_PUT_ATT( ntrajid, nssuid , 'long_name', 'ocean u component' ) 151 iret = NF90_PUT_ATT( ntrajid, nssuid , 'units' , 'm/s' ) 152 iret = NF90_PUT_ATT( ntrajid, nssvid , 'long_name', 'ocean v component' ) 153 iret = NF90_PUT_ATT( ntrajid, nssvid , 'units' , 'm/s' ) 154 154 iret = NF90_PUT_ATT( ntrajid, nuaid , 'long_name', 'atmosphere u component' ) 155 155 iret = NF90_PUT_ATT( ntrajid, nuaid , 'units' , 'm/s' ) … … 231 231 iret = NF90_PUT_VAR( ntrajid, nuvelid , pt%uvel , (/ jn /) ) 232 232 iret = NF90_PUT_VAR( ntrajid, nvvelid , pt%vvel , (/ jn /) ) 233 iret = NF90_PUT_VAR( ntrajid, n uoid , pt%uo, (/ jn /) )234 iret = NF90_PUT_VAR( ntrajid, n void , pt%vo, (/ jn /) )233 iret = NF90_PUT_VAR( ntrajid, nssuid , pt%ssu , (/ jn /) ) 234 iret = NF90_PUT_VAR( ntrajid, nssvid , pt%ssv , (/ jn /) ) 235 235 iret = NF90_PUT_VAR( ntrajid, nuaid , pt%ua , (/ jn /) ) 236 236 iret = NF90_PUT_VAR( ntrajid, nvaid , pt%va , (/ jn /) ) -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/ICB/icbutl.F90
r13281 r14031 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 -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/OCE/SBC/sbcmod.F90
r14015 r14031 475 475 476 476 IF( ln_icebergs ) THEN 477 CALL icb_stp( kt ) ! compute icebergs477 CALL icb_stp( kt, Kmm ) ! compute icebergs 478 478 ! Icebergs do not melt over the haloes. 479 479 ! So emp values over the haloes are no more consistent with the inner domain values. -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/SAS/nemogcm.F90
r14015 r14031 34 34 USE diu_layers ! diurnal bulk SST and coolskin 35 35 USE step_diu ! diurnal bulk SST timestepping (called from here if run offline) 36 USE icb_oce ! icebergs 36 37 ! 37 38 USE prtctl ! Print control … … 393 394 ! ==> 394 395 CALL icb_init( rn_Dt, nit000) ! initialise icebergs instance 396 397 ! compatibility check 398 IF( ln_icebergs .AND. ln_M2016 ) THEN 399 IF( lwp ) WRITE(numout,*) ' ==>>> ln_iceberg and ln_M2016 not compatible with SAS (need 3d data)' 400 CALL ctl_stop('ln_iceberg and ln_M2016 not compatible with SAS (need 3d data)') 401 END IF 395 402 ! 396 403 IF(lwp) WRITE(numout,cform_aaa) ! Flag AAAAAAA -
NEMO/branches/2020/dev_r13541_TOP-01_rlod_Antarctic_ice_Sheet_Fe_Source/src/SAS/sbcssm.F90
r13286 r14031 21 21 USE zpshde ! z-coord. with partial steps: horizontal derivatives 22 22 USE closea ! for ln_closea 23 USE icb_oce ! for icebergs 23 24 ! 24 25 USE in_out_manager ! I/O manager … … 226 227 ln_closea = .false. 227 228 ENDIF 228 229 IF( ln_icebergs .AND. ln_M2016 ) THEN 230 IF( lwp ) WRITE(numout,*) ' ==>>> ln_iceberg and ln_M2016 not compatible with SAS (need 3d data)' 231 CALL ctl_stop('ln_iceberg and ln_M2016 not compatible with SAS (need 3d data)') 232 END IF 229 233 ! 230 234 IF( l_sasread ) THEN ! store namelist information in an array
Note: See TracChangeset
for help on using the changeset viewer.