Changeset 13357 for NEMO/branches
- Timestamp:
- 2020-07-30T13:20:57+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/tickets_icb_1900
- Files:
-
- 10 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/namelist_ref
r13216 r13357 614 614 ln_use_calving = .false. ! Use calving data even when nn_test_icebergs > 0 615 615 rn_speed_limit = 0. ! CFL speed limit for a berg 616 616 ! 617 ln_M2016 = .false. ! use Merino et al. (2016) modification (use of 3d ocean data instead of only sea surface data) 618 ! 617 619 cn_dir = './' ! root directory for the calving data location 618 620 !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! -
NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icb_oce.F90
r13265 r13357 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! properties of iceberg environment63 REAL(wp) :: ssu, ssv, ui, vi, ua, va, ssh_x, ssh_y, sst, 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 :: tt_e, fr_e88 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: ssu_e, ssv_e 89 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: sst_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 92 93 REAl(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: rlon_e, rlat_e, ff_e 93 ! prepration to Nacho work 94 ! tt3d_e 95 ! uo3d_e 96 ! vo3d_e 97 ! 94 REAl(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE :: uoce_e, voce_e, toce_e, e3t_e 95 ! 98 96 #if defined key_si3 || defined key_cice 99 97 REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: hi_e, ui_e, vi_e … … 123 121 INTEGER , PUBLIC :: nn_verbose_write !: timesteps between verbose messages 124 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) 125 124 REAL(wp), PUBLIC :: rn_LoW_ratio !: Initial ratio L/W for newly calved icebergs 126 125 REAL(wp), PUBLIC :: rn_bits_erosion_fraction !: Fraction of erosion melt flux to divert to bergy bits … … 130 129 LOGICAL , PUBLIC :: ln_time_average_weight !: Time average the weight on the ocean !!gm I don't understand that ! 131 130 REAL(wp), PUBLIC :: rn_speed_limit !: CFL speed limit for a berg 131 LOGICAL , PUBLIC :: ln_M2016 !: use Nacho's Merino 2016 work 132 132 ! 133 133 ! restart … … 141 141 REAL(wp), DIMENSION(nclasses), PUBLIC :: rn_initial_thickness ! Single instance of an icebergs type initialised in icebergs_init and updated in icebergs_run 142 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 143 144 INTEGER , PUBLIC , SAVE :: numicb !: iceberg IO 144 145 INTEGER , PUBLIC , SAVE, DIMENSION(nkounts) :: num_bergs !: iceberg counter … … 180 181 ! 181 182 ! expanded arrays for bilinear interpolation 182 ALLOCATE( uo_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) ,&183 & vo_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) ,&183 ALLOCATE( ssu_e(0:jpi+1,0:jpj+1) , ua_e(0:jpi+1,0:jpj+1) , & 184 & ssv_e(0:jpi+1,0:jpj+1) , va_e(0:jpi+1,0:jpj+1) , & 184 185 #if defined key_si3 || defined key_cice 185 186 & ui_e(0:jpi+1,0:jpj+1) , & … … 188 189 #endif 189 190 & fr_e(0:jpi+1,0:jpj+1) , & 190 & tt_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) ,&191 & sst_e(0:jpi+1,0:jpj+1) , ssh_e(0:jpi+1,0:jpj+1) , & 191 192 & first_width(nclasses) , first_length(nclasses) , & 192 193 & src_calving (jpi,jpj) , & … … 194 195 icb_alloc = icb_alloc + ill 195 196 196 ! prepration to Nacho work 197 ! tt3d_e 198 ! uo3d_e 199 ! vo3d_e 197 IF ( ln_M2016 ) THEN 198 ALLOCATE( uoce_e(0:jpi+1,0:jpj+1,jpk), voce_e(0:jpi+1,0:jpj+1,jpk), & 199 & toce_e(0:jpi+1,0:jpj+1,jpk), e3t_e(0:jpi+1,0:jpj+1,jpk) , STAT=ill ) 200 icb_alloc = icb_alloc + ill 201 END IF 200 202 ! 201 203 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), & -
NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbclv.F90
r13273 r13357 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) -
NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbdyn.F90
r13265 r13357 197 197 ij = mj1( ij ) 198 198 ! 199 ! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0 200 !IF ( ln_icb_ground ) THEN 201 ! ! interpol needed data 202 ! CALL icb_utl_interp( pxi, pyj, pe3t=ze3t ) ! 3d velocities 203 ! 204 ! !compute bottom level 205 ! CALL icb_utl_getkb( ikb, ze3t, zD ) 206 ! 207 ! IF( tmask(ii,ij,ikb) /= 0._wp ) RETURN ! berg reach a new t-cell, but an ocean one 208 !ELSE 199 209 IF( tmask(ii,ij,1) /= 0._wp ) RETURN ! berg reach a new t-cell, but an ocean one 210 !END IF 200 211 ! 201 212 ! From here, berg have reach land: treat grounding/bouncing … … 254 265 REAL(wp), PARAMETER :: pp_Cr0 = 0.06_wp ! 255 266 ! 256 INTEGER :: itloop 257 REAL(wp) :: zuo, z ui, zua, zuwave, zssh_x, zsst, zcn, zhi258 REAL(wp) :: zvo, z vi, zva, zvwave, zssh_y267 INTEGER :: itloop, ikb, jk 268 REAL(wp) :: zuo, zssu, zui, zua, zuwave, zssh_x, zcn, zhi 269 REAL(wp) :: zvo, zssv, zvi, zva, zvwave, zssh_y 259 270 REAL(wp) :: zff, zT, zD, zW, zL, zM, zF 260 271 REAL(wp) :: zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad 261 REAL(wp) :: z_ocn, z_atm, z_ice 272 REAL(wp) :: z_ocn, z_atm, z_ice, zdep 262 273 REAL(wp) :: zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop 263 274 REAL(wp) :: zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi 264 275 REAL(wp) :: zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new 276 REAL(wp), DIMENSION(jpk) :: zuoce, zvoce, ze3t, zdepw 265 277 !!---------------------------------------------------------------------- 266 278 267 279 ! Interpolate gridded fields to berg 268 280 nknberg = berg%number(1) 269 CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2, & ! scale factor 270 & puo=zuo, pui=zui, pua=zua, & ! oce/ice/atm velocities 271 & pvo=zvo, pvi=zvi, pva=zva, & ! oce/ice/atm velocities 272 & pssh_i=zssh_x, pssh_j=zssh_y, & ! ssh gradient 273 & phi=zhi, pff=zff ) ! ice thickness and coriolis 274 275 IF (lwp) THEN 276 WRITE(numout,*) 'icbdyn ', pxi, pyj 277 WRITE(numout,*) pe1, zuo, zui, zua, zssh_x 278 WRITE(numout,*) pe2, zvo, zvi, zva, zssh_y 279 WRITE(numout,*) zhi, zff 280 WRITE(numout,*) '' 281 END IF 282 281 CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2, & ! scale factor 282 & pssu=zssu, pui=zui, pua=zua, & ! oce/ice/atm velocities 283 & pssv=zssv, pvi=zvi, pva=zva, & ! oce/ice/atm velocities 284 & pssh_i=zssh_x, pssh_j=zssh_y, & ! ssh gradient 285 & phi=zhi, pff=zff) ! ice thickness and coriolis 283 286 284 287 zM = berg%current_point%mass 285 288 zT = berg%current_point%thickness ! total thickness 286 zD = ( rn_rho_bergs / pp_rho_seawater ) * zT! draught (keel depth)289 zD = rho_berg_1_oce * zT ! draught (keel depth) 287 290 zF = zT - zD ! freeboard 288 291 zW = berg%current_point%width … … 291 294 zhi = MIN( zhi , zD ) 292 295 zD_hi = MAX( 0._wp, zD-zhi ) 293 294 295 zuwave = zua - z uo ; zvwave = zva - zvo! Use wind speed rel. to ocean for wave model296 297 ! Wave radiation 298 zuwave = zua - zssu ; zvwave = zva - zssv ! Use wind speed rel. to ocean for wave model 296 299 zwmod = zuwave*zuwave + zvwave*zvwave ! The wave amplitude and length depend on the current; 297 300 ! ! wind speed relative to the ocean. Actually wmod is wmod**2 here. … … 318 321 IF( abs(zui) + abs(zvi) == 0._wp ) z_ice = 0._wp 319 322 323 ! lateral velocities 324 ! default ssu and ssv 325 ! ln_M2016: mean velocity along the profile 326 IF ( ln_M2016 ) THEN 327 ! interpol needed data 328 CALL icb_utl_interp( pxi, pyj, puoce=zuoce, pvoce=zvoce, pe3t=ze3t ) ! 3d velocities 329 330 !compute bottom level 331 CALL icb_utl_getkb( ikb, ze3t, zD ) 332 333 ! compute mean velocity 334 zuo = 0.0 ; zvo = 0.0; zdep = 0.0 335 DO jk = 1, ikb-1 336 zuo = zuo + zuoce(jk) * ze3t(jk) 337 zvo = zvo + zvoce(jk) * ze3t(jk) 338 zdep = zdep + ze3t(jk) 339 END DO 340 zuo = (zuo + zuoce(ikb) * (zD - zdep)) / zD 341 zvo = (zvo + zvoce(ikb) * (zD - zdep)) / zD 342 ELSE 343 zuo = zssu 344 zvo = zssv 345 END IF 346 320 347 zuveln = puvel ; zvveln = pvvel ! Copy starting uvel, vvel 321 348 ! … … 330 357 ! Explicit accelerations 331 358 !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave & 332 ! -zdrag_ocn*(puvel-z uo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)359 ! -zdrag_ocn*(puvel-zssu) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) 333 360 !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave & 334 ! -zdrag_ocn*(pvvel-z vo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)361 ! -zdrag_ocn*(pvvel-zssv) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) 335 362 zaxe = -grav * zssh_x + zwave_rad * zuwave 336 363 zaye = -grav * zssh_y + zwave_rad * zvwave -
NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbini.F90
r13276 r13357 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 ; 84 ! prepration to Nacho work 85 ! e3t_e 86 ! tt3d_e 87 ! uo3d_e 88 ! vo3d_e 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 ; 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 89 91 ! 90 92 #if defined key_si3 … … 106 108 first_width (:) = SQRT( rn_initial_mass(:) / ( rn_LoW_ratio * rn_rho_bergs * rn_initial_thickness(:) ) ) 107 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 108 116 109 117 berg_grid%calving (:,:) = 0._wp … … 365 373 localpt%uvel = 0._wp 366 374 localpt%vvel = 0._wp 375 localpt%kb = 1 367 376 CALL icb_utl_incr() 368 377 localberg%number(:) = num_bergs(:) … … 398 407 & rn_bits_erosion_fraction , rn_sicn_shift , ln_passive_mode , & 399 408 & ln_time_average_weight , nn_test_icebergs , rn_test_box , & 400 & ln_use_calving , rn_speed_limit , cn_dir, sn_icb , 409 & ln_use_calving , rn_speed_limit , cn_dir, sn_icb , ln_M2016 , & 401 410 & cn_icbrst_indir, cn_icbrst_in , cn_icbrst_outdir , cn_icbrst_out 402 411 !!---------------------------------------------------------------------- -
NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbstp.F90
r11536 r13357 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 … … 92 93 ! 93 94 ! !* copy nemo forcing arrays into iceberg versions with extra halo 94 CALL icb_utl_copy( ) ! only necessary for variables not on T points95 CALL icb_utl_copy( Kmm ) ! only necessary for variables not on T points 95 96 ! 96 97 ! -
NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbthm.F90
r13265 r13357 48 48 INTEGER, INTENT(in) :: kt ! timestep number, just passed to icb_utl_print_berg 49 49 ! 50 INTEGER :: ii, ij 51 REAL(wp) :: zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn 52 REAL(wp) :: zMv, zMe, zMb, zmelt, zdvo, zdv a, zdM, zSs, zdMe, zdMb, zdMv50 INTEGER :: ii, ij, jk, ikb 51 REAL(wp) :: zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn, zD, zvb, zub, ztb 52 REAL(wp) :: zMv, zMe, zMb, zmelt, zdvo, zdvob, zdva, zdM, zSs, zdMe, zdMb, zdMv 53 53 REAL(wp) :: zMnew, zMnew1, zMnew2, zheat_hcflux, zheat_latent, z1_12 54 54 REAL(wp) :: zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb 55 REAL(wp) :: zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2 55 REAL(wp) :: zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2, zdepw 56 REAL(wp), DIMENSION(jpk) :: ztoce, zuoce, zvoce, ze3t, ztmp 56 57 TYPE(iceberg), POINTER :: this, next 57 58 TYPE(point) , POINTER :: pt … … 83 84 pt => this%current_point 84 85 nknberg = this%number(1) 86 87 CALL icb_utl_interp( pt%xi, pt%yj, & ! position 88 & pssu=pt%ssu, pua=pt%ua, & ! oce/atm velocities 89 & pssv=pt%ssv, pva=pt%va, & ! oce/atm velocities 90 & psst=pt%sst, pcn=pt%cn ) 91 85 92 IF ( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) THEN 86 93 CALL icb_utl_interp( pt%xi, pt%yj, pe1=pt%e1, pe2=pt%e2, & 87 & puo=pt%uo, pui=pt%ui, pua=pt%ua, pssh_i=pt%ssh_x, & 88 & pvo=pt%vo, pvi=pt%vi, pva=pt%va, pssh_j=pt%ssh_y, & 89 & psst=pt%sst, pcn=pt%cn, phi=pt%hi, & 90 & plat=pt%lat, plon=pt%lon ) 91 ELSE 92 CALL icb_utl_interp( pt%xi, pt%yj, & ! position 93 & puo=pt%uo, pua=pt%ua, & ! oce/atm velocities 94 & pvo=pt%vo, pva=pt%va, & ! oce/atm velocities 95 & psst=pt%sst, pcn=pt%cn ) ! sst and ice concentration 96 ! preparation Nacho add t3d and uo, vo, to basal 94 & pui=pt%ui, pssh_i=pt%ssh_x, & 95 & pvi=pt%vi, pssh_j=pt%ssh_y, & 96 & phi=pt%hi, & 97 & plat=pt%lat, plon=pt%lon ) 97 98 END IF 98 99 ! … … 101 102 zM = pt%mass 102 103 zT = pt%thickness ! total thickness 103 ! D = (rn_rho_bergs/pp_rho_seawater)*zT! draught (keel depth)104 zD = rho_berg_1_oce * zT ! draught (keel depth) 104 105 ! F = zT - D ! freeboard 105 106 zW = pt%width … … 114 115 115 116 ! Environment 116 zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 ) 117 zdva = SQRT( (pt%ua -pt%uo)**2 + (pt%va -pt%vo)**2 ) 118 zSs = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva ! Sea state (eqn M.A9) 119 117 ! default sst, ssu and ssv 118 ! ln_M2016: use temp, u and v profile 119 IF ( ln_M2016 ) THEN 120 121 ! load t, u, v and e3 profile at icb position 122 CALL icb_utl_interp( pt%xi, pt%yj, ptoce=ztoce, puoce=zuoce, pvoce=zvoce, pe3t=ze3t ) 123 124 !compute bottom level 125 CALL icb_utl_getkb( pt%kb, ze3t, zD ) 126 127 ikb = MIN(pt%kb,mbkt(ii,ij)) 128 ztb = ztoce(ikb) ! basal temperature 129 zub = zuoce(ikb) 130 zvb = zvoce(ikb) 131 ELSE 132 ztb = pt%sst 133 zub = pt%ssu 134 zvb = pt%ssv 135 END IF 136 137 zdvob = SQRT( (pt%uvel-zub)**2 + (pt%vvel-zvb)**2 ) ! relative basal velocity 138 zdva = SQRT( (pt%ua -pt%ssu)**2 + (pt%va -pt%ssv)**2 ) ! relative wind 139 zSs = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva ! Sea state (eqn M.A9) 140 ! 120 141 ! Melt rates in m/s (i.e. division by rday) 121 zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2) , 0._wp ) * z1_rday ! Buoyant convection at sides (eqn M.A10) 122 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 ) 123 zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3))) , 0._wp ) * z1_rday ! Wave erosion (eqn M.A8 ) 142 IF ( ln_M2016 ) THEN 143 ! Buoyant convection at sides (eqn M.A10) but averaging along all the iceberg draft 144 ztmp(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday 145 zMv = 0.0; zdepw = 0.0 146 DO jk = 1,ikb-1 147 zMv = zMv + ze3t(jk)*ztmp(jk) 148 zdepw = zdepw + ze3t(jk) 149 END DO 150 zMv = (zMv + (zD - zdepw)*ztmp(ikb)) / zD 151 ELSE 152 zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday 153 END IF 154 155 ! Basal turbulent melting (eqn M.A7 ) but using the basal temperature 156 zMb = MAX( 0.58_wp*(zdvob**0.8_wp)*(ztb+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday 157 158 ! Wave erosion (eqn M.A8 ) 159 zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3))) , 0._wp ) * z1_rday 124 160 125 161 IF( ln_operator_splitting ) THEN ! Operator split update of volume/mass … … 209 245 210 246 ! Rolling 211 zDn = ( rn_rho_bergs / pp_rho_seawater )* zTn ! draught (keel depth)247 zDn = rho_berg_1_oce * zTn ! draught (keel depth) 212 248 IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN 213 249 zT = zTn -
NEMO/branches/2020/tickets_icb_1900/src/OCE/ICB/icbtrj.F90
r13062 r13357 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/tickets_icb_1900/src/OCE/ICB/icbutl.F90
r13265 r13357 19 19 !!---------------------------------------------------------------------- 20 20 USE par_oce ! ocean parameters 21 USE oce, ONLY: ts, uu, vv 21 22 USE dom_oce ! ocean domain 22 23 USE in_out_manager ! IO parameters … … 34 35 PRIVATE 35 36 37 INTERFACE icb_utl_bilin_h 38 MODULE PROCEDURE icb_utl_bilin_2d_h, icb_utl_bilin_3d_h 39 END INTERFACE 40 36 41 PUBLIC icb_utl_copy ! routine called in icbstp module 42 PUBLIC icb_utl_getkb ! routine called in icbdyn and icbthm modules 37 43 PUBLIC icb_utl_interp ! routine called in icbdyn, icbthm modules 38 44 PUBLIC icb_utl_bilin_h ! routine called in icbdyn module … … 56 62 CONTAINS 57 63 58 SUBROUTINE icb_utl_copy( )64 SUBROUTINE icb_utl_copy( Kmm ) 59 65 !!---------------------------------------------------------------------- 60 66 !! *** ROUTINE icb_utl_copy *** … … 64 70 !! ** Method : - blah blah 65 71 !!---------------------------------------------------------------------- 72 REAL(wp), DIMENSION(0:jpi+1,0:jpj+1) :: ztmp 66 73 #if defined key_si3 67 74 REAL(wp), DIMENSION(jpi,jpj) :: zssh_lead_m ! ocean surface (ssh_m) if ice is not embedded 68 75 ! ! ocean surface in leads if ice is embedded 69 76 #endif 77 INTEGER :: jk ! vertical loop index 78 INTEGER :: Kmm ! ocean time levelindex 79 ! 70 80 ! copy nemo forcing arrays into iceberg versions with extra halo 71 81 ! only necessary for variables not on T points 72 82 ! and ssh which is used to calculate gradients 73 74 uo_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 75 vo_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 76 tt_e(1:jpi,1:jpj) = sst_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 ! 84 ! surface forcing 85 ! 86 ssu_e(1:jpi,1:jpj) = ssu_m(:,:) * umask(:,:,1) 87 ssv_e(1:jpi,1:jpj) = ssv_m(:,:) * vmask(:,:,1) 88 sst_e(1:jpi,1:jpj) = sst_m(:,:) 89 fr_e (1:jpi,1:jpj) = fr_i (:,:) 90 ua_e (1:jpi,1:jpj) = utau (:,:) * umask(:,:,1) ! maybe mask useless because mask applied in sbcblk 91 va_e (1:jpi,1:jpj) = vtau (:,:) * vmask(:,:,1) ! maybe mask useless because mask applied in sbcblk 92 ! 93 CALL lbc_lnk_icb( 'icbutl', ssu_e, 'U', -1._wp, 1, 1 ) 94 CALL lbc_lnk_icb( 'icbutl', ssv_e, 'V', -1._wp, 1, 1 ) 83 95 CALL lbc_lnk_icb( 'icbutl', ua_e, 'U', -1._wp, 1, 1 ) 84 96 CALL lbc_lnk_icb( 'icbutl', va_e, 'V', -1._wp, 1, 1 ) … … 98 110 #endif 99 111 ! 112 ! (PM) could be improve with a 3d lbclnk gathering both variables 113 ! should be done once extra haloe generalised 114 IF ( ln_M2016 ) THEN 115 DO jk = 1,jpk 116 ! uoce 117 ztmp(1:jpi,1:jpj) = uu(:,:,jk,Kmm) 118 CALL lbc_lnk_icb( 'icbutl', ztmp, 'U', -1._wp, 1, 1 ) 119 uoce_e(:,:,jk) = ztmp(:,:) 120 ! 121 ! voce 122 ztmp(1:jpi,1:jpj) = vv(:,:,jk,Kmm) 123 CALL lbc_lnk_icb( 'icbutl', ztmp, 'V', -1._wp, 1, 1 ) 124 voce_e(:,:,jk) = ztmp(:,:) 125 END DO 126 toce_e(1:jpi,1:jpj,1:jpk) = ts(:,:,:,1,Kmm) 127 e3t_e (1:jpi,1:jpj,1:jpk) = e3t(:,:,:,Kmm) 128 END IF 129 ! 100 130 END SUBROUTINE icb_utl_copy 101 131 102 132 103 SUBROUTINE icb_utl_interp( pi, pj, pe1, puo, pui, pua, pssh_i, & 104 & pe2, pvo, pvi, pva, pssh_j, & 105 & psst, pcn, phi, pff, plon, plat) 133 SUBROUTINE icb_utl_interp( pi, pj, pe1, pssu, pui, pua, pssh_i, & 134 & pe2, pssv, pvi, pva, pssh_j, & 135 & psst, pcn, phi, pff, plon, plat, & 136 & ptoce, puoce, pvoce, pe3t ) 106 137 !!---------------------------------------------------------------------- 107 138 !! *** ROUTINE icb_utl_interp *** … … 122 153 REAL(wp), INTENT(in ) :: pi , pj ! position in (i,j) referential 123 154 REAL(wp), INTENT( out), OPTIONAL :: pe1, pe2 ! i- and j scale factors 124 REAL(wp), INTENT( out), OPTIONAL :: p uo, pvo, pui, pvi, pua, pva ! ocean, ice and wind speeds155 REAL(wp), INTENT( out), OPTIONAL :: pssu, pssv, pui, pvi, pua, pva ! ocean, ice and wind speeds 125 156 REAL(wp), INTENT( out), OPTIONAL :: pssh_i, pssh_j ! ssh i- & j-gradients 126 157 REAL(wp), INTENT( out), OPTIONAL :: psst, pcn, phi, pff ! SST, ice concentration, ice thickness, Coriolis 127 158 REAL(wp), INTENT( out), OPTIONAL :: plat, plon ! position 159 REAL(wp), DIMENSION(jpk), INTENT( out), OPTIONAL :: ptoce, puoce, pvoce, pe3t ! 3D variables 128 160 ! 129 161 REAL(wp), DIMENSION(4) :: zwT , zwU , zwV , zwF ! interpolation weight … … 147 179 IF ( PRESENT(plat) ) plat= icb_utl_bilin_h( rlat_e, iiT, ijT, zwT, .false. ) 148 180 ! 149 IF ( PRESENT(p uo ) ) puo = icb_utl_bilin_h( uo_e, iiU, ijU, zwU , .false. ) ! ocean velocities150 IF ( PRESENT(p vo ) ) pvo = icb_utl_bilin_h( vo_e, iiV, ijV, zwV , .false. ) !151 IF ( PRESENT(psst) ) psst = icb_utl_bilin_h( tt_e, iiT, ijT, zwT * zmskT, .false. ) ! sst152 IF ( PRESENT(pcn ) ) pcn = icb_utl_bilin_h( fr_e , iiT, ijT, zwT * zmskT, .false. ) ! ice concentration153 IF ( PRESENT(pff ) ) pff = icb_utl_bilin_h( ff_e , iiF, ijF, zwF , .false. ) ! Coriolis parameter181 IF ( PRESENT(pssu) ) pssu = icb_utl_bilin_h( ssu_e, iiU, ijU, zwU , .false. ) ! ocean velocities 182 IF ( PRESENT(pssv) ) pssv = icb_utl_bilin_h( ssv_e, iiV, ijV, zwV , .false. ) ! 183 IF ( PRESENT(psst) ) psst = icb_utl_bilin_h( sst_e, iiT, ijT, zwT * zmskT, .false. ) ! sst 184 IF ( PRESENT(pcn ) ) pcn = icb_utl_bilin_h( fr_e , iiT, ijT, zwT * zmskT, .false. ) ! ice concentration 185 IF ( PRESENT(pff ) ) pff = icb_utl_bilin_h( ff_e , iiF, ijF, zwF , .false. ) ! Coriolis parameter 154 186 ! 155 187 IF ( PRESENT(pua) .AND. PRESENT(pva) ) THEN … … 188 220 & icb_utl_bilin_h( ssh_e, iiTm, ijTm, zwTm*zmskTm, .false. ) ) / ( 0.2_wp * pe2 ) 189 221 END IF 222 ! 223 ! 3d interpolation 224 IF ( PRESENT(puoce) .AND. PRESENT(pvoce) ) THEN 225 puoce(:) = icb_utl_bilin_h( uoce_e , iiU, ijU, zwU ) 226 pvoce(:) = icb_utl_bilin_h( voce_e , iiV, ijV, zwV ) 227 END IF 228 IF ( PRESENT(ptoce) ) ptoce(:) = icb_utl_bilin_h( toce_e , iiT, ijT, zwT * zmskT ) 229 IF ( PRESENT(pe3t) ) pe3t(:) = e3t_e(iiT,ijT,:) ! as in Nacho tarball need to be fix once we are able to reproduce Nacho results 190 230 ! 191 231 END SUBROUTINE icb_utl_interp … … 291 331 END SUBROUTINE icb_utl_pos 292 332 293 REAL(wp) FUNCTION icb_utl_bilin_ h( pfld, pii, pij, pw, pllon )333 REAL(wp) FUNCTION icb_utl_bilin_2d_h( pfld, pii, pij, pw, pllon ) 294 334 !!---------------------------------------------------------------------- 295 335 !! *** FUNCTION icb_utl_bilin *** … … 321 361 ! 322 362 ! compute interpolated value 323 icb_utl_bilin_h = ( zdat(1)*pw(1) + zdat(2)*pw(2) + zdat(3)*pw(3) + zdat(4)*pw(4) ) / MAX(1.e-20, pw(1)+pw(2)+pw(3)+pw(4)) 324 ! 325 IF( pllon .AND. icb_utl_bilin_h > 180._wp ) icb_utl_bilin_h = icb_utl_bilin_h - 360._wp 326 ! 327 END FUNCTION icb_utl_bilin_h 363 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)) 364 ! 365 IF( pllon .AND. icb_utl_bilin_2d_h > 180._wp ) icb_utl_bilin_2d_h = icb_utl_bilin_2d_h - 360._wp 366 ! 367 END FUNCTION icb_utl_bilin_2d_h 368 369 FUNCTION icb_utl_bilin_3d_h( pfld, pii, pij, pw ) 370 !!---------------------------------------------------------------------- 371 !! *** FUNCTION icb_utl_bilin *** 372 !! 373 !! ** Purpose : bilinear interpolation at berg location depending on the grid-point type 374 !! this version deals with extra halo points 375 !! 376 !! !!gm CAUTION an optional argument should be added to handle 377 !! the slip/no-slip conditions ==>>> to be done later 378 !! 379 !!---------------------------------------------------------------------- 380 REAL(wp), DIMENSION(0:jpi+1,0:jpj+1, jpk), INTENT(in) :: pfld ! field to be interpolated 381 REAL(wp), DIMENSION(4) , INTENT(in) :: pw ! weight 382 INTEGER , INTENT(in) :: pii, pij ! bottom left corner 383 REAL(wp), DIMENSION(jpk) :: icb_utl_bilin_3d_h 384 ! 385 REAL(wp), DIMENSION(4,jpk) :: zdat ! input data 386 !!---------------------------------------------------------------------- 387 ! 388 ! data 389 zdat(1,:) = pfld(pii ,pij ,:) 390 zdat(2,:) = pfld(pii+1,pij ,:) 391 zdat(3,:) = pfld(pii ,pij+1,:) 392 zdat(4,:) = pfld(pii+1,pij+1,:) 393 ! 394 ! compute interpolated value 395 icb_utl_bilin_3d_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)) 396 ! 397 END FUNCTION icb_utl_bilin_3d_h 328 398 329 399 REAL(wp) FUNCTION icb_utl_bilin_e( pet, peu, pev, pef, pi, pj ) … … 405 475 END FUNCTION icb_utl_bilin_e 406 476 477 SUBROUTINE icb_utl_getkb( kb, pe3, pD ) 478 !!---------------------------------------------------------------------- 479 !! *** ROUTINE icb_utl_getkb *** 480 !! 481 !! ** Purpose : compute the latest level affected by icb 482 !! 483 !!---------------------------------------------------------------------- 484 INTEGER, INTENT(out) :: kb 485 REAL(wp), DIMENSION(:), INTENT(in) :: pe3 486 REAL(wp), INTENT(in) :: pD 487 !! 488 INTEGER :: jk 489 REAL(wp) :: zdepw 490 !! 491 zdepw = 0.0 492 kb = 1 493 DO WHILE ( zdepw < pD) 494 zdepw = zdepw + pe3(kb) 495 kb = kb + 1 496 END DO 497 kb = kb - 1 498 END SUBROUTINE 407 499 408 500 SUBROUTINE icb_utl_add( bergvals, ptvals ) … … 593 685 WRITE(numicb, 9200) kt, berg%number(1), & 594 686 pt%xi, pt%yj, pt%lon, pt%lat, pt%uvel, pt%vvel, & 595 pt% uo, pt%vo, pt%ua, pt%va, pt%ui, pt%vi687 pt%ssu, pt%ssv, pt%ua, pt%va, pt%ui, pt%vi 596 688 CALL flush( numicb ) 597 689 9200 FORMAT(5x,i5,2x,i10,6(2x,2f10.4)) … … 619 711 WRITE(numicb,'(a," pe=(",i3,")")' ) cd_label, narea 620 712 WRITE(numicb,'(a8,4x,a6,12x,a5,15x,a7,19x,a3,17x,a5,17x,a5,17x,a5)' ) & 621 & 'timestep', 'number', 'xi,yj','lon,lat','u,v',' uo,vo','ua,va','ui,vi'713 & 'timestep', 'number', 'xi,yj','lon,lat','u,v','ssu,ssv','ua,va','ui,vi' 622 714 ENDIF 623 715 DO WHILE( ASSOCIATED(this) ) -
NEMO/branches/2020/tickets_icb_1900/src/OCE/SBC/sbcmod.F90
r13226 r13357 456 456 457 457 IF( ln_icebergs ) THEN 458 CALL icb_stp( kt ) ! compute icebergs458 CALL icb_stp( kt, Kmm ) ! compute icebergs 459 459 ! Icebergs do not melt over the haloes. 460 460 ! So emp values over the haloes are no more consistent with the inner domain values.
Note: See TracChangeset
for help on using the changeset viewer.