| 1 | MODULE icbthm |
|---|
| 2 | !!====================================================================== |
|---|
| 3 | !! *** MODULE icbthm *** |
|---|
| 4 | !! Icebergs: thermodynamics routines for icebergs |
|---|
| 5 | !!====================================================================== |
|---|
| 6 | !! History : 3.3.1 ! 2010-01 (Martin&Adcroft) Original code |
|---|
| 7 | !! - ! 2011-03 (Madec) Part conversion to NEMO form |
|---|
| 8 | !! - ! Removal of mapping from another grid |
|---|
| 9 | !! - ! 2011-04 (Alderson) Split into separate modules |
|---|
| 10 | !! - ! 2011-05 (Alderson) Use tmask instead of tmask_i |
|---|
| 11 | !!---------------------------------------------------------------------- |
|---|
| 12 | !!---------------------------------------------------------------------- |
|---|
| 13 | !! icb_thm : initialise |
|---|
| 14 | !! reference for equations - M = Martin + Adcroft, OM 34, 2010 |
|---|
| 15 | !!---------------------------------------------------------------------- |
|---|
| 16 | USE par_oce ! NEMO parameters |
|---|
| 17 | USE dom_oce ! NEMO domain |
|---|
| 18 | USE in_out_manager ! NEMO IO routines, numout in particular |
|---|
| 19 | USE lib_mpp ! NEMO MPI routines, ctl_stop in particular |
|---|
| 20 | USE phycst ! NEMO physical constants |
|---|
| 21 | USE sbc_oce |
|---|
| 22 | USE eosbn2 ! equation of state |
|---|
| 23 | USE lib_fortran, ONLY : DDPDD |
|---|
| 24 | |
|---|
| 25 | USE icb_oce ! define iceberg arrays |
|---|
| 26 | USE icbutl ! iceberg utility routines |
|---|
| 27 | USE icbdia ! iceberg budget routines |
|---|
| 28 | |
|---|
| 29 | IMPLICIT NONE |
|---|
| 30 | PRIVATE |
|---|
| 31 | |
|---|
| 32 | PUBLIC icb_thm ! routine called in icbstp.F90 module |
|---|
| 33 | |
|---|
| 34 | !! * Substitutions |
|---|
| 35 | # include "single_precision_substitute.h90" |
|---|
| 36 | |
|---|
| 37 | !!---------------------------------------------------------------------- |
|---|
| 38 | !! NEMO/OCE 4.0 , NEMO Consortium (2018) |
|---|
| 39 | !! $Id$ |
|---|
| 40 | !! Software governed by the CeCILL license (see ./LICENSE) |
|---|
| 41 | !!---------------------------------------------------------------------- |
|---|
| 42 | CONTAINS |
|---|
| 43 | |
|---|
| 44 | SUBROUTINE icb_thm( kt ) |
|---|
| 45 | !!---------------------------------------------------------------------- |
|---|
| 46 | !! *** ROUTINE icb_thm *** |
|---|
| 47 | !! |
|---|
| 48 | !! ** Purpose : compute the iceberg thermodynamics. |
|---|
| 49 | !! |
|---|
| 50 | !! ** Method : - See Martin & Adcroft, Ocean Modelling 34, 2010 |
|---|
| 51 | !!---------------------------------------------------------------------- |
|---|
| 52 | INTEGER, INTENT(in) :: kt ! timestep number, just passed to icb_utl_print_berg |
|---|
| 53 | ! |
|---|
| 54 | INTEGER :: ii, ij, jk, ikb |
|---|
| 55 | REAL(wp) :: zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn, zD, zvb, zub, ztb |
|---|
| 56 | REAL(wp) :: zMv, zMe, zMb, zmelt, zdvo, zdvob, zdva, zdM, zSs, zdMe, zdMb, zdMv |
|---|
| 57 | REAL(wp) :: zSSS, zfzpt |
|---|
| 58 | REAL(wp) :: zMnew, zMnew1, zMnew2, zheat_hcflux, zheat_latent, z1_12 |
|---|
| 59 | REAL(wp) :: zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb |
|---|
| 60 | REAL(wp) :: zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2, zdepw |
|---|
| 61 | REAL(wp), DIMENSION(jpk) :: ztoce, zuoce, zvoce, ze3t, zzMv |
|---|
| 62 | TYPE(iceberg), POINTER :: this, next |
|---|
| 63 | TYPE(point) , POINTER :: pt |
|---|
| 64 | ! |
|---|
| 65 | COMPLEX(dp), DIMENSION(jpi,jpj) :: cicb_melt, cicb_hflx |
|---|
| 66 | !!---------------------------------------------------------------------- |
|---|
| 67 | ! |
|---|
| 68 | !! initialiaze cicb_melt and cicb_heat |
|---|
| 69 | cicb_melt = CMPLX( 0.e0, 0.e0, dp ) |
|---|
| 70 | cicb_hflx = CMPLX( 0.e0, 0.e0, dp ) |
|---|
| 71 | ! |
|---|
| 72 | z1_rday = 1._wp / rday |
|---|
| 73 | z1_12 = 1._wp / 12._wp |
|---|
| 74 | zdt = berg_dt |
|---|
| 75 | z1_dt = 1._wp / zdt |
|---|
| 76 | ! |
|---|
| 77 | ! we're either going to ignore berg fresh water melt flux and associated heat |
|---|
| 78 | ! or we pass it into the ocean, so at this point we set them both to zero, |
|---|
| 79 | ! accumulate the contributions to them from each iceberg in the while loop following |
|---|
| 80 | ! and then pass them (or not) to the ocean |
|---|
| 81 | ! |
|---|
| 82 | berg_grid%floating_melt(:,:) = 0._wp |
|---|
| 83 | ! calving_hflx re-used here as temporary workspace for the heat flux associated with melting |
|---|
| 84 | berg_grid%calving_hflx(:,:) = 0._wp |
|---|
| 85 | ! |
|---|
| 86 | this => first_berg |
|---|
| 87 | DO WHILE( ASSOCIATED(this) ) |
|---|
| 88 | ! |
|---|
| 89 | pt => this%current_point |
|---|
| 90 | nknberg = this%number(1) |
|---|
| 91 | |
|---|
| 92 | CALL icb_utl_interp( pt%xi, pt%yj, & ! position |
|---|
| 93 | & pssu=pt%ssu, pua=pt%ua, & ! oce/atm velocities |
|---|
| 94 | & pssv=pt%ssv, pva=pt%va, & ! oce/atm velocities |
|---|
| 95 | & psst=pt%sst, pcn=pt%cn, & |
|---|
| 96 | & psss=pt%sss ) |
|---|
| 97 | |
|---|
| 98 | IF ( nn_sample_rate > 0 .AND. MOD(kt-1,nn_sample_rate) == 0 ) THEN |
|---|
| 99 | CALL icb_utl_interp( pt%xi, pt%yj, pe1=pt%e1, pe2=pt%e2, & |
|---|
| 100 | & pui=pt%ui, pssh_i=pt%ssh_x, & |
|---|
| 101 | & pvi=pt%vi, pssh_j=pt%ssh_y, & |
|---|
| 102 | & phi=pt%hi, & |
|---|
| 103 | & plat=pt%lat, plon=pt%lon ) |
|---|
| 104 | END IF |
|---|
| 105 | ! |
|---|
| 106 | zSST = pt%sst |
|---|
| 107 | zSSS = pt%sss |
|---|
| 108 | CALL eos_fzp(zSSS,zfzpt) ! freezing point |
|---|
| 109 | zIC = MIN( 1._wp, pt%cn + rn_sicn_shift ) ! Shift sea-ice concentration !!gm ??? |
|---|
| 110 | zM = pt%mass |
|---|
| 111 | zT = pt%thickness ! total thickness |
|---|
| 112 | zD = rho_berg_1_oce * zT ! draught (keel depth) |
|---|
| 113 | zW = pt%width |
|---|
| 114 | zL = pt%length |
|---|
| 115 | zxi = pt%xi ! position in (i,j) referential |
|---|
| 116 | zyj = pt%yj |
|---|
| 117 | ii = INT( zxi + 0.5 ) ! T-cell of the berg |
|---|
| 118 | ii = mi1( ii ) |
|---|
| 119 | ij = INT( zyj + 0.5 ) |
|---|
| 120 | ij = mj1( ij ) |
|---|
| 121 | zVol = zT * zW * zL |
|---|
| 122 | |
|---|
| 123 | ! Environment |
|---|
| 124 | ! default sst, ssu and ssv |
|---|
| 125 | ! ln_M2016: use temp, u and v profile |
|---|
| 126 | IF ( ln_M2016 ) THEN |
|---|
| 127 | |
|---|
| 128 | ! load t, u, v and e3 profile at icb position |
|---|
| 129 | CALL icb_utl_interp( pt%xi, pt%yj, ptoce=ztoce, puoce=zuoce, pvoce=zvoce, pe3t=ze3t ) |
|---|
| 130 | |
|---|
| 131 | !compute bottom level |
|---|
| 132 | CALL icb_utl_getkb( pt%kb, ze3t, zD ) |
|---|
| 133 | |
|---|
| 134 | ikb = MIN(pt%kb,mbkt(ii,ij)) ! limit pt%kb by mbkt |
|---|
| 135 | ! => bottom temperature used to fill ztoce(mbkt:jpk) |
|---|
| 136 | ztb = ztoce(ikb) ! basal temperature |
|---|
| 137 | zub = zuoce(ikb) |
|---|
| 138 | zvb = zvoce(ikb) |
|---|
| 139 | ELSE |
|---|
| 140 | ztb = pt%sst |
|---|
| 141 | zub = pt%ssu |
|---|
| 142 | zvb = pt%ssv |
|---|
| 143 | END IF |
|---|
| 144 | |
|---|
| 145 | zdvob = SQRT( (pt%uvel-zub)**2 + (pt%vvel-zvb)**2 ) ! relative basal velocity |
|---|
| 146 | zdva = SQRT( (pt%ua -pt%ssu)**2 + (pt%va -pt%ssv)**2 ) ! relative wind |
|---|
| 147 | zSs = 1.5_wp * SQRT( zdva ) + 0.1_wp * zdva ! Sea state (eqn M.A9) |
|---|
| 148 | ! |
|---|
| 149 | ! Melt rates in m/s (i.e. division by rday) |
|---|
| 150 | ! Buoyant convection at sides (eqn M.A10) |
|---|
| 151 | IF ( ln_M2016 ) THEN |
|---|
| 152 | ! averaging along all the iceberg draft |
|---|
| 153 | zzMv(:) = MAX( 7.62d-3*ztoce(:)+1.29d-3*(ztoce(:)**2), 0._wp ) * z1_rday |
|---|
| 154 | CALL icb_utl_zavg(zMv, zzMv, ze3t, zD, ikb ) |
|---|
| 155 | ELSE |
|---|
| 156 | zMv = MAX( 7.62d-3*zSST+1.29d-3*(zSST**2), 0._wp ) * z1_rday |
|---|
| 157 | END IF |
|---|
| 158 | ! |
|---|
| 159 | ! Basal turbulent melting (eqn M.A7 ) |
|---|
| 160 | IF ( zSST > zfzpt ) THEN ! Calculate basal melting only if SST above freezing point |
|---|
| 161 | zMb = MAX( 0.58_wp*(zdvob**0.8_wp)*(ztb+4.0_wp)/(zL**0.2_wp) , 0._wp ) * z1_rday |
|---|
| 162 | ELSE |
|---|
| 163 | zMb = 0._wp ! No basal melting if SST below freezing point |
|---|
| 164 | ENDIF |
|---|
| 165 | ! |
|---|
| 166 | ! Wave erosion (eqn M.A8 ) |
|---|
| 167 | zMe = MAX( z1_12*(zSST+2.)*zSs*(1._wp+COS(rpi*(zIC**3))) , 0._wp ) * z1_rday |
|---|
| 168 | |
|---|
| 169 | IF( ln_operator_splitting ) THEN ! Operator split update of volume/mass |
|---|
| 170 | zTn = MAX( zT - zMb*zdt , 0._wp ) ! new total thickness (m) |
|---|
| 171 | znVol = zTn * zW * zL ! new volume (m^3) |
|---|
| 172 | zMnew1 = ( znVol / zVol ) * zM ! new mass (kg) |
|---|
| 173 | zdMb = zM - zMnew1 ! mass lost to basal melting (>0) (kg) |
|---|
| 174 | ! |
|---|
| 175 | zLn = MAX( zL - zMv*zdt , 0._wp ) ! new length (m) |
|---|
| 176 | zWn = MAX( zW - zMv*zdt , 0._wp ) ! new width (m) |
|---|
| 177 | znVol = zTn * zWn * zLn ! new volume (m^3) |
|---|
| 178 | zMnew2 = ( znVol / zVol ) * zM ! new mass (kg) |
|---|
| 179 | zdMv = zMnew1 - zMnew2 ! mass lost to buoyant convection (>0) (kg) |
|---|
| 180 | ! |
|---|
| 181 | zLn = MAX( zLn - zMe*zdt , 0._wp ) ! new length (m) |
|---|
| 182 | zWn = MAX( zWn - zMe*zdt , 0._wp ) ! new width (m) |
|---|
| 183 | znVol = zTn * zWn * zLn ! new volume (m^3) |
|---|
| 184 | zMnew = ( znVol / zVol ) * zM ! new mass (kg) |
|---|
| 185 | zdMe = zMnew2 - zMnew ! mass lost to erosion (>0) (kg) |
|---|
| 186 | zdM = zM - zMnew ! mass lost to all erosion and melting (>0) (kg) |
|---|
| 187 | ! |
|---|
| 188 | ELSE ! Update dimensions of berg |
|---|
| 189 | zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp ) ! (m) |
|---|
| 190 | zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp ) ! (m) |
|---|
| 191 | zTn = MAX( zT - zMb *zdt ,0._wp ) ! (m) |
|---|
| 192 | ! Update volume and mass of berg |
|---|
| 193 | znVol = zTn*zWn*zLn ! (m^3) |
|---|
| 194 | zMnew = (znVol/zVol)*zM ! (kg) |
|---|
| 195 | zdM = zM - zMnew ! (kg) |
|---|
| 196 | zdMb = (zM/zVol) * (zW* zL ) *zMb*zdt ! approx. mass loss to basal melting (kg) |
|---|
| 197 | zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt ! approx. mass lost to erosion (kg) |
|---|
| 198 | zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt ! approx. mass loss to buoyant convection (kg) |
|---|
| 199 | ENDIF |
|---|
| 200 | |
|---|
| 201 | IF( rn_bits_erosion_fraction > 0._wp ) THEN ! Bergy bits |
|---|
| 202 | ! |
|---|
| 203 | zMbits = pt%mass_of_bits ! mass of bergy bits (kg) |
|---|
| 204 | zdMbitsE = rn_bits_erosion_fraction * zdMe ! change in mass of bits (kg) |
|---|
| 205 | znMbits = zMbits + zdMbitsE ! add new bergy bits to mass (kg) |
|---|
| 206 | zLbits = MIN( zL, zW, zT, 40._wp ) ! assume bergy bits are smallest dimension or 40 meters |
|---|
| 207 | zAbits = ( zMbits / rn_rho_bergs ) / zLbits ! Effective bottom area (assuming T=Lbits) |
|---|
| 208 | zMbb = MAX( 0.58_wp*(zdvo**0.8_wp)*(zSST+2._wp) / & |
|---|
| 209 | & ( zLbits**0.2_wp ) , 0._wp ) * z1_rday ! Basal turbulent melting (for bits) |
|---|
| 210 | zMbb = rn_rho_bergs * zAbits * zMbb ! in kg/s |
|---|
| 211 | zdMbitsM = MIN( zMbb*zdt , znMbits ) ! bergy bits mass lost to melting (kg) |
|---|
| 212 | znMbits = znMbits-zdMbitsM ! remove mass lost to bergy bits melt |
|---|
| 213 | IF( zMnew == 0._wp ) THEN ! if parent berg has completely melted then |
|---|
| 214 | zdMbitsM = zdMbitsM + znMbits ! instantly melt all the bergy bits |
|---|
| 215 | znMbits = 0._wp |
|---|
| 216 | ENDIF |
|---|
| 217 | ELSE ! No bergy bits |
|---|
| 218 | zAbits = 0._wp |
|---|
| 219 | zdMbitsE = 0._wp |
|---|
| 220 | zdMbitsM = 0._wp |
|---|
| 221 | znMbits = pt%mass_of_bits ! retain previous value incase non-zero |
|---|
| 222 | ENDIF |
|---|
| 223 | |
|---|
| 224 | ! use tmask rather than tmask_i when dealing with icebergs |
|---|
| 225 | IF( tmask(ii,ij,1) /= 0._wp ) THEN ! Add melting to the grid and field diagnostics |
|---|
| 226 | z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling |
|---|
| 227 | z1_dt_e1e2 = z1_dt * z1_e1e2 |
|---|
| 228 | ! |
|---|
| 229 | ! iceberg melt |
|---|
| 230 | !! the use of DDPDD function for the cumulative sum is needed for reproducibility |
|---|
| 231 | zmelt = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt ! kg/s |
|---|
| 232 | CALL DDPDD( CMPLX( zmelt * z1_e1e2, 0.e0, dp ), cicb_melt(ii,ij) ) |
|---|
| 233 | ! |
|---|
| 234 | ! iceberg heat flux |
|---|
| 235 | !! the use of DDPDD function for the cumulative sum is needed for reproducibility |
|---|
| 236 | !! NB. The src_calving_hflx field is currently hardwired to zero in icb_stp, which means that the |
|---|
| 237 | !! heat density of the icebergs is zero and the heat content flux to the ocean from iceberg |
|---|
| 238 | !! melting is always zero. Leaving the term in the code until such a time as this is fixed. DS. |
|---|
| 239 | zheat_hcflux = zmelt * pt%heat_density ! heat content flux : kg/s x J/kg = J/s |
|---|
| 240 | zheat_latent = - zmelt * rLfus ! latent heat flux: kg/s x J/kg = J/s |
|---|
| 241 | CALL DDPDD( CMPLX( ( zheat_hcflux + zheat_latent ) * z1_e1e2, 0.e0, dp ), cicb_hflx(ii,ij) ) |
|---|
| 242 | ! |
|---|
| 243 | ! diagnostics |
|---|
| 244 | CALL icb_dia_melt( ii, ij, zMnew, zheat_hcflux, zheat_latent, this%mass_scaling, & |
|---|
| 245 | & zdM, zdMbitsE, zdMbitsM, zdMb, zdMe, & |
|---|
| 246 | & zdMv, z1_dt_e1e2 ) |
|---|
| 247 | ELSE |
|---|
| 248 | WRITE(numout,*) 'icb_thm: berg ',this%number(:),' appears to have grounded at ',narea,ii,ij |
|---|
| 249 | CALL icb_utl_print_berg( this, kt ) |
|---|
| 250 | WRITE(numout,*) 'msk=',tmask(ii,ij,1), e1e2t(ii,ij) |
|---|
| 251 | CALL ctl_stop('icb_thm', 'berg appears to have grounded!') |
|---|
| 252 | ENDIF |
|---|
| 253 | |
|---|
| 254 | ! Rolling |
|---|
| 255 | zDn = rho_berg_1_oce * zTn ! draught (keel depth) |
|---|
| 256 | IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN |
|---|
| 257 | zT = zTn |
|---|
| 258 | zTn = zWn |
|---|
| 259 | zWn = zT |
|---|
| 260 | ENDIF |
|---|
| 261 | |
|---|
| 262 | ! Store the new state of iceberg (with L>W) |
|---|
| 263 | pt%mass = zMnew |
|---|
| 264 | pt%mass_of_bits = znMbits |
|---|
| 265 | pt%thickness = zTn |
|---|
| 266 | pt%width = MIN( zWn , zLn ) |
|---|
| 267 | pt%length = MAX( zWn , zLn ) |
|---|
| 268 | |
|---|
| 269 | next=>this%next |
|---|
| 270 | |
|---|
| 271 | !!gm add a test to avoid over melting ? |
|---|
| 272 | !!pm I agree, over melting could break conservation (more melt than calving) |
|---|
| 273 | |
|---|
| 274 | IF( zMnew <= 0._wp ) THEN ! Delete the berg if completely melted |
|---|
| 275 | CALL icb_utl_delete( first_berg, this ) |
|---|
| 276 | ! |
|---|
| 277 | ELSE ! Diagnose mass distribution on grid |
|---|
| 278 | z1_e1e2 = r1_e1e2t(ii,ij) * this%mass_scaling |
|---|
| 279 | CALL icb_dia_size( ii, ij, zWn, zLn, zAbits, & |
|---|
| 280 | & this%mass_scaling, zMnew, znMbits, z1_e1e2 ) |
|---|
| 281 | ENDIF |
|---|
| 282 | ! |
|---|
| 283 | this=>next |
|---|
| 284 | ! |
|---|
| 285 | END DO |
|---|
| 286 | ! |
|---|
| 287 | berg_grid%floating_melt = CASTDP(cicb_melt) ! kg/m2/s |
|---|
| 288 | berg_grid%calving_hflx = CASTDP(cicb_hflx) |
|---|
| 289 | ! |
|---|
| 290 | ! now use melt and associated heat flux in ocean (or not) |
|---|
| 291 | ! |
|---|
| 292 | IF(.NOT. ln_passive_mode ) THEN |
|---|
| 293 | emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:) |
|---|
| 294 | qns (:,:) = qns (:,:) + berg_grid%calving_hflx (:,:) |
|---|
| 295 | ENDIF |
|---|
| 296 | ! |
|---|
| 297 | END SUBROUTINE icb_thm |
|---|
| 298 | |
|---|
| 299 | !!====================================================================== |
|---|
| 300 | END MODULE icbthm |
|---|
