- Timestamp:
- 2012-04-18T12:42:56+02:00 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/ICB/icbthm.F90
r3339 r3359 29 29 PRIVATE 30 30 31 PUBLIC thermodynamics ! routine called in xxx.F90 module31 PUBLIC thermodynamics ! routine called in icbrun.F90 module 32 32 33 33 CONTAINS … … 43 43 INTEGER :: kt ! timestep number, just passed to print_berg 44 44 ! 45 REAL(wp) :: M, T, W, L, SST, Vol, Ln, Wn, Tn, nVol, IC,Dn46 REAL(wp) :: Mv, Me, Mb, melt, dvo, dva, dM, Ss, dMe, dMb,dMv47 REAL(wp) :: Mnew, Mnew1, Mnew2,heat48 REAL(wp) :: Mbits, nMbits, dMbitsE, dMbitsM, Lbits, Abits,Mbb49 REAL(wp) :: xi,yj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e245 REAL(wp) :: zM, zT, zW, zL, zSST, zVol, zLn, zWn, zTn, znVol, zIC, zDn 46 REAL(wp) :: zMv, zMe, zMb, zmelt, zdvo, zdva, zdM, zSs, zdMe, zdMb, zdMv 47 REAL(wp) :: zMnew, zMnew1, zMnew2, zheat 48 REAL(wp) :: zMbits, znMbits, zdMbitsE, zdMbitsM, zLbits, zAbits, zMbb 49 REAL(wp) :: zxi, zyj, zff, z1_rday, z1_e1e2, zdt, z1_dt, z1_dt_e1e2 50 50 INTEGER :: ii, ij 51 51 TYPE(iceberg) , POINTER :: this, next … … 67 67 ! 68 68 pt => this%current_point 69 knberg = this%number(1)69 nknberg = this%number(1) 70 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 71 & pt%yj, pt%e2, pt%vo, pt%vi, pt%va, pt%ssh_y, & 72 & pt%sst, pt%cn, pt%hi, zff ) 73 ! 74 zSST = pt%sst 75 zIC = MIN( 1._wp, pt%cn + rn_sicn_shift ) ! Shift sea-ice concentration !!gm ??? 76 zM = pt%mass 77 zT = pt%thickness ! total thickness 78 ! D = (rn_rho_bergs/pp_rho_seawater)*zT ! draught (keel depth) 79 ! F = zT - D ! freeboard 80 zW = pt%width 81 zL = pt%length 82 zxi = pt%xi ! position in (i,j) referential 83 zyj = pt%yj 84 ii = INT( zxi + 0.5 ) - nimpp + 1 ! t-cell of the berg 85 ij = INT( zyj + 0.5 ) - njmpp + 1 86 zVol = zT * zW * zL 86 87 zdt = berg_dt ; z1_dt = 1._wp / zdt 87 88 88 89 ! 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)90 zdvo = SQRT( (pt%uvel-pt%uo)**2 + (pt%vvel-pt%vo)**2 ) 91 zdva = SQRT( (pt%ua -pt%uo)**2 + (pt%va -pt%vo)**2 ) 92 zSs = 1.5 * SQRT( zdva ) + 0.1 * zdva ! Sea state (eqn M.A9) 92 93 93 94 ! 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 )95 zMv = MAX( 7.62e-3*zSST+1.29e-3*(zSST**2) , 0._wp ) * z1_rday ! Buoyant convection at sides (eqn M.A10) 96 zMb = MAX( 0.58*(zdvo**0.8)*(zSST+4.0)/(zL**0.2) , 0._wp ) * z1_rday ! Basal turbulent melting (eqn M.A7 ) 97 zMe = MAX( 1./12.*(zSST+2.)*zSs*(1+cos(rpi*(zIC**3))) , 0._wp ) * z1_rday ! Wave erosion (eqn M.A8 ) 97 98 98 99 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)100 zTn = MAX( zT - zMb*zdt , 0._wp ) ! new total thickness (m) 101 znVol = zTn * zW * zL ! new volume (m^3) 102 zMnew1 = (znVol/zVol) * zM ! new mass (kg) 103 zdMb = zM - zMnew1 ! mass lost to basal melting (>0) (kg) 104 ! 105 zLn = MAX( zL - zMv*zdt , 0._wp ) ! new length (m) 106 zWn = MAX( zW - zMv*zdt , 0._wp ) ! new width (m) 107 znVol = zTn * zWn * zLn ! new volume (m^3) 108 zMnew2 = (znVol/zVol) * zM ! new mass (kg) 109 zdMv = zMnew1 - zMnew2 ! mass lost to buoyant convection (>0) (kg) 110 ! 111 zLn = MAX( zLn - zMe*zdt , 0._wp ) ! new length (m) 112 zWn = MAX( zWn - zMe*zdt , 0._wp ) ! new width (m) 113 znVol = zTn * zWn * zLn ! new volume (m^3) 114 zMnew = ( znVol / zVol ) * zM ! new mass (kg) 115 zdMe = zMnew2 - zMnew ! mass lost to erosion (>0) (kg) 116 zdM = zM - zMnew ! mass lost to all erosion and melting (>0) (kg) 116 117 ! 117 118 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)119 zLn = MAX( zL -(zMv+zMe)*zdt ,0._wp ) ! (m) 120 zWn = MAX( zW -(zMv+zMe)*zdt ,0._wp ) ! (m) 121 zTn = MAX( zT - zMb *zdt ,0._wp ) ! (m) 121 122 ! 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)123 znVol = zTn*zWn*zLn ! (m^3) 124 zMnew = (znVol/zVol)*zM ! (kg) 125 zdM = zM - zMnew ! (kg) 126 zdMb = (zM/zVol) * (zW* zL ) *zMb*zdt ! approx. mass loss to basal melting (kg) 127 zdMe = (zM/zVol) * (zT*(zW+zL)) *zMe*zdt ! approx. mass lost to erosion (kg) 128 zdMv = (zM/zVol) * (zT*(zW+zL)) *zMv*zdt ! approx. mass loss to buoyant convection (kg) 128 129 ENDIF 129 130 130 131 IF( rn_bits_erosion_fraction > 0._wp ) THEN ! Bergy bits 131 132 ! 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 meters136 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/s139 dMbitsM = MIN( Mbb*zdt ,nMbits ) ! bergy bits mass lost to melting (kg)140 nMbits = nMbits-dMbitsM ! remove mass lost to bergy bits melt141 IF( Mnew == 0._wp ) THEN ! if parent berg has completely melted then142 dMbitsM = dMbitsM +nMbits ! instantly melt all the bergy bits143 nMbits = 0._wp133 zMbits = pt%mass_of_bits ! mass of bergy bits (kg) 134 zdMbitsE = rn_bits_erosion_fraction * zdMe ! change in mass of bits (kg) 135 znMbits = zMbits + zdMbitsE ! add new bergy bits to mass (kg) 136 zLbits = MIN( zL, zW, zT, 40._wp ) ! assume bergy bits are smallest dimension or 40 meters 137 zAbits = ( zMbits / rn_rho_bergs ) / zLbits ! Effective bottom area (assuming T=Lbits) 138 zMbb = MAX( 0.58*(zdvo**0.8)*(zSST+2.0)/(zLbits**0.2), 0.) * z1_rday ! Basal turbulent melting (for bits) 139 zMbb = rn_rho_bergs * zAbits * zMbb ! in kg/s 140 zdMbitsM = MIN( zMbb*zdt , znMbits ) ! bergy bits mass lost to melting (kg) 141 znMbits = znMbits-zdMbitsM ! remove mass lost to bergy bits melt 142 IF( zMnew == 0._wp ) THEN ! if parent berg has completely melted then 143 zdMbitsM = zdMbitsM + znMbits ! instantly melt all the bergy bits 144 znMbits = 0._wp 144 145 ENDIF 145 146 ELSE ! No bergy bits 146 Abits = 0._wp147 dMbitsE = 0._wp148 dMbitsM = 0._wp149 nMbits = pt%mass_of_bits ! retain previous value incase non-zero147 zAbits = 0._wp 148 zdMbitsE = 0._wp 149 zdMbitsM = 0._wp 150 znMbits = pt%mass_of_bits ! retain previous value incase non-zero 150 151 ENDIF 151 152 … … 154 155 z1_e1e2 = 1._wp / e1e2t(ii,ij) * this%mass_scaling 155 156 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 ) 157 zmelt = ( zdM - ( zdMbitsE - zdMbitsM ) ) * z1_dt ! kg/s 158 berg_grid%floating_melt(ii,ij) = berg_grid%floating_melt(ii,ij) + zmelt * z1_e1e2 ! kg/m2/s 159 zheat = zmelt * pt%heat_density ! kg/s x J/kg = J/s 160 berg_grid%calving_hflx (ii,ij) = berg_grid%calving_hflx (ii,ij) + zheat * z1_e1e2 ! W/m2 161 CALL melt_budget( ii, ij, zMnew, zheat, this%mass_scaling, & 162 & zdM, zdMbitsE, zdMbitsM, zdMb, zdMe, & 163 & zdMv, z1_dt_e1e2 ) 161 164 ELSE 162 165 WRITE(numout,*) 'thermodynamics: berg ',this%number(:),' appears to have grounded at ',narea,ii,ij … … 167 170 168 171 ! 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 ) ) THEN171 T =Tn172 Tn =Wn173 Wn =T172 zDn = ( rn_rho_bergs / pp_rho_seawater ) * zTn ! draught (keel depth) 173 IF( zDn > 0._wp .AND. MAX(zWn,zLn) < SQRT( 0.92*(zDn**2) + 58.32*zDn ) ) THEN 174 zT = zTn 175 zTn = zWn 176 zWn = zT 174 177 endif 175 178 176 179 ! Store the new state of iceberg (with L>W) 177 pt%mass = Mnew178 pt%mass_of_bits = nMbits179 pt%thickness = Tn180 pt%width = min( Wn,Ln)181 pt%length = max( Wn,Ln)180 pt%mass = zMnew 181 pt%mass_of_bits = znMbits 182 pt%thickness = zTn 183 pt%width = min(zWn,zLn) 184 pt%length = max(zWn,zLn) 182 185 183 186 next=>this%next … … 185 188 !!gm add a test to avoid over melting ? 186 189 187 IF( Mnew <= 0._wp ) THEN ! Delete the berg if completely melted190 IF( zMnew <= 0._wp ) THEN ! Delete the berg if completely melted 188 191 CALL delete_iceberg_from_list( first_berg, this ) 189 192 ! 190 193 ELSE ! Diagnose mass distribution on grid 191 194 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) 195 CALL size_budget( ii, ij, zWn, zLn, zAbits, & 196 & this%mass_scaling, zMnew, znMbits, z1_e1e2) 193 197 ENDIF 194 198 !
Note: See TracChangeset
for help on using the changeset viewer.
