| 1 | MODULE icbthm |
|---|
| 2 | |
|---|
| 3 | !!====================================================================== |
|---|
| 4 | !! *** MODULE icbthm *** |
|---|
| 5 | !! Ocean physics: thermodynamics routines for icebergs |
|---|
| 6 | !!====================================================================== |
|---|
| 7 | !! History : 3.3.1 ! 2010-01 (Martin&Adcroft) Original code |
|---|
| 8 | !! - ! 2011-03 (Madec) Part conversion to NEMO form |
|---|
| 9 | !! - ! Removal of mapping from another grid |
|---|
| 10 | !! - ! 2011-04 (Alderson) Split into separate modules |
|---|
| 11 | !! - ! 2011-05 (Alderson) Use tmask instead of tmask_i |
|---|
| 12 | !!---------------------------------------------------------------------- |
|---|
| 13 | !!---------------------------------------------------------------------- |
|---|
| 14 | !! thermodynamics : initialise |
|---|
| 15 | !! reference for equations - M = Martin + Adcroft, OM 34, 2010 |
|---|
| 16 | !!---------------------------------------------------------------------- |
|---|
| 17 | USE par_oce ! NEMO parameters |
|---|
| 18 | USE dom_oce ! NEMO domain |
|---|
| 19 | USE in_out_manager ! NEMO IO routines, numout in particular |
|---|
| 20 | USE lib_mpp ! NEMO MPI routines, ctl_stop in particular |
|---|
| 21 | USE phycst ! NEMO physical constants |
|---|
| 22 | USE sbc_oce |
|---|
| 23 | |
|---|
| 24 | USE icb_oce ! define iceberg arrays |
|---|
| 25 | USE icbutl ! iceberg utility routines |
|---|
| 26 | USE icbdia ! iceberg budget routines |
|---|
| 27 | |
|---|
| 28 | IMPLICIT NONE |
|---|
| 29 | PRIVATE |
|---|
| 30 | |
|---|
| 31 | PUBLIC thermodynamics ! routine called in xxx.F90 module |
|---|
| 32 | |
|---|
| 33 | CONTAINS |
|---|
| 34 | |
|---|
| 35 | SUBROUTINE thermodynamics( kt ) |
|---|
| 36 | !!---------------------------------------------------------------------- |
|---|
| 37 | !! *** ROUTINE thermodynamics *** |
|---|
| 38 | !! |
|---|
| 39 | !! ** Purpose : compute the iceberg thermodynamics. |
|---|
| 40 | !! |
|---|
| 41 | !! ** Method : - blah blah |
|---|
| 42 | !!---------------------------------------------------------------------- |
|---|
| 43 | INTEGER :: kt ! timestep number, just passed to print_berg |
|---|
| 44 | ! |
|---|
| 45 | REAL(wp) :: M, T, W, L, SST, Vol, Ln, Wn, Tn, nVol, IC, Dn |
|---|
| 46 | REAL(wp) :: Mv, Me, Mb, melt, dvo, dva, dM, Ss, dMe, dMb, dMv |
|---|
| 47 | REAL(wp) :: Mnew, Mnew1, Mnew2, heat |
|---|
| 48 | REAL(wp) :: Mbits, nMbits, dMbitsE, dMbitsM, Lbits, Abits, Mbb |
|---|
| 49 | REAL(wp) :: xi, yj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2 |
|---|
| 50 | INTEGER :: ii, ij |
|---|
| 51 | TYPE(iceberg) , POINTER :: this, next |
|---|
| 52 | TYPE(point) , POINTER :: pt |
|---|
| 53 | !!---------------------------------------------------------------------- |
|---|
| 54 | ! |
|---|
| 55 | z1_rday = 1._wp / rday |
|---|
| 56 | |
|---|
| 57 | ! we're either going to ignore berg fresh water melt flux and associated heat |
|---|
| 58 | ! or we pass it into the ocean, so at this point we set them both to zero, |
|---|
| 59 | ! accumulate the contributions to them from each iceberg in the while loop following |
|---|
| 60 | ! and then pass them (or not) to the ocean |
|---|
| 61 | ! |
|---|
| 62 | berg_grid%floating_melt(:,:) = 0._wp |
|---|
| 63 | berg_grid%calving_hflx(:,:) = 0._wp |
|---|
| 64 | |
|---|
| 65 | this => first_berg |
|---|
| 66 | DO WHILE( associated(this) ) |
|---|
| 67 | ! |
|---|
| 68 | pt => this%current_point |
|---|
| 69 | knberg = this%number(1) |
|---|
| 70 | CALL interp_flds( pt%xi, pt%e1, pt%uo, pt%ui, pt%ua, pt%ssh_x, & |
|---|
| 71 | & pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y, pt%sst, pt%cn, pt%hi, zff ) |
|---|
| 72 | ! |
|---|
| 73 | SST = pt%sst |
|---|
| 74 | IC = MIN( 1._wp, pt%cn + rn_sicn_shift ) ! Shift sea-ice concentration !!gm ??? |
|---|
| 75 | M = pt%mass |
|---|
| 76 | T = pt%thickness ! total thickness |
|---|
| 77 | ! D = (rn_rho_bergs/rho_seawater)*T ! draught (keel depth) |
|---|
| 78 | ! F = T - D ! freeboard |
|---|
| 79 | W = pt%width |
|---|
| 80 | L = pt%length |
|---|
| 81 | xi = pt%xi ! position in (i,j) referential |
|---|
| 82 | yj = pt%yj |
|---|
| 83 | ii = INT( xi + 0.5 ) - nimpp + 1 ! t-cell of the berg |
|---|
| 84 | ij = INT( yj + 0.5 ) - njmpp + 1 |
|---|
| 85 | Vol = T * W * L |
|---|
| 86 | zdt = berg_dt ; z1_dt = 1._wp / zdt |
|---|
| 87 | |
|---|
| 88 | ! Environment |
|---|
| 89 | dvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 ) |
|---|
| 90 | dva = SQRT( (pt%ua -pt%uo)**2 + (pt%va -pt%vo)**2 ) |
|---|
| 91 | Ss = 1.5 * SQRT( dva ) + 0.1 * dva ! Sea state (eqn M.A9) |
|---|
| 92 | |
|---|
| 93 | ! Melt rates in m/s (i.e. division by rday) |
|---|
| 94 | Mv = MAX( 7.62e-3*SST+1.29e-3*(SST**2) , 0._wp ) * z1_rday ! Buoyant convection at sides (eqn M.A10) |
|---|
| 95 | Mb = MAX( 0.58*(dvo**0.8)*(SST+4.0)/(L**0.2) , 0._wp ) * z1_rday ! Basal turbulent melting (eqn M.A7 ) |
|---|
| 96 | Me = MAX( 1./12.*(SST+2.)*Ss*(1+cos(rpi*(IC**3))) , 0._wp ) * z1_rday ! Wave erosion (eqn M.A8 ) |
|---|
| 97 | |
|---|
| 98 | IF( ln_operator_splitting ) THEN ! Operator split update of volume/mass |
|---|
| 99 | Tn = MAX( T - Mb*zdt , 0._wp ) ! new total thickness (m) |
|---|
| 100 | nVol = Tn * W * L ! new volume (m^3) |
|---|
| 101 | Mnew1 = (nVol/Vol) * M ! new mass (kg) |
|---|
| 102 | dMb = M - Mnew1 ! mass lost to basal melting (>0) (kg) |
|---|
| 103 | ! |
|---|
| 104 | Ln = MAX( L - Mv*zdt , 0._wp ) ! new length (m) |
|---|
| 105 | Wn = MAX( W - Mv*zdt , 0._wp ) ! new width (m) |
|---|
| 106 | nVol = Tn * Wn * Ln ! new volume (m^3) |
|---|
| 107 | Mnew2 = (nVol/Vol) * M ! new mass (kg) |
|---|
| 108 | dMv = Mnew1 - Mnew2 ! mass lost to buoyant convection (>0) (kg) |
|---|
| 109 | ! |
|---|
| 110 | Ln = MAX( Ln - Me*zdt , 0._wp ) ! new length (m) |
|---|
| 111 | Wn = MAX( Wn - Me*zdt , 0._wp ) ! new width (m) |
|---|
| 112 | nVol = Tn * Wn * Ln ! new volume (m^3) |
|---|
| 113 | Mnew = ( nVol / Vol ) * M ! new mass (kg) |
|---|
| 114 | dMe = Mnew2 - Mnew ! mass lost to erosion (>0) (kg) |
|---|
| 115 | dM = M - Mnew ! mass lost to all erosion and melting (>0) (kg) |
|---|
| 116 | ! |
|---|
| 117 | ELSE ! Update dimensions of berg |
|---|
| 118 | Ln = MAX( L -(Mv+Me)*zdt ,0._wp ) ! (m) |
|---|
| 119 | Wn = MAX( W -(Mv+Me)*zdt ,0._wp ) ! (m) |
|---|
| 120 | Tn = MAX( T - Mb *zdt ,0._wp ) ! (m) |
|---|
| 121 | ! Update volume and mass of berg |
|---|
| 122 | nVol = Tn*Wn*Ln ! (m^3) |
|---|
| 123 | Mnew = (nVol/Vol)*M ! (kg) |
|---|
| 124 | dM = M - Mnew ! (kg) |
|---|
| 125 | dMb = (M/Vol) * (W* L ) *Mb*zdt ! approx. mass loss to basal melting (kg) |
|---|
| 126 | dMe = (M/Vol) * (T*(W+L)) *Me*zdt ! approx. mass lost to erosion (kg) |
|---|
| 127 | dMv = (M/Vol) * (T*(W+L)) *Mv*zdt ! approx. mass loss to buoyant convection (kg) |
|---|
| 128 | ENDIF |
|---|
| 129 | |
|---|
| 130 | IF( rn_bits_erosion_fraction > 0._wp ) THEN ! Bergy bits |
|---|
| 131 | ! |
|---|
| 132 | Mbits = pt%mass_of_bits ! mass of bergy bits (kg) |
|---|
| 133 | dMbitsE = rn_bits_erosion_fraction * dMe ! change in mass of bits (kg) |
|---|
| 134 | nMbits = Mbits + dMbitsE ! add new bergy bits to mass (kg) |
|---|
| 135 | Lbits = MIN( L, W, T, 40._wp ) ! assume bergy bits are smallest dimension or 40 meters |
|---|
| 136 | Abits = ( Mbits / rn_rho_bergs ) / Lbits ! Effective bottom area (assuming T=Lbits) |
|---|
| 137 | Mbb = MAX( 0.58*(dvo**0.8)*(SST+2.0)/(Lbits**0.2), 0.) * z1_rday ! Basal turbulent melting (for bits) |
|---|
| 138 | Mbb = rn_rho_bergs * Abits * Mbb ! in kg/s |
|---|
| 139 | dMbitsM = MIN( Mbb*zdt , nMbits ) ! bergy bits mass lost to melting (kg) |
|---|
| 140 | nMbits = nMbits-dMbitsM ! remove mass lost to bergy bits melt |
|---|
| 141 | IF( Mnew == 0._wp ) THEN ! if parent berg has completely melted then |
|---|
| 142 | dMbitsM = dMbitsM + nMbits ! instantly melt all the bergy bits |
|---|
| 143 | nMbits = 0._wp |
|---|
| 144 | ENDIF |
|---|
| 145 | ELSE ! No bergy bits |
|---|
| 146 | Abits = 0._wp |
|---|
| 147 | dMbitsE = 0._wp |
|---|
| 148 | dMbitsM = 0._wp |
|---|
| 149 | nMbits = pt%mass_of_bits ! retain previous value incase non-zero |
|---|
| 150 | ENDIF |
|---|
| 151 | |
|---|
| 152 | ! use tmask rather than tmask_i when dealing with icebergs |
|---|
| 153 | IF( tmask(ii,ij,1) /= 0._wp ) THEN ! Add melting to the grid and field diagnostics |
|---|
| 154 | z1_e1e2 = 1._wp / e1e2t(ii,ij) * this%mass_scaling |
|---|
| 155 | z1_dt_e1e2 = z1_dt * z1_e1e2 |
|---|
| 156 | melt = ( dM - ( dMbitsE - dMbitsM ) ) * z1_dt ! kg/s |
|---|
| 157 | berg_grid%floating_melt(ii,ij) = berg_grid%floating_melt(ii,ij) + melt * z1_e1e2 ! kg/m2/s |
|---|
| 158 | heat = melt * pt%heat_density ! kg/s x J/kg = J/s |
|---|
| 159 | berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + heat * z1_e1e2 ! W/m2 |
|---|
| 160 | CALL melt_budget(ii, ij, Mnew, heat, this%mass_scaling, dM, dMbitsE, dMbitsM, dMb, dMe, dMv, z1_dt_e1e2 ) |
|---|
| 161 | ELSE |
|---|
| 162 | WRITE(numout,*) 'thermodynamics: berg ',this%number(:),' appears to have grounded at ',narea,ii,ij |
|---|
| 163 | CALL print_berg( this, kt ) |
|---|
| 164 | WRITE(numout,*) 'msk=',tmask(ii,ij,1), e1e2t(ii,ij) |
|---|
| 165 | CALL ctl_stop('thermodynamics', 'berg appears to have grounded!') |
|---|
| 166 | ENDIF |
|---|
| 167 | |
|---|
| 168 | ! Rolling |
|---|
| 169 | Dn = ( rn_rho_bergs / rho_seawater ) * Tn ! draught (keel depth) |
|---|
| 170 | IF( Dn > 0._wp .AND. MAX(Wn,Ln) < SQRT( 0.92*(Dn**2) + 58.32*Dn ) ) THEN |
|---|
| 171 | T = Tn |
|---|
| 172 | Tn = Wn |
|---|
| 173 | Wn = T |
|---|
| 174 | endif |
|---|
| 175 | |
|---|
| 176 | ! Store the new state of iceberg (with L>W) |
|---|
| 177 | pt%mass = Mnew |
|---|
| 178 | pt%mass_of_bits = nMbits |
|---|
| 179 | pt%thickness = Tn |
|---|
| 180 | pt%width = min(Wn,Ln) |
|---|
| 181 | pt%length = max(Wn,Ln) |
|---|
| 182 | |
|---|
| 183 | next=>this%next |
|---|
| 184 | |
|---|
| 185 | !!gm add a test to avoid over melting ? |
|---|
| 186 | |
|---|
| 187 | IF( Mnew <= 0._wp ) THEN ! Delete the berg if completely melted |
|---|
| 188 | CALL delete_iceberg_from_list( first_berg, this ) |
|---|
| 189 | ! |
|---|
| 190 | ELSE ! Diagnose mass distribution on grid |
|---|
| 191 | z1_e1e2 = 1._wp / e1e2t(ii,ij) * this%mass_scaling |
|---|
| 192 | CALL size_budget(ii, ij, Wn, Ln, Abits, this%mass_scaling, Mnew, nMbits, z1_e1e2) |
|---|
| 193 | ENDIF |
|---|
| 194 | ! |
|---|
| 195 | this=>next |
|---|
| 196 | ! |
|---|
| 197 | END DO |
|---|
| 198 | |
|---|
| 199 | ! now use melt and associated heat flux in ocean (or not) |
|---|
| 200 | ! |
|---|
| 201 | IF(.NOT. ln_passive_mode ) THEN |
|---|
| 202 | emp (:,:) = emp (:,:) - berg_grid%floating_melt(:,:) |
|---|
| 203 | emps(:,:) = emps(:,:) - berg_grid%floating_melt(:,:) |
|---|
| 204 | !! qns (:,:) = qns (:,:) + berg_grid%calving_hflx(:,:) !! heat flux not yet properly coded |
|---|
| 205 | ENDIF |
|---|
| 206 | ! |
|---|
| 207 | END SUBROUTINE thermodynamics |
|---|
| 208 | |
|---|
| 209 | END MODULE icbthm |
|---|
