| 1 |
cswdice |
| 2 |
c |
| 3 |
#include "CPP_OPTIONS.h" |
| 4 |
|
| 5 |
SUBROUTINE ICE_THERM( |
| 6 |
I i,j,bi,bj,qleft,fsalt,fresh,compact,myThid) |
| 7 |
C /==========================================================\ |
| 8 |
C | S/R ICE_MODEL | |
| 9 |
C |==========================================================| |
| 10 |
C | Calculate thermodynamic changes to ice column | |
| 11 |
C \==========================================================/ |
| 12 |
IMPLICIT NONE |
| 13 |
|
| 14 |
C == Global data == |
| 15 |
#include "SIZE.h" |
| 16 |
#include "EEPARAMS.h" |
| 17 |
#include "PARAMS.h" |
| 18 |
#include "GRID.h" |
| 19 |
#include "DYNVARS.h" |
| 20 |
#include "FFIELDS.h" |
| 21 |
#ifdef ALLOW_THERM_SEAICE |
| 22 |
#include "ICE.h" |
| 23 |
#include "BULKF_ICE_CONSTANTS.h" |
| 24 |
#endif |
| 25 |
#ifdef ALLOW_BULK_FORCE |
| 26 |
#include "BULKF.h" |
| 27 |
#endif |
| 28 |
|
| 29 |
C === Routine arguments === |
| 30 |
C myThid :: Thread no. that called this routine. |
| 31 |
INTEGER myThid |
| 32 |
INTEGER i,j, bi,bj |
| 33 |
_RL qleft ! net heat flux to ocean (W/m^2) |
| 34 |
_RL fsalt ! salt flux to ocean (kg/m^2/s) |
| 35 |
_RL fresh ! fresh water flux to ocean (kg/m^2/s) |
| 36 |
_RL compact ! fraction of grid area covered in ice |
| 37 |
|
| 38 |
c ADAPTED FROM: |
| 39 |
c LANL CICE.v2.0.2 |
| 40 |
c----------------------------------------------------------------------- |
| 41 |
c.. thermodynamics (vertical physics) based on M. Winton 3-layer model |
| 42 |
c.. See Bitz, C. M. and W. H. Lipscomb, 1999: "An energy-conserving |
| 43 |
c.. thermodynamic sea ice model for climate study." J. Geophys. |
| 44 |
c.. Res., 104, 15669 - 15677. |
| 45 |
c.. Winton, M., 1999: "A reformulated three-layer sea ice model." |
| 46 |
c.. Submitted to J. Atmos. Ocean. Technol. |
| 47 |
|
| 48 |
c.. authors Elizabeth C. Hunke and William Lipscomb |
| 49 |
c.. Fluid Dynamics Group, Los Alamos National Laboratory |
| 50 |
c----------------------------------------------------------------------- |
| 51 |
cc****subroutine thermo_winton(n,fice,fsnow,dqice,dTsfc) |
| 52 |
c.. Compute temperature change using Winton model with 2 ice layers, of |
| 53 |
c.. which only the top layer has a variable heat capacity. |
| 54 |
|
| 55 |
#ifdef ALLOW_THERM_SEAICE |
| 56 |
|
| 57 |
c |
| 58 |
c Local variables |
| 59 |
c |
| 60 |
integer k, nitmax |
| 61 |
parameter(nitmax=20) |
| 62 |
|
| 63 |
_RL area ! fraction of grid covered by ice (either 0 or 1) |
| 64 |
_RL hi ! ice height |
| 65 |
_RL hs ! snow height |
| 66 |
_RL frsnow ! fractional snow cover |
| 67 |
|
| 68 |
_RL fswabs ! net SW down at surface (W m-2) |
| 69 |
_RL fswpen ! SW penetrating beneath surface (W m-2) |
| 70 |
_RL fswdn ! SW absorbed at surface (W m-2) |
| 71 |
_RL fswint ! SW absorbed in ice (W m-2) |
| 72 |
_RL fswocn ! SW passed through ice to ocean (W m-2) |
| 73 |
_RL albedo ! surface albedo |
| 74 |
|
| 75 |
_RL flwup ! upward LW at surface (W m-2) |
| 76 |
_RL fsh ! surface downward sensible heat (W m-2) |
| 77 |
_RL flh ! surface downward latent heat (W m-2) |
| 78 |
_RL flx0 ! net surf flux, w/out conduction (W m-2) |
| 79 |
|
| 80 |
_RL df0dT ! deriv of flx0 wrt Tsf (W m-2 deg-1) |
| 81 |
|
| 82 |
_RL fct ! heat conducted to top surface |
| 83 |
_RL fcb ! heat conducted to bottom surface |
| 84 |
_RL Fbot ! oceanic heat flux |
| 85 |
_RL frzmlt ! freezing/melting potential |
| 86 |
|
| 87 |
_RL Tice(nlyr) ! internal ice temperatures |
| 88 |
|
| 89 |
_RL k12, k32 ! thermal conductivity terms |
| 90 |
_RL a10, b10 ! coefficients in quadratic eqn for T1 |
| 91 |
_RL a1, b1, c1 ! coefficients in quadratic eqn for T1 |
| 92 |
_RL Tsf ! surface (ice or snow) temperature |
| 93 |
_RL dTsf ! change in Tsf |
| 94 |
_RL Tsf_start ! old value of Tsf |
| 95 |
_RL Tf ! freezing temperature (C) |
| 96 |
|
| 97 |
_RL etop ! energy for top melting (J m-2) |
| 98 |
_RL ebot ! energy for bottom melting (J m-2) |
| 99 |
_RL hnew(nlyr) ! new ice layer thickness (m) |
| 100 |
_RL hlyr ! individual ice layer thickness (m) |
| 101 |
_RL dhi ! change in ice thickness |
| 102 |
_RL dhs ! change in snow thickness |
| 103 |
_RL rq ! rho * q for a layer |
| 104 |
_RL rqh ! rho * q * h for a layer |
| 105 |
_RL qbot ! q for new ice at bottom surf (J m-3) |
| 106 |
_RL dt ! timestep |
| 107 |
_RL esurp ! surplus energy from melting (J m-2) |
| 108 |
_RL qicen (nlyr) ! ice enthalpy (J m-3) |
| 109 |
_RL mwater0 ! fresh water mass gained/lost (kg/m^2) |
| 110 |
_RL msalt0 ! salt gained/lost (kg/m^2) |
| 111 |
|
| 112 |
_RL ustar, cphm, cpchr |
| 113 |
_RL ust, vst, evap, ssq |
| 114 |
|
| 115 |
_RL hilim, hslim, chi, chs |
| 116 |
integer iceornot |
| 117 |
|
| 118 |
_RL test |
| 119 |
|
| 120 |
dt = deltaTtracer |
| 121 |
hilim = 5. _d 0 |
| 122 |
hslim = 1.5 _d 0 |
| 123 |
|
| 124 |
hi = 0. ! initialize ice state |
| 125 |
hs = 0. |
| 126 |
Tsf = 0. |
| 127 |
do k=1,nlyr |
| 128 |
qicen(k) = 0. |
| 129 |
Tice(k) = 0. |
| 130 |
enddo |
| 131 |
flh = 0. ! initialize fluxes |
| 132 |
fsh = 0. |
| 133 |
Fbot = 0.d0 |
| 134 |
flwup = 0. |
| 135 |
fswocn = 0. |
| 136 |
esurp = 0. ! initialize energies |
| 137 |
etop = 0. |
| 138 |
ebot = 0. |
| 139 |
mwater0 = 0. ! initialize water, salt mass |
| 140 |
msalt0 = 0. |
| 141 |
frsnow = 0. |
| 142 |
chi = 0. |
| 143 |
chs = 0. |
| 144 |
flx0 = 0. |
| 145 |
df0dT = 0. |
| 146 |
evap = 0. |
| 147 |
ssq = 0. |
| 148 |
ust = 0. |
| 149 |
|
| 150 |
|
| 151 |
|
| 152 |
cBB debugging |
| 153 |
cBB print*,' ' |
| 154 |
cBB print*,'Q *** START *** i=',i,' j=',j |
| 155 |
cBB print*,'Q1 ti1,ti2,tsf,Ta', |
| 156 |
cBB & Tice1(i,j,bi,bj), Tice2(i,j,bi,bj), Tsrf(i,j,bi,bj), |
| 157 |
cBB & Tair(i,j,bi,bj) |
| 158 |
cBB print*,'Q qice1, qice2',qice1(i,j,bi,bj), qice2(i,j,bi,bj) |
| 159 |
cBB print*,'Q1 hi, hs, rain',iceHeight(i,j,bi,bj), |
| 160 |
cBB & snowHeight(i,j,bi,bj), rain(i,j,bi,bj) |
| 161 |
cBB print*,'Q1 theta, salt', theta(i,j,1,bi,bj), |
| 162 |
cBB & salt(i,j,1,bi,bj) |
| 163 |
|
| 164 |
c.. specific melting temperature for grid cell |
| 165 |
Tf = -mu_Tf*salt(i,j,1,bi,bj) ! melt temperaure |
| 166 |
|
| 167 |
cQQ from V3.0 |
| 168 |
cphm = cpwater*rhosw*delz(1) |
| 169 |
frzmlt = (Tf-theta(i,j,1,bi,bj))*cphm/dt |
| 170 |
cQQ from v2 |
| 171 |
cQQ frzmlt = (Tf-theta(i,j,1,bi,bj))*cpwater*30.d0/dt |
| 172 |
cQQQQQQQQQQQQQQQQQQQQQQQQQQQQ reduced to 1000 re v3QQQQQQQQQQQQQQQQQQQQ |
| 173 |
cQQ frzmlt = min(max(frzmlt,-1.d4),1.d4) |
| 174 |
cBB debugging |
| 175 |
cBB print*,'frzmlt, Tf',frzmlt, Tf |
| 176 |
|
| 177 |
if (iceheight(i,j,bi,bj).ge.himin) then |
| 178 |
|
| 179 |
if (frzmlt.ge. 0.) then |
| 180 |
!----------------------------------------------------------------- |
| 181 |
! freezing conditions |
| 182 |
!----------------------------------------------------------------- |
| 183 |
Fbot = 0.d0 |
| 184 |
else |
| 185 |
!----------------------------------------------------------------- |
| 186 |
! melting conditions |
| 187 |
!----------------------------------------------------------------- |
| 188 |
cSSSSSSSS try get this done |
| 189 |
cQQQQ what is this |
| 190 |
ustar = 5.e-2 !for no currents |
| 191 |
c frictional velocity between ice and water |
| 192 |
ustar = sqrt(0.00536*(uvel(i,j,1,bi,bj)**2 + |
| 193 |
& vvel(i,j,1,bi,bj)**2) ) |
| 194 |
ustar=max(5e-3,ustar) |
| 195 |
cpchr =cpwater*rhosw*transcoef |
| 196 |
Fbot = cpchr*(Tf-theta(i,j,1,bi,bj))*ustar ! < 0 |
| 197 |
Fbot = max(Fbot,frzmlt) ! frzmlt < Fbot < 0 |
| 198 |
Fbot = min(Fbot,0.d0) |
| 199 |
cBB debugging |
| 200 |
cBB print*,'Fbot,frzmlt',Fbot, frzmlt |
| 201 |
!!! use all frzmlt for standalone runs |
| 202 |
cQQQQ used in this case ?? |
| 203 |
cQQQQ Fbot = min(0.d0,frzmlt) ! frzmlt = Fbot < 0 |
| 204 |
endif |
| 205 |
|
| 206 |
hi = iceheight(i,j,bi,bj) |
| 207 |
hs = snowheight(i,j,bi,bj) |
| 208 |
area = icemask(i,j,bi,bj) |
| 209 |
Tsf = Tsrf(i,j,bi,bj) |
| 210 |
qicen(1) = qice1(i,j,bi,bj) |
| 211 |
qicen(2) = qice2(i,j,bi,bj) |
| 212 |
c qicen(2) = qicen(2) + Lfresh !!!! note !!!! only if using transp |
| 213 |
hlyr = hi / rnlyr |
| 214 |
|
| 215 |
c mass of fresh water and salt initially present in ice |
| 216 |
mwater0 = rhos*hs + rhoi*hi |
| 217 |
msalt0 = rhoi*hi*saltice |
| 218 |
|
| 219 |
c fractional snow cover |
| 220 |
frsnow = 0. |
| 221 |
if (hs .gt. 0.) frsnow = 1. |
| 222 |
c |
| 223 |
c Compute SW flux absorbed at surface and penetrating to layer 1. |
| 224 |
c |
| 225 |
call sfc_albedo(hi,hs,Tsf,sage(i,j,bi,bj),albedo) |
| 226 |
c positive for dasilva, negative for ncep |
| 227 |
fswabs = -solar(i,j,bi,bj) * (1. - albedo) |
| 228 |
fswpen = fswabs * (1. - frsnow) * i0 |
| 229 |
fswocn = fswpen * exp(-ksolar*hi) |
| 230 |
fswint = fswpen - fswocn |
| 231 |
|
| 232 |
fswdn = fswabs - fswpen |
| 233 |
c |
| 234 |
c Compute conductivity terms at layer interfaces. |
| 235 |
c |
| 236 |
k12 = 4*kice*ksnow / (ksnow*hi+4*kice*hs) |
| 237 |
k32 = 2*kice / hi |
| 238 |
c |
| 239 |
c compute ice temperatures |
| 240 |
a1 = cpice |
| 241 |
b1 = qicen(1) + (cpwater-cpice )*Tmlt1 - Lfresh |
| 242 |
c1 = Lfresh * Tmlt1 |
| 243 |
Tice(1) = 0.5 * (-b1 - sqrt(b1*b1-4.*a1*c1)) / a1 |
| 244 |
Tice(2) = (Lfresh-qicen(2)) / cpice |
| 245 |
|
| 246 |
cBB debugging |
| 247 |
if (Tice(1).gt.0.or.Tice(2).gt.0.) then |
| 248 |
write (6,*) 'BBerr Tice(1) > 0 = ',Tice(1) |
| 249 |
write (6,*) 'BBerr Tice(2) > 0 = ',Tice(2) |
| 250 |
endif |
| 251 |
c |
| 252 |
c Compute coefficients used in quadratic formula. |
| 253 |
c |
| 254 |
a10 = rhoi*cpice *hi/(2.d0*dt) + |
| 255 |
$ k32 * (4.d0*dt*k32 + rhoi*cpice *hi) |
| 256 |
$ / (6.d0*dt*k32 + rhoi*cpice *hi) |
| 257 |
b10 = -hi* |
| 258 |
$ (rhoi*cpice *Tice(1)+rhoi*Lfresh*Tmlt1/Tice(1))/(2.d0*dt) |
| 259 |
$ - k32 * (4.d0*dt*k32*Tf+rhoi*cpice *hi*Tice(2)) |
| 260 |
$ / (6.d0*dt*k32 + rhoi*cpice *hi) - fswint |
| 261 |
c1 = rhoi*Lfresh*hi*Tmlt1 / (2.d0*dt) |
| 262 |
c |
| 263 |
c Compute new surface and internal temperatures; iterate until |
| 264 |
c Tsfc converges. |
| 265 |
c |
| 266 |
c ----- begin iteration ----- |
| 267 |
do 100 k = 1,nitmax |
| 268 |
c |
| 269 |
c Save temperatures at start of iteration. |
| 270 |
Tsf_start = Tsf |
| 271 |
c |
| 272 |
|
| 273 |
c Compute top surface flux. |
| 274 |
if (hs.gt.3.d-1) then |
| 275 |
iceornot=2 |
| 276 |
else |
| 277 |
iceornot=1 |
| 278 |
endif |
| 279 |
call bulkf_formula_lanl(uwind(i,j,bi,bj), vwind(i,j,bi,bj), |
| 280 |
& wspeed(i,j,bi,bj), |
| 281 |
& Tair(i,j,bi,bj), Qair(i,j,bi,bj), cloud(i,j,bi,bj), Tsf, |
| 282 |
& flwup, flh, fsh, df0dT, ust, vst, evap, ssq, iceornot, |
| 283 |
& readwindstress) |
| 284 |
|
| 285 |
cQQ use lw data |
| 286 |
flwup=(flwup-flw(i,j,bi,bj)) |
| 287 |
|
| 288 |
flx0 = (fswdn + flwup + fsh + flh) |
| 289 |
|
| 290 |
c |
| 291 |
c Compute new top layer and surface temperatures. |
| 292 |
c If Tsfc is computed to be > 0 C, fix Tsfc = 0 and recompute T1 |
| 293 |
c with different coefficients. |
| 294 |
c |
| 295 |
a1 = a10 - k12*df0dT / (k12-df0dT) |
| 296 |
b1 = b10 - k12*(flx0-df0dT*Tsf) / (k12-df0dT) |
| 297 |
Tice(1) = -(b1 + sqrt(b1*b1-4.*a1*c1)) / (2.*a1) |
| 298 |
cBB |
| 299 |
cBB print*,'Tice(1)',Tice(1), a1, b1, c1, k12, df0dT |
| 300 |
cBB |
| 301 |
dTsf = (flx0 + k12*(Tice(1)-Tsf)) / (k12-df0dT) |
| 302 |
Tsf = Tsf + dTsf |
| 303 |
cBB |
| 304 |
cBB print*,'Tsf,dTsf',Tsf,dTsf |
| 305 |
if (Tsf .gt. 0.) then |
| 306 |
a1 = a10 + k12 |
| 307 |
b1 = b10 ! note b1 = b10 - k12*Tf0 |
| 308 |
Tice(1) = (-b1 - sqrt(b1*b1-4.*a1*c1)) / (2.*a1) |
| 309 |
cBB |
| 310 |
cBB print*,'Tsf>0, Tice(1)',Tice(1),a1,b1,c1 |
| 311 |
cBB |
| 312 |
Tsf = 0. |
| 313 |
if (hs.gt.3.d-1) then |
| 314 |
iceornot=2 |
| 315 |
else |
| 316 |
iceornot=1 |
| 317 |
endif |
| 318 |
call bulkf_formula_lanl(uwind(i,j,bi,bj), vwind(i,j,bi,bj), |
| 319 |
& wspeed(i,j,bi,bj), |
| 320 |
& Tair(i,j,bi,bj), Qair(i,j,bi,bj), cloud(i,j,bi,bj), |
| 321 |
& Tsf, |
| 322 |
& flwup, flh, fsh, df0dT, ust, vst, evap, ssq, iceornot, |
| 323 |
& readwindstress) |
| 324 |
cQQ use lw data |
| 325 |
flwup=(flwup-flw(i,j,bi,bj)) |
| 326 |
|
| 327 |
cBB |
| 328 |
cBB print*,'(b)heat flux, ',fswdn,flwup,fsh,flh |
| 329 |
flx0 = (fswdn + flwup + fsh + flh) |
| 330 |
dTsf = 0.d0 |
| 331 |
Tsf = 0.d0 |
| 332 |
endif |
| 333 |
c |
| 334 |
c Check for convergence. If no convergence, then repeat. |
| 335 |
c |
| 336 |
c Convergence test: Make sure Tsfc has converged, within prescribed error. |
| 337 |
c (Energy conservation is guaranteed within machine roundoff, even |
| 338 |
c if Tsfc has not converged.) |
| 339 |
c If no convergence, then repeat. |
| 340 |
c |
| 341 |
if (abs(dTsf).lt.Terrmax) then |
| 342 |
goto 150 |
| 343 |
elseif (k.eq.nitmax) then |
| 344 |
write (6,*) 'BB: thermw conv err ',i,j,dTsf |
| 345 |
write (6,*) 'BB: thermw conv err, iceheight ', |
| 346 |
& iceheight(i,j,bi,bj) |
| 347 |
write (6,*) 'BB: thermw conv err: Tsf, flx0', Tsf,flx0 |
| 348 |
if (Tsf.lt.-70.d0) stop !QQQQ fails |
| 349 |
endif |
| 350 |
|
| 351 |
100 continue ! surface temperature iteration |
| 352 |
150 continue |
| 353 |
c ------ end iteration ------------ |
| 354 |
|
| 355 |
cQQ |
| 356 |
if (Tsf.lt.-70.d0) then |
| 357 |
cQQ print*,'QQQQQ Tsf =',Tsf |
| 358 |
cQQ stop !QQQQ fails |
| 359 |
endif |
| 360 |
c |
| 361 |
c Compute new bottom layer temperature. |
| 362 |
c |
| 363 |
Tice(2) = (2.d0*dt*k32*(Tice(1)+2.*Tf) |
| 364 |
$ + rhoi*cpice *hi*Tice(2)) |
| 365 |
$ / (6.d0*dt*k32 + rhoi*cpice *hi) |
| 366 |
|
| 367 |
cBB |
| 368 |
cBB print*,'**********************************' |
| 369 |
cBB print*,'Q final Tice1,Tice2,Tsf,sst', |
| 370 |
cBB & Tice(1),Tice(2),Tsf,theta(i,j,1,bi,bj) |
| 371 |
|
| 372 |
c |
| 373 |
c Compute final flux values at surfaces. |
| 374 |
c |
| 375 |
fct = k12*(Tsf-Tice(1)) |
| 376 |
fcb = 4.d0*kice *(Tice(2)-Tf)/hi |
| 377 |
if (Tsf.lt.0.) flx0 = flx0 + df0dT*dTsf |
| 378 |
c |
| 379 |
c Compute new enthalpy for each layer. |
| 380 |
c |
| 381 |
qicen(1) = -cpwater*Tmlt1 + cpice *(Tmlt1-Tice(1)) + |
| 382 |
$ Lfresh*(1.-Tmlt1/Tice(1)) |
| 383 |
qicen(2) = -cpice *Tice(2) + Lfresh |
| 384 |
c |
| 385 |
cBB debugging |
| 386 |
C Make sure internal ice temperatures do not exceed Tmlt. |
| 387 |
c (This should not happen for reasonable values of i0.) |
| 388 |
c |
| 389 |
if (Tice(1) .ge. Tmlt1) then |
| 390 |
write (6,*) 'BBerr - Bug: IceT(1) > Tmlt', |
| 391 |
& i,j,Tice(1),Tmlt1 |
| 392 |
end if |
| 393 |
if (Tice(2) .ge. 0.0) then |
| 394 |
write (6,*) 'BBerr - Bug: IceT(2) > 0', |
| 395 |
& i,j,Tice(2) |
| 396 |
end if |
| 397 |
|
| 398 |
|
| 399 |
c ******************************************************************* |
| 400 |
c...................................................................... |
| 401 |
c |
| 402 |
c.. Compute growth and/or melting at the top and bottom surfaces....... |
| 403 |
c...................................................................... |
| 404 |
c |
| 405 |
c Compute enthalpy of new ice growing at bottom surface. |
| 406 |
c |
| 407 |
qbot = -cpice *Tf + Lfresh |
| 408 |
c |
| 409 |
c Compute energy available for melting/growth. |
| 410 |
c |
| 411 |
if (Tsf .eq. 0.) etop = (flx0 - fct) * dt |
| 412 |
ebot = (fcb-Fbot) * dt |
| 413 |
cBBB |
| 414 |
cBBB if (ebot.gt.0) |
| 415 |
cBBB & print*, 'fcb,fbot,ebot',fcb,fbot,ebot,frzmlt |
| 416 |
c |
| 417 |
c Initialize layer thicknesses. |
| 418 |
c Make sure internal ice temperatures do not exceed Tmlt. |
| 419 |
c If they do, then eliminate the layer. (Dont think this will happen |
| 420 |
c for reasonable values of i0.) |
| 421 |
c |
| 422 |
do k = 1, nlyr |
| 423 |
hnew(k) = hlyr |
| 424 |
enddo |
| 425 |
c |
| 426 |
c Top melt: snow, then ice. |
| 427 |
c |
| 428 |
if (etop .gt. 0.) then |
| 429 |
if (hs. gt. 0.) then |
| 430 |
rq = rhos * qsnow |
| 431 |
rqh = rq * hs |
| 432 |
if (etop .lt. rqh) then |
| 433 |
hs = hs - etop/rq |
| 434 |
etop = 0. |
| 435 |
else |
| 436 |
etop = etop - rqh |
| 437 |
hs = 0. |
| 438 |
endif |
| 439 |
endif |
| 440 |
|
| 441 |
do k = 1, nlyr |
| 442 |
if (etop .gt. 0.) then |
| 443 |
rq = rhoi * qicen(k) |
| 444 |
rqh = rq * hnew(k) |
| 445 |
if (etop .lt. rqh) then |
| 446 |
hnew(k) = hnew(k) - etop / rq |
| 447 |
etop = 0. |
| 448 |
else |
| 449 |
etop = etop - rqh |
| 450 |
hnew(k) = 0. |
| 451 |
endif |
| 452 |
endif |
| 453 |
enddo |
| 454 |
endif |
| 455 |
c |
| 456 |
c If ice is gone and melting energy remains |
| 457 |
c |
| 458 |
if (etop .gt. 0.) then |
| 459 |
write (6,*) 'All ice melts from top ', i,j |
| 460 |
go to 200 |
| 461 |
endif |
| 462 |
c |
| 463 |
c Bottom melt/growth. |
| 464 |
c |
| 465 |
if (ebot .lt. 0.) then |
| 466 |
dhi = -ebot / (qbot * rhoi) |
| 467 |
ebot = 0. |
| 468 |
k = nlyr |
| 469 |
qicen(k) = (hnew(k)*qicen(k)+dhi*qbot) / (hnew(k)+dhi) |
| 470 |
hnew(k) = hnew(k) + dhi |
| 471 |
else |
| 472 |
do k = nlyr, 1, -1 |
| 473 |
if (ebot .gt. 0. .and. hnew(k) .gt. 0.) then |
| 474 |
rq = rhoi * qicen(k) |
| 475 |
rqh = rq * hnew(k) |
| 476 |
if (ebot .lt. rqh) then |
| 477 |
hnew(k) = hnew(k) - ebot / rq |
| 478 |
ebot = 0. |
| 479 |
else |
| 480 |
ebot = ebot - rqh |
| 481 |
hnew(k) = 0. |
| 482 |
endif |
| 483 |
endif |
| 484 |
enddo |
| 485 |
c |
| 486 |
c If ice melts completely and snow is left, remove the snow with |
| 487 |
c energy from the mixed layer |
| 488 |
c |
| 489 |
if (ebot .gt. 0. .and. hs .gt. 0.) then |
| 490 |
rq = rhos * qsnow |
| 491 |
rqh = rq * hs |
| 492 |
if (ebot .lt. rqh) then |
| 493 |
hs = hs - ebot / rq |
| 494 |
esurp = esurp - rhos*qsnow*hs |
| 495 |
ebot = 0. |
| 496 |
else |
| 497 |
ebot = ebot - rqh |
| 498 |
hs = 0. |
| 499 |
endif |
| 500 |
endif |
| 501 |
|
| 502 |
if (ebot .gt. 0.) then |
| 503 |
cBB write (6,*) 'All ice melts from bottom', i,j |
| 504 |
go to 200 |
| 505 |
endif |
| 506 |
|
| 507 |
endif |
| 508 |
c |
| 509 |
c Compute new total ice thickness. |
| 510 |
c |
| 511 |
hi = 0. |
| 512 |
do k = 1, nlyr |
| 513 |
hi = hi + hnew(k) |
| 514 |
enddo |
| 515 |
|
| 516 |
c If hi < himin, melt the ice. |
| 517 |
if (hi .lt. himin) then |
| 518 |
esurp = esurp - rhos*qsnow*hs |
| 519 |
do k = 1, nlyr |
| 520 |
esurp = esurp - rhoi*qicen(k)*hnew(k) |
| 521 |
enddo |
| 522 |
hi = 0. |
| 523 |
hs = 0. |
| 524 |
Tsf=0. |
| 525 |
do k = 1, nlyr |
| 526 |
qicen(k) = 0. |
| 527 |
enddo |
| 528 |
go to 200 |
| 529 |
endif |
| 530 |
|
| 531 |
c |
| 532 |
c Let it snow |
| 533 |
c |
| 534 |
cQQm check |
| 535 |
frsnow = 0. ! temporarily indicates whether fresh snow accumulates |
| 536 |
if (Tsf.lt.0.and.Tair(i,j,bi,bj).le.Tf0kel.and.hi.gt.0.d0) then |
| 537 |
frsnow = 1. |
| 538 |
cQQ change if reading in snow -- including left over value of rain |
| 539 |
sage(i,j,bi,bj)=(dt+sage(i,j,bi,bj)*(1-dt/(50*24*60*60))) |
| 540 |
& *exp(rain(i,j,bi,bj)) |
| 541 |
snow(i,j,bi,bj)=-rain(i,j,bi,bj)*1000/rhos |
| 542 |
rain(i,j,bi,bj)=0.0 |
| 543 |
hs = hs + snow(i,j,bi,bj)*dt |
| 544 |
else |
| 545 |
snow(i,j,bi,bj)=0.0 |
| 546 |
sage(i,j,bi,bj)=(dt+sage(i,j,bi,bj)*(1-dt/(50*24*60*60))) |
| 547 |
endif |
| 548 |
cBB |
| 549 |
if (sage(i,j,bi,bj).gt.50*24*60*60) then |
| 550 |
print*,'sage',i,j,sage(i,j,bi,bj) |
| 551 |
endif |
| 552 |
|
| 553 |
cswd QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ |
| 554 |
c If ice evap is used to sublimate surface snow/ice or |
| 555 |
c if no ice pass on to ocean |
| 556 |
if (hs.gt.0.d0) then |
| 557 |
if (-evap*1000/rhos *dt.gt.hs) then |
| 558 |
evap=evap+hs*rhos/1000/dt |
| 559 |
hs=0.d0 |
| 560 |
else |
| 561 |
hs = hs + evap*1000/rhos *dt |
| 562 |
evap=0.d0 |
| 563 |
endif |
| 564 |
endif |
| 565 |
if (hi.gt.0.d0.and.evap.lt.0.d0) then |
| 566 |
if (-evap*1000/rhoi *dt.gt.hi) then |
| 567 |
evap=evap+hi*rhoi/1000/dt |
| 568 |
hi=0.d0 |
| 569 |
else |
| 570 |
hi = hi + evap*1000/rhoi *dt |
| 571 |
evap=0.d0 |
| 572 |
endif |
| 573 |
endif |
| 574 |
|
| 575 |
|
| 576 |
c If there is enough snow to lower the ice/snow interface below |
| 577 |
c freeboard, convert enough snow to ice to bring the interface back |
| 578 |
c to sea-level. Adjust enthalpy of top ice layer accordingly. |
| 579 |
c |
| 580 |
if ( hs .gt. hi*rhoiw/rhos ) then |
| 581 |
cBB write (6,*) 'Freeboard adjusts' |
| 582 |
dhi = (hs * rhos - hi * rhoiw) / rhosw |
| 583 |
dhs = dhi * rhoi / rhos |
| 584 |
rqh = rhoi*qicen(1)*hnew(1) + rhos*qsnow*dhs |
| 585 |
hnew(1) = hnew(1) + dhi |
| 586 |
qicen(1) = rqh / (rhoi*hnew(1)) |
| 587 |
hi = hi + dhi |
| 588 |
hs = hs - dhs |
| 589 |
end if |
| 590 |
|
| 591 |
cswd |
| 592 |
c If ice evap is used to sublimate surface snow/ice or |
| 593 |
c if no ice pass on to ocean |
| 594 |
|
| 595 |
|
| 596 |
cQQQQQQ |
| 597 |
c if (hi.gt.hilim) then |
| 598 |
c print*,'BBerr, hi>hilim',i,j,hi |
| 599 |
c chi=hi-hilim |
| 600 |
c hi=hilim |
| 601 |
c endif |
| 602 |
c if (hs.gt.hslim) then |
| 603 |
c print*,'BBerr, hs>hslim',i,j,hs |
| 604 |
c chs=hs-hslim |
| 605 |
c hs=hslim |
| 606 |
c endif |
| 607 |
cQQQQQQ |
| 608 |
|
| 609 |
|
| 610 |
hlyr = hi/rnlyr |
| 611 |
call new_layers_winton(hs,hi,hlyr,hnew,qicen) |
| 612 |
|
| 613 |
200 continue |
| 614 |
c |
| 615 |
c Compute surplus energy left over from melting. |
| 616 |
c |
| 617 |
esurp = esurp + etop + ebot |
| 618 |
c.. heat fluxes left over for ocean |
| 619 |
cQQ qleft = -fswocn |
| 620 |
cQQ if (frzmlt.le.0.) qleft =qleft -(Fbot+esurp/dt) ! W/m^2 |
| 621 |
else ! 0 < hi < himin |
| 622 |
c If hi < himin, melt the ice. |
| 623 |
esurp = esurp - rhos*qsnow*snowheight(i,j,bi,bj) |
| 624 |
do k = 1, nlyr |
| 625 |
esurp = esurp - rhoi*qicen(k)*iceheight(i,j,bi,bj)*0.5 |
| 626 |
qicen(k) = 0. |
| 627 |
enddo |
| 628 |
hi = 0. |
| 629 |
hs = 0. |
| 630 |
Tsf = 0. |
| 631 |
cQQ call bulkf_formula_lanl(uwind(i,j,bi,bj), vwind(i,j,bi,bj), |
| 632 |
cQQ & wspeed(i,j,bi,bj), |
| 633 |
cQQ & Tair(i,j,bi,bj), Qair(i,j,bi,bj), cloud(i,j,bi,bj), Tsf, |
| 634 |
cQQ & flwup, flh, fsh, df0dT, ust, vst, evap, ssq, iceornot, |
| 635 |
cQQ & readwindstress) |
| 636 |
cQQ use lw data |
| 637 |
cQQ flwup=(flwup-flw(i,j,bi,bj)) |
| 638 |
cQQ albedo=0.1 !QQQQ need to unhardwire |
| 639 |
cQQ flx0 = (-solar(i,j,bi,bj)*(1.d0-albedo) + flwup + fsh + flh) |
| 640 |
c.. heat fluxes left over for ocean |
| 641 |
cQQ qleft=-(flx0+esurp/dt) ! W/m^2 |
| 642 |
endif ! hi |
| 643 |
c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 644 |
c.. fresh and salt fluxes |
| 645 |
c QQQQQ |
| 646 |
fresh=0.0 |
| 647 |
fsalt=0.0 |
| 648 |
fresh = (mwater0 - (rhos*(hs+chs) + rhoi*(hi+chi)))/dt+evap |
| 649 |
fsalt = (msalt0 - rhoi*hi*saltice)/35.d0/dt ! for same units as fresh |
| 650 |
c.. heat fluxes left over for ocean |
| 651 |
qleft = -fswocn |
| 652 |
if (frzmlt.le.0.) qleft =qleft -(Fbot+esurp/dt) ! W/m^2 |
| 653 |
cBB debug |
| 654 |
cBB print*,'Q qnet, qleft, Fbot,esurp/dt', qnet(i,j,bi,bj), |
| 655 |
cBB & qleft, Fbot, esurp/dt |
| 656 |
cQQQQQQQQQQQQ |
| 657 |
c if (Tsf.lt.-50.d0) then |
| 658 |
c Tsf=-50.d0 |
| 659 |
c endif |
| 660 |
cQQ if (Tsf.lt.Tair(i,j,bi,bj)-Tf0kel-3.d0) then |
| 661 |
cQQ Tsf=Tair(i,j,bi,bj)-Tf0kel-3.d0 |
| 662 |
cQQ endif |
| 663 |
cQQQQQQ |
| 664 |
if (hi.gt.hilim) then |
| 665 |
cBBB print*,'BBerr, hi>hilim',i,j,hi |
| 666 |
chi=hi-hilim |
| 667 |
fresh=fresh+rhoi*chi/dt |
| 668 |
hi=hilim |
| 669 |
endif |
| 670 |
if (hs.gt.hslim) then |
| 671 |
cBB print*,'BBerr, hs>hslim',i,j,hs |
| 672 |
chs=hs-hslim |
| 673 |
fresh=fresh+rhos*chs/dt |
| 674 |
hs=hslim |
| 675 |
endif |
| 676 |
cQQQQQQ |
| 677 |
|
| 678 |
|
| 679 |
c |
| 680 |
c .. store new values |
| 681 |
iceheight(i,j,bi,bj)=hi |
| 682 |
c if (hi.gt.0) then |
| 683 |
c icemask(i,j,bi,bj)=1 |
| 684 |
c else |
| 685 |
c icemask(i,j,bi,bj)=0 |
| 686 |
c endif |
| 687 |
if (hi.gt.0.d0) then |
| 688 |
compact=iceMask(i,j,bi,bj) |
| 689 |
else |
| 690 |
compact=0.d0 |
| 691 |
endif |
| 692 |
snowheight(i,j,bi,bj)=hs |
| 693 |
tsrf(i,j,bi,bj)=Tsf |
| 694 |
tice1(i,j,bi,bj)=Tice(1) |
| 695 |
tice2(i,j,bi,bj)=Tice(2) |
| 696 |
qice1(i,j,bi,bj)=qicen(1) |
| 697 |
qice2(i,j,bi,bj)=qicen(2) |
| 698 |
|
| 699 |
cBB debugging |
| 700 |
cBB print*,'Q4 i,j,hi,hs,Tsf',i,j,hi,hs,Tsf |
| 701 |
cBB print*,'Q4 qleft,fsalt,fresh', qleft,fsalt,fresh |
| 702 |
|
| 703 |
#endif |
| 704 |
|
| 705 |
return |
| 706 |
end ! main thermo routine |
| 707 |
c -------------------------------------------------------------------- |
| 708 |
|
| 709 |
#include "CPP_OPTIONS.h" |
| 710 |
|
| 711 |
c -------------------------------------------------------------------- |
| 712 |
|
| 713 |
CStartofinterface |
| 714 |
subroutine new_layers_winton(hs,hi,hlyr,hnew,qicen) |
| 715 |
c |
| 716 |
c Repartition into equal-thickness layers, conserving energy. |
| 717 |
c This is the 2-layer version. |
| 718 |
c |
| 719 |
implicit none |
| 720 |
C == Global data == |
| 721 |
#include "SIZE.h" |
| 722 |
#include "EEPARAMS.h" |
| 723 |
#include "PARAMS.h" |
| 724 |
#include "GRID.h" |
| 725 |
#include "DYNVARS.h" |
| 726 |
#include "FFIELDS.h" |
| 727 |
#ifdef ALLOW_THERM_SEAICE |
| 728 |
#include "ICE.h" |
| 729 |
#include "BULKF_ICE_CONSTANTS.h" |
| 730 |
#endif |
| 731 |
|
| 732 |
#ifndef ALLOW_THERM_SEAICE |
| 733 |
integer nlyr |
| 734 |
parameter(nlyr=1) |
| 735 |
#endif |
| 736 |
|
| 737 |
_RL hnew(nlyr) ! new ice layer thickness (m) |
| 738 |
_RL hs |
| 739 |
_RL hi |
| 740 |
_RL hlyr ! individual ice layer thickness (m) |
| 741 |
_RL qicen (nlyr) ! ice enthalpy (J m-3) |
| 742 |
c |
| 743 |
#ifdef ALLOW_THERM_SEAICE |
| 744 |
c Local variables |
| 745 |
_RL f1 ! Fraction of upper layer ice in new layer |
| 746 |
_RL qh1, qh2 ! qice*h for layers 1 and 2 |
| 747 |
_RL qhtot ! qh1 + qh2 |
| 748 |
_RL q2tmp ! Temporary value of qice for layer 2 |
| 749 |
|
| 750 |
if (hnew(1).gt.hnew(2)) then ! Layer 1 gives ice to layer 2 |
| 751 |
f1 = (hnew(1)-hlyr)/hlyr |
| 752 |
q2tmp = f1*qicen(1) + (1.-f1)*qicen(2) |
| 753 |
if (q2tmp.gt.Lfresh) then |
| 754 |
qicen(2) = q2tmp |
| 755 |
else ! Keep q2 fixed to avoid q2<Lfresh and T2>0 |
| 756 |
qh2 = hlyr*qicen(2) |
| 757 |
qhtot = hnew(1)*qicen(1) + hnew(2)*qicen(2) |
| 758 |
qh1 = qhtot - qh2 |
| 759 |
qicen(1) = qh1/hlyr |
| 760 |
endif |
| 761 |
else ! Layer 2 gives ice to layer 1 |
| 762 |
f1 = hnew(1)/hlyr |
| 763 |
qicen(1) = f1*qicen(1) + (1.-f1)*qicen(2) |
| 764 |
endif |
| 765 |
|
| 766 |
#endif |
| 767 |
|
| 768 |
return |
| 769 |
end |
| 770 |
|
| 771 |
c -------------------------------------------------------------------- |
| 772 |
CStartofinterface |
| 773 |
subroutine sfc_albedo(hi,hs,Tsf,age,albedo) |
| 774 |
c.. Compute surface albedo |
| 775 |
implicit none |
| 776 |
C == Global data == |
| 777 |
#include "SIZE.h" |
| 778 |
#include "EEPARAMS.h" |
| 779 |
#include "PARAMS.h" |
| 780 |
#include "GRID.h" |
| 781 |
#include "DYNVARS.h" |
| 782 |
#include "FFIELDS.h" |
| 783 |
#ifdef ALLOW_THERM_SEAICE |
| 784 |
#include "ICE.h" |
| 785 |
#include "BULKF_ICE_CONSTANTS.h" |
| 786 |
#endif |
| 787 |
|
| 788 |
_RL frsnow ! fractional snow cover |
| 789 |
_RL hi ! ice height |
| 790 |
_RL hs ! snow height |
| 791 |
_RL Tsf ! surface temperature |
| 792 |
_RL albedo ! surface albedo |
| 793 |
|
| 794 |
_RL frsnalb, salb |
| 795 |
_RL albsno, albice, age |
| 796 |
|
| 797 |
#ifdef ALLOW_THERM_SEAICE |
| 798 |
c frsnow = 0. |
| 799 |
c if (hs .gt. 0.) frsnow = 1. |
| 800 |
|
| 801 |
c if (Tsf .lt. 0.) then |
| 802 |
c albedo = frsnow*albsnodry + |
| 803 |
c $ (1.-frsnow)*(albicemin + (albicemax - albicemin) |
| 804 |
c $ *(1.-exp(-hi/halb))) |
| 805 |
c else |
| 806 |
c albedo = frsnow*albsnowet + |
| 807 |
c $ (1.-frsnow)*(albicemin + (albicemax - albicemin) |
| 808 |
c $ *(1.-exp(-hi/halb))) |
| 809 |
c endif |
| 810 |
c GISS model albedo calculation |
| 811 |
if (Tsf.lt.-10.d0) then |
| 812 |
frsnalb=0.3 |
| 813 |
elseif (Tsf.le.0.d0) then |
| 814 |
frsnalb=0.3-0.015*(Tsf+10) |
| 815 |
else |
| 816 |
frsnalb=0.15 |
| 817 |
endif |
| 818 |
|
| 819 |
albice=0.7 |
| 820 |
albsno=0.55+frsnalb*exp(-0.2*age/(24*60*60)) !make age units into days |
| 821 |
albedo=albice*exp(-hs/0.3)+albsno*(1-exp(-hs/0.3)) |
| 822 |
|
| 823 |
if (albedo.gt.1.d0.or.albedo.lt..5d0) then |
| 824 |
print*,'QQ - albedo problem', albedo, age, hs, albsno |
| 825 |
stop |
| 826 |
endif |
| 827 |
#endif |
| 828 |
|
| 829 |
return |
| 830 |
end |
| 831 |
c -------------------------------------------------------------------- |
Parent Directory
| 
Revision Log
| 
Revision Graph