SUBROUTINE SFC_SICE(IM,KM,PS,U1,V1,T1,Q1, & HICE ,FICE ,TICE, SFCDSW, & SHELEG,SNWDPH,TSKIN,QSURF,TPRCP,SRFLAG,STC,DM, & DLWFLX,SLRAD,SNOWMT,DELT,GFLUX,CM,CH, & RCL,PRSL1,PRSLKI,SLIMSK, + CMM,CHH,ZLVL, Clu_q2m_iter [-1L/+1L]: add flag_iter !* & EVAP,HFLX,EP,DDVEL) & EVAP,HFLX,EP,DDVEL,flag_iter,MOM4ICE,lsm) ! USE MACHINE , ONLY : kind_phys USE FUNCPHYS, ONLY : fpvs USE PHYSCONS, SBC => con_sbc, HVAP => con_HVAP, TGICE => con_tice &, CP => con_CP, HFUS => con_HFUS, JCAL => con_JCAL &, EPS => con_eps, EPSM1 => con_epsm1, grav => con_g &, RVRDM1 => con_FVirt, T0C => con_T0C, RD => con_RD implicit none ! ! include 'constant.h' ! integer IM, km, kmi, lsm ! real(kind=kind_phys), parameter :: cpinv=1.0/cp, HVAPI=1.0/HVAP real(kind=kind_phys) DELT real(kind=kind_phys) PS(IM), U1(IM), V1(IM), & T1(IM), Q1(IM), SHELEG(IM), & TSKIN(IM), QSURF(IM), STC(IM,KM), & DM(IM), DLWFLX(IM), SLRAD(IM), & SNOWMT(IM), GFLUX(IM), & CM(IM), CH(IM), & RCL(IM), PRSL1(IM), PRSLKI(IM), & SLIMSK(IM), EVAP(IM), HFLX(IM), & RNET(IM), EP(IM), DDVEL(IM) &, TPRCP(IM), SRFLAG(IM) &, SFCDSW(IM), FICE(IM), HICE(IM) &, CMM(IM),CHH(IM),ZLVL(IM) &, SNWDPH(IM) real(kind=kind_phys) TICE(IM), FFW(IM), EVAPI(IM), & EVAPW(IM), HFLXI(IM), HFLXW(IM), & SNETI(IM), SNETW(IM), QSSI(IM), & QSSW(IM) real(kind=kind_phys) hf(IM), hfd(IM), hfi(IM), & hfw(IM), focn(IM), snof(IM) real(kind=kind_phys) hi_save(IM), hs_save(IM) Clu_q2m_iter [+1L]: add flag_iter logical flag_iter(im) ! ! Locals ! integer k,i ! real(kind=kind_phys) & PSURF(IM), Q0(IM), QS1(IM), & QSS(IM), RCAP(IM), RCH(IM), & RHO(IM), SLWD(IM), SNET(IM), & SNOEVP(IM), SNOWD(IM),THETA1(IM), & TSURF(IM), TV1(IM), XRCL(IM), & PS1(IM), WIND(IM) ! real(kind=kind_phys) & elocp, sigma, cimin, himin, & himax, hsmax, timin, & t12, t14, tem real(kind=kind_phys) albfw, emissiv ! cc PARAMETER (kmi=2) ! 2-layer of ice PARAMETER (sigma=sbc) PARAMETER (ELOCP=HVAP/CP) ! PARAMETER (CIMIN=0.15) ! minimum ice concentration required PARAMETER (HIMAX=8.0) ! maximum ice thickness allowed PARAMETER (HIMIN=0.1) ! minimum ice thickness required PARAMETER (HSMAX=2.) ! maximum snow depth allowed PARAMETER (TIMIN=173.) ! minimum temperature allowed for snow/ice ! PARAMETER (CICE=1880.*917.) PARAMETER (ALBFW=0.06) ! Albedo for lead ! PARAMETER (EMISSIV=0.95) ! emissivity of snow and ice PARAMETER (EMISSIV=1.0) ! emissivity of snow and ice real, PARAMETER :: DSI=1.0/0.33 ! LOGICAL FLAG(IM), FLAGSNW(IM) LOGICAL MOM4ICE real(kind=kind_phys) STSICE(IM,KMI) IF (MOM4ICE) then CIMIN=0.15 ! MOM4ICE and MASK ELSE CIMIN=0.50 ! GFS ONLY ENDIF ! C C FLAG for sea-ice C DO I = 1, IM FLAG(I) = SLIMSK(I).GE.2. .AND. flag_iter(i) ENDDO IF (MOM4ICE) then DO I = 1, IM IF(FLAG(I)) THEN HI_save(I)=HICE(I) HS_save(I)=SHELEG(I)*0.001 ENDIF ENDDO ENDIF DO I = 1, IM IF (flag_iter(i).AND.SLIMSK(I).LT.1.5) THEN HICE(I )= 0. FICE(I )= 0. ENDIF ENDDO C C snow-rain detection C if (.not. mom4ice .and. lsm > 0) then DO I=1,IM IF(FLAG(I)) THEN IF(SRFLAG(I) .EQ. 1.) THEN EP(I) = 0. SHELEG(i) = SHELEG(i) + 1.E3*TPRCP(i) TPRCP(i) = 0. ENDIF ENDIF ENDDO endif C C INITIALIZE VARIABLES. ALL UNITS ARE SUPPOSEDLY M.K.S. UNLESS SPECIFIE C PSURF IS IN PASCALS C WIND IS WIND SPEED, THETA1 IS ADIABATIC SURFACE TEMP FROM LEVEL 1 C RHO IS DENSITY, QS1 IS SAT. HUM. AT LEVEL1 AND QSS IS SAT. HUM. AT C SURFACE C CONVERT SLRAD TO THE CIVILIZED UNIT FROM LANGLEY MINUTE-1 K-4 C ! qs1 = fpvs(t1) ! qss = fpvs(tskin) DO I=1,IM IF(FLAG(I)) THEN XRCL(I) = SQRT(RCL(I)) PSURF(I) = 1000. * PS(I) PS1(I) = 1000. * PRSL1(I) SLWD(I) = SLRAD(I) c c DLWFLX has been given a negative sign for downward longwave c snet is the net shortwave flux c SNET(I) = -SLWD(I) - DLWFLX(I) WIND(I) = XRCL(I) * SQRT(U1(I) * U1(I) + V1(I) * V1(I)) & + MAX(0.0, MIN(DDVEL(I), 30.0)) WIND(I) = MAX(WIND(I),1.) Q0(I) = MAX(Q1(I),1.E-8) TSURF(I) = TSKIN(I) THETA1(I) = T1(I) * PRSLKI(I) TV1(I) = T1(I) * (1. + RVRDM1 * Q0(I)) RHO(I) = PS1(I) / (RD * TV1(I)) cjfe QS1(I) = 1000. * FPVS(T1(I)) qs1(i) = fpvs(t1(i)) QS1(I) = EPS * QS1(I) / (PS1(I) + EPSM1 * QS1(I)) QS1(I) = MAX(QS1(I), 1.E-8) Q0(I) = min(QS1(I),Q0(I)) cjfe QSS(I) = 1000. * FPVS(TSURF(I)) ! qss(i) = fpvs(tskin(i)) ! QSS(I) = EPS * QSS(I) / (PSURF(I) + EPSM1 * QSS(I)) FFW(I) = 1.0-FICE(I) ! IF (FICE(I) .GE. cimin) THEN ! TICE(I) = (TSKIN(I)-TGICE*FFW(I))/FICE(I) ! ELSE IF (FICE(I) .LT. cimin) THEN PRINT *,'WARNING: ice fraction is low:', FICE(I) FICE(I )= cimin FFW(I) = 1.0-FICE(I) TICE(I) = tgice TSKIN(I)= tgice PRINT *,'Fix ice fraction: reset it to:', FICE(I) ENDIF qssi(i) = fpvs(tice(i)) QSSI(I) = EPS * QSSI(I) / (PSURF(I) + EPSM1 * QSSI(I)) qssw(i) = fpvs(tgice) QSSW(I) = EPS * QSSW(I) / (PSURF(I) + EPSM1 * QSSW(I)) C C SNOW DEPTH IN WATER EQUIVALENT IS CONVERTED FROM MM TO M UNIT C IF (MOM4ICE) then SNOWD(I) = SHELEG(I) * 0.001 / FICE(I) ELSE SNOWD(I) = SHELEG(I) / 1000. ENDIF FLAGSNW(I) = .FALSE. C C WHEN SNOW DEPTH IS LESS THAN 1 MM, A PATCHY SNOW IS ASSUMED AND C SOIL IS ALLOWED TO INTERACT WITH THE ATMOSPHERE. C WE SHOULD EVENTUALLY MOVE TO A LINEAR COMBINATION OF SOIL AND C SNOW UNDER THE CONDITION OF PATCHY SNOW. C IF(SNOWD(I).GT..001) FLAGSNW(I) = .TRUE. ENDIF ENDDO !! C C RCP = RHO CP CH V C DO I = 1, IM IF(FLAG(I)) THEN RCH(I) = RHO(I) * CP * CH(I) * WIND(I) Cwei added 10/24/2006 CMM(I)=CM(I)* WIND(I) CHH(I)=RHO(I)*CH(I)* WIND(I) ZLVL(I) = -RD * TV1(I) * LOG(PS1(I)/PSURF(I)) / GRAV ENDIF ENDDO C C SENSIBLE AND LATENT HEAT FLUX OVER OPEN WATER & SEA ICE C DO I = 1, IM IF(FLAG(I)) THEN EVAPI(I) = ELOCP * RCH(I) * (QSSI(I) - Q0(I)) EVAPW(I) = ELOCP * RCH(I) * (QSSW(I) - Q0(I)) ! EVAP(I) = FICE(I)*EVAPI(I) + FFW(I)*EVAPW(I) ENDIF ENDDO C C UPDATE SEA ICE TEMPERATURE C DO K = 1, KMI DO I=1,IM IF(FLAG(I)) THEN STSICE(I,K) = STC(I,K) ENDIF ENDDO ENDDO DO I=1,IM IF(FLAG(I)) THEN SNETW(I) = SFCDSW(I)*(1.0-ALBFW) SNETW(I) = min(3.*SNET(I)/(1.+2.*FFW(I)),SNETW(I)) SNETI(I) = (SNET(I)-FFW(I)*SNETW(I))/FICE(I) t12=tice(i)*tice(i) t14=t12*t12 C hfi = Net non-solar and upIR heat flux @ ice surface hfi(i) =-dlwflx(i)+emissiv*sigma*t14 + evapi(i) & +rch(i)*(tice(i)-theta1(i)) C hfd = Heat flux derivat @ surface C = 4 emissivity T**3 + RCH [ 1 + L/CP * DQS/DT) ] hfd(i) = 4.*emissiv*sigma*tice(i)*t12 & +(1.+elocp*eps*hvap*qs1(i)/(rd*t12))*rch(i) t12=tgice*tgice t14=t12*t12 C hfw = Net heat flux @ water surface (within ice) hfw(i) =-dlwflx(i)+emissiv*sigma*t14 + evapw(i) & +rch(i)*(tgice-theta1(i))-snetw(i) focn(i) = 2. ! heat flux from ocean - should be from ocn model snof(i) = 0. ! snowfall rate - snow accumulates in gbphys CCC...QC if (hice(i) .GT. himax) hice(i)=himax if (hice(i) .LT. himin) hice(i)=himin if (snowd(i) .GT. hsmax) snowd(i)=hsmax if (snowd(i) .GT. (2.*hice(i))) then print *,'WARNING: too much snow :',snowd(i) snowd(i)=2.*hice(i) print *,'FIX: decrease snow depth to:',snowd(i) endif CCC...QC END ENDIF ENDDO call ice3lay(IM,kmi,snowd,hice,fice,flag, & stsice,tice,hfi,sneti,hfd,hfw,focn, & snof,snowmt,gflux,delt) IF (MOM4ICE) then DO I = 1, IM IF(FLAG(I)) THEN HICE(I) = HI_save(I) SNOWD(I) = HS_save(I) ENDIF ENDDO ENDIF DO I=1,IM IF(FLAG(I)) THEN if (tice(i).LT.timin) then print *,'WARNING: snow/ice temperature is too low:',tice(i) tice(i)=timin PRINT *,'Fix snow/ice temperature: reset it to:',tice(i) endif if (stsice(i,1).LT.timin) then print *,'WARNING: Layer 1 ice temp is too low:',stsice(i,1) stsice(i,1)=timin PRINT *,'Fix Layer 1 ice temp: reset it to:',stsice(i,1) endif if (stsice(i,2).LT.timin) then print *,'WARNING: Layer 2 ice temp is too low:',stsice(i,2) stsice(i,2)=timin PRINT *,'Fix Layer 2 ice temp: reset it to:',stsice(i,2) endif tskin(i) = tice(i)*fice(i) + tgice*ffw(i) ENDIF ENDDO DO k=1,KMI DO I=1,IM IF(FLAG(I)) THEN stc(i,k) = min(stsice(i,k),T0C) ENDIF ENDDO ENDDO C C CALCULATE SENSIBLE HEAT FLUX (& EVAP over SEA ICE) C DO I = 1, IM IF (FLAG(I)) THEN HFLXI(I) = RCH(I) * (TICE(I) - THETA1(I)) HFLXW(I) = RCH(I) * (TGICE - THETA1(I)) HFLX(I) = FICE(I)*HFLXI(I) + FFW(I)*HFLXW(I) EVAP(I) = FICE(I)*EVAPI(I) + FFW(I)*EVAPW(I) ENDIF ENDDO C C THE REST OF THE OUTPUT C DO I = 1, IM IF(FLAG(I)) THEN QSURF(I) = Q1(I) + EVAP(I) / (ELOCP * RCH(I)) DM(I) = 1. C C CONVERT SNOW DEPTH BACK TO MM OF WATER EQUIVALENT C SHELEG(I) = SNOWD(I) * 1000. SNWDPH(I) = SHELEG(I) * DSI ! Snow depth in mm ENDIF ENDDO ! do i=1,im IF(FLAG(I)) THEN tem = 1.0 / rho(i) hflx(i) = hflx(i) * tem * cpinv evap(i) = evap(i) * tem * hvapi ENDIF enddo ! RETURN END SUBROUTINE ICE3LAY(IM,kmi,hs,hi,fice,flag, & TI,TS,HF,SOL,HFD,FRZMLT,FB, & SNOW,SNOWMT,GFLUX,DT) C C************************************************************************** C * C THREE-LAYER SEA ICE VERTICAL THERMODYNAMICS * C * C Based on: M. Winton, "A reformulated three-layer sea ice model", * C Journal of Atmospheric and Oceanic Technology, 2000 * C * C * C -> +---------+ <- ts - diagnostic surface temperature ( <= 0C ) * C / | | * C hs | snow | <- 0-heat capacity snow layer * C \ | | * C => +---------+ * C / | | * C / | | <- t1 - upper 1/2 ice temperature; this layer has * C / | | a variable (T/S dependent) heat capacity * C hi |...ice...| * C \ | | * C \ | | <- t2 - lower 1/2 ice temp. (fixed heat capacity) * C \ | | * C -> +---------+ <- base of ice fixed at seawater freezing temp. * C * C************************************************************************** C USE MACHINE , ONLY : kind_phys implicit none C C *************************************************************** C * C Argument Description Units Changed? * C -------- ----------- ----- -------- * C * C INPUT/OUTPUT: * C * C hs Snow Thickness m y * C hi Ice Thickness m y * C ts Surface Temperature deg C y * C ti(1) Temp @ Midpt of Ice1 deg C y * C ti(2) Temp @ Midpt of Ice2 deg C y * C hf(+) Net non-solar and upIR * C heat flux @ surface watt/m^2 n * C sol(+) Net solar incoming top watt/m^2 n * C hfd(+) Heat flux derivat @ sfc watt/(m^2deg-C) n * C fb(+) Heat Flux from Ocean watt/m^2 n * C snow Snowfall Rate m/sec n * C dt timestep sec n * C * C ADDITIONAL: * C * C snowmt snow melt during dt m y * C gflux conductive heat flux watt/m^2 y * C flag Ice mask flag n * C * C LOCAL: * C * C hdi ice-water interface m y * C hsni snow-ice m y * C *************************************************************** C integer IM, kmi real (kind=kind_phys) KS, DS, I0, KI, DI, CI, DICI, & LI, SI, MU, TFI, TFW, DILI, DSLI real (kind=kind_phys) DW,tffresh,TFI0 real (kind=kind_phys) ki4, dt, dt2, dt4, dt6 real (kind=kind_phys) r0,r1,r2,r4,p5 C C variables for temperature calculation [see Winton (2000) 2.a] C real (kind=kind_phys) A(IM), B(IM), Ip(IM), & A1(IM),B1(IM),C1(IM), & A10(IM),B10(IM), & K12(IM),K32(IM), & TSF(IM),SNOWD(IM), & H1(IM),H2(IM),DH(IM), & F1(IM),TMELT(IM),BMELT(IM) real (kind=kind_phys) HF(IM),HFD(IM), & HS(IM),TS(IM),SOL(IM), & GFLUX(IM),SNOWMT(IM), & FB(IM),SNOW(IM), & HI(IM),TI(IM,KMI) real (kind=kind_phys) FICE(IM),FRZMLT(IM) real (kind=kind_phys) HSNI(IM),HDI(IM) LOGICAL FLAG(IM) integer I C properties of ice, snow, and seawater (local) PARAMETER (DS=330.0) ! snow (over sea ice) density - 330 kg/(m^3) PARAMETER (DW=1000.0) ! fresh water density - 1000 kg/(m^3) PARAMETER (TFFRESH=273.15) ! Freezing temp of fresh ice (K) PARAMETER (KS = 0.31) ! conductivity of snow - 0.31 W/(mK) PARAMETER (I0 = 0.3) ! ice surface penetrating solar fraction PARAMETER (KI = 2.03) ! conductivity of ice - 2.03 W/(mK) PARAMETER (DI = 917.0) ! density of ice - 917 kg/(m^3) PARAMETER (CI = 2054.0) ! heat capacity of fresh ice - 2054 J/(kg K) PARAMETER (LI = 3.34e5) ! latent heat of fusion - 334e3 J/(kg-ice) PARAMETER (SI = 1.0) ! salinity of sea ice PARAMETER (MU = 0.054) ! relates freezing temp. to salinity PARAMETER (TFI = -MU*SI) ! sea ice freezing temp. = -mu*salinity PARAMETER (TFW = -1.8) ! TFW - seawater freezing temperature -1.8 C PARAMETER (r0 = 0.) ! r0=0 PARAMETER (r1 = 1.) ! r1=1 PARAMETER (r2 = 2.) ! r2=2 PARAMETER (r4 = 4.) ! r4=4 PARAMETER (p5 = r1/r2) ! p5=1/2 PARAMETER (TFI0 = TFI-0.0001) ! TFI-0.0001 PARAMETER (DICI = DI*CI) PARAMETER (DILI = DI*LI) PARAMETER (DSLI = DS*LI) PARAMETER (KI4 = KI*r4) dt2 = dt*r2 dt4 = dt*r4 dt6 = dt*6. DO I=1,IM IF(FLAG(I)) THEN hs(i) = hs(i)*DW/DS hdi(i) = (DS*hs(i)+DI*hi(i))/DW if (hi(i).LT.hdi(i)) then hs(i)=hs(i)+hi(i)-hdi(i) hsni(i)=(hdi(i)-hi(i))*DS/DI hi(i)=hi(i)+hsni(i) endif snow(i) = snow(i)*DW/DS ts(i) = TS(I)-TFFRESH ! degC ti(i,1) = min(TI(I,1)-TFFRESH,TFI0) ! degC ti(i,2) = min(TI(I,2)-TFFRESH,TFI0) ! degC Ip(i) = I0*sol(i) ! Ip +v (in Winton Ip=-I0*sol as sol -v) if (hs(i) .gt. r0) then tsf(i) = r0 Ip(i) = r0 else tsf(i) = TFI Ip(i) = I0*sol(i) ! Ip +v here (in Winton Ip=-I0*sol) endif ts(i) = min(TS(I),TSF(I)) C C Compute ice temperature C B(i) = hfd(i) A(i) = hf(i)-sol(i)+Ip(i)-ts(i)*B(i) ! +v sol input here K12(i) = KI4*KS/(KS*hi(i)+KI4*hs(i)) K32(i) = r2*KI/hi(i) A10(i) = DICI*hi(i)/dt2 + K32(i)*(dt4*K32(i)+DICI*hi(i)) & /(dt6*K32(i)+DICI*hi(i)) B10(i) = -DI*hi(i)*(CI*ti(i,1)+LI*TFI/ti(i,1))/dt2 - Ip(i) & -K32(i)*(dt4*K32(i)*TFW+DICI*hi(i)*ti(i,2)) & /(dt6*K32(i)+DICI*hi(i)) A1(i) = A10(i)+K12(i)*B(i)/(K12(i)+B(i)) B1(i) = B10(i)+A(i)*K12(i)/(K12(i)+B(i)) C1(i) = DILI*TFI/dt2*hi(i) ti(i,1) = -(sqrt(B1(i)*B1(i)-r4*A1(i)*C1(i))+B1(i))/(r2*A1(i)) ts(i) = (K12(i)*ti(i,1)-A(i))/(K12(i)+B(i)) if (ts(i) .gt. tsf(i)) then A1(i) = A10(i)+K12(i) B1(i) = B10(i)-K12(i)*tsf(i) ti(i,1) = -(sqrt(B1(i)*B1(i)-r4*A1(i)*C1(i))+B1(i))/(r2*A1(i)) ts(i) = tsf(i) tmelt(i) = (K12(i)*(ti(i,1)-tsf(i))-(A(i)+B(i)*tsf(i)))*dt else tmelt(i) = r0 hs(i) = hs(i) + snow(i)*dt endif ti(i,2) = (dt2*K32(i)*(ti(i,1)+r2*TFW)+DICI*hi(i)*ti(i,2)) & /(dt6*K32(i)+DICI*hi(i)) bmelt(i) = (fb(i)+KI4*(ti(i,2)-TFW)/hi(i))*dt C C Resize the ice ... C h1(i) = p5*hi(i) h2(i) = p5*hi(i) C C ... top ... C if (tmelt(i) .le. hs(i)*DSLI) then snowmt(i) = tmelt(i)/DSLI hs(i) = hs(i) - tmelt(i)/DSLI else snowmt(i) = hs(i) h1(i)=h1(i) & -(tmelt(i)-hs(i)*DSLI)/(DI*(CI-LI/ti(i,1))*(TFI-ti(i,1))) hs(i) = r0 endif C C ... and bottom C if (bmelt(i) .lt. r0) then dh(i) = -bmelt(i)/(DILI+DICI*(TFI-TFW)) ti(i,2) = (h2(i)*ti(i,2)+dh(i)*TFW)/(h2(i)+dh(i)) h2(i) = h2(i)+dh(i) else h2(i) = h2(i)-bmelt(i)/(DILI+DICI*(TFI-ti(i,2))) endif C C if ice remains, even up 2 layers, else, pass negative energy back in snow C hi(i) = h1(i) + h2(i) if (hi(i) .gt. r0) then if (h1(i) .gt. p5*hi(i)) then f1(i) = r1-r2*h2(i)/hi(i) ti(i,2)=f1(i)*(ti(i,1)+LI*TFI/(CI*ti(i,1))) & +(r1-f1(i))*ti(i,2) if (ti(i,2).GT.TFI) then hi(i)=hi(i)-h2(i)*CI*(TI(i,2)-TFI)/(LI*dt) ti(i,2)=TFI endif else f1(i) = r2*h1(i)/hi(i) ti(i,1)=f1(i)*(ti(i,1)+LI*TFI/(CI*ti(i,1))) & +(r1-f1(i))*ti(i,2) ti(i,1) = (ti(i,1)-sqrt(ti(i,1)*ti(i,1)-r4*TFI*LI/CI))/r2 endif K12(i) = KI4*KS/(KS*hi(i)+KI4*hs(i)) gflux(i) = k12(i) * (ti(i,1) - ts(i)) else hs(i) = hs(i) + (h1(i)*(CI*(ti(i,1)-TFI)-LI*(r1-TFI/ti(i,1))) & +h2(i)*(CI*(ti(i,2)-TFI)-LI))/LI hi(i) = max(r0,hs(i)*DS/DI) hs(i) = r0 ti(i,1) = TFW ti(i,2) = TFW gflux(i) = r0 endif gflux(i) = fice(i)*gflux(i) snowmt(i)=snowmt(i)*DS/DW hs(i)=hs(i)*DS/DW ts(i)=ts(i)+tffresh ti(i,1)=ti(i,1)+tffresh ti(i,2)=ti(i,2)+tffresh ENDIF ! FLAG ENDDO ! IM return end