! Generated by TAPENADE (INRIA, Tropics team)
! Version 2.1 - (Id: 1.33 vmp Stable - Wed Jan 12 13:35:07 MET 2005)
!
! $Id: ice_atmo.F,v 1.10 2004/02/09 17:58:46 lipscomb Exp $
!=======================================================================
!BOP
!
! !MODULE: ice_atmo - atm-ice interface: stability based flux calculations
!
! !DESCRIPTION:
!
! Atmospheric boundary interface (stability based flux calculations)
!
! !REVISION HISTORY:
!
! author: Elizabeth C. Hunke, LANL
!
! Vectorized by Clifford Chen (Fujitsu) and William H. Lipscomb (LANL)
!
! !INTERFACE:
!
MODULE
ICE_ATMO
USE
ice_constantsUSE
ice_domainUSE
ice_fluxUSE
ice_gridUSE
ice_kinds_modUSE
ice_model_sizeUSE
ice_stateREAL(KIND=DBL_KIND)
::umin
CONTAINS
SUBROUTINE
ATMO_BOUNDARY_LAYER
(n
,sfctype
,tsf
,strx
,stry
,trf
,qrf
,&
&
delt
,delq
)IMPLICIT NONE
!ech & - c2*atan(xd) + 1.571_dbl_kind
!
!
! !DESCRIPTION:
!
! Compute coefficients for atm/ice fluxes, stress, and reference
! temperature. NOTE: \\
! (1) all fluxes are positive downward, \\
! (2) net heat flux = fswabs + flwup + (1-emissivity)flwdn + fsh + flh, \\
! (3) here, tstar = (WT)/U*, and qstar = (WQ)/U*, \\
! (4) wind speeds should all be above a minimum speed (eg. 1.0 m/s). \\
!
! ASSUME: \\
! The saturation humidity of air at T(K): qsat(T) (kg/m**3) \\
!
! Code originally based on CSM1 \\
!
! !REVISION HISTORY:
!
! author: Elizabeth C. Hunke, LANL
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
! thickness category index
!
! ice or ocean
!
! surface temperature of ice or ocean
!
! x surface stress (N)
! y surface stress (N)
! reference height temperature (K)
! reference height specific humidity (kg/kg)
! potential T difference (K)
! humidity difference (kg/kg)
!
!EOP
!
! iteration index
! horizontal indices
!
! surface temperature in Kelvin (K)
! temporary variable
! stability function at zlvl (momentum)
! stress at zlvl
! interpolation factor
! ln(z10 /zTrf)
! stability function at zTrf (heat and water)
! stable profile
! sat surface humidity (kg/kg)
! for qsat, dqsatdt
! for qsat, dqsatdt
! the saturation humidity of air (kg/m^3)
! Lvap or Lsub, depending on surface type
!
! ustar (m/s)
! tstar
! qstar
! sqrt of neutral exchange coefficient (momentum)
! sqrt of neutral exchange coefficient (heat)
! sqrt of neutral exchange coefficient (water)
! sqrt of exchange coefficient (momentum)
! sqrt of exchange coefficient (water)
! sqrt of exchange coefficient (heat)
! surface wind magnitude (m/s)
! ln(zlvl /z10)
! virtual temperature (K)
! specific heat of moist air
! H (at zlvl ) over L
! stability factor
! stability function at zlvl (heat and water)
!
! compressed index in i-direction
! compressed index in j-direction
!
! number of cells with aicen > puny or tmask true
! combined ij index
!
! defined as cp_wv/cp_air - 1.
! reference height for air temp (m)
! local functions
!jkim &, umin = c1 ! minimum wind speed (m/s)
!
! dummy argument
! unstable part of psimh
! unstable part of psimx
!------------------------------------------------------------
! Define functions
!------------------------------------------------------------
!
!
INCLUDE
'DIFFSIZES.inc
'! Hint: nbdirsmax should be the maximum number of differentiation directions
REAL(KIND=DBL_KIND)
, DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT)
::delq
REAL(KIND=DBL_KIND)
, DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT)
::delt
INTEGER(KIND=INT_KIND),INTENT(IN)
::n
REAL(KIND=DBL_KIND)
, DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT)
::qrf
CHARACTER
(len
=3
),INTENT(IN)
::sfctype
REAL(KIND=DBL_KIND)
, DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT)
::strx
REAL(KIND=DBL_KIND)
, DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT)
::stry
REAL(KIND=DBL_KIND)
, DIMENSION(ilo:ihi, jlo:jhi),INTENT(OUT)
::trf
REAL(KIND=DBL_KIND)
, DIMENSION(ilo:ihi, jlo:jhi),INTENT(IN)
::tsf
REAL(KIND=DBL_KIND)
,PARAMETER
::cpvir
=cp_wv/cp_air-c1REAL(KIND=DBL_KIND)
,PARAMETER
::ztrf
=c2REAL
::abs1
REAL
::abs2
REAL(KIND=DBL_KIND)
::arg1
INTEGER(KIND=INT_KIND)
::icells
,ij
INTEGER(KIND=INT_KIND)
:: indxi((ihi-ilo+1
)*(jhi-jlo+1
)), indxj((ihi&
&
-ilo+1
)*(jhi-jlo+1
))INTEGER(KIND=INT_KIND)
::i
,j
,k
REAL
::min1
,x1
REAL(KIND=DBL_KIND)
::psimhu
,psixhu
,xd
REAL(KIND=DBL_KIND)
::al2
,fac
,lheat
,psimh
,psimhs
,psix2
,qqq
,&
&
qsat
,ssq
,tau
,tsfk
,ttt
,xqq
REAL(KIND=DBL_KIND)
:: alz((ihi-ilo+1
)*(jhi-jlo+1
)), cp((ihi-ilo+1
)*&
&
(jhi-jlo+1
)), hol((ihi-ilo+1
)*(jhi-jlo+1
)), psixh((ihi-ilo+1
)*(jhi-&
&
jlo+1
)), qstar((ihi-ilo+1
)*(jhi-jlo+1
)), rd((ihi-ilo+1
)*(jhi-jlo+1
)&
&
), rdn((ihi-ilo+1
)*(jhi-jlo+1
)), re((ihi-ilo+1
)*(jhi-jlo+1
)), ren((&
&
ihi-ilo+1
)*(jhi-jlo+1
)), rh((ihi-ilo+1
)*(jhi-jlo+1
)), rhn((ihi-ilo+&
&
1
)*(jhi-jlo+1
)), stable((ihi-ilo+1
)*(jhi-jlo+1
)), thva((ihi-ilo+1
)*&
&
(jhi-jlo+1
)), tstar((ihi-ilo+1
)*(jhi-jlo+1
)), ustar((ihi-ilo+1
)*(&
&
jhi-jlo+1
)), vmag((ihi-ilo+1
)*(jhi-jlo+1
))REAL
::x2
REAL
::y1
INTRINSIC
EXP
,MAX
,SIGN
,ABS
,LOG
,ATAN
,MIN
,SQRT
!
!------------------------------------------------------------
! Initialize
!------------------------------------------------------------
arg1 = zref/ztrf al2 =LOG
(arg1)!
!
DO
j=jlo,jhiDO
i=ilo,ihi lhcoef(i, j) = c0 shcoef(i, j) = c0 strx(i, j) = c0 stry(i, j) = c0 trf(i, j) = c0 qrf(i, j) = c0 delt(i, j) = c0 delq(i, j) = c0END DO
END DO
!------------------------------------------------------------
! Compute turbulent flux coefficients, wind stress, and
! reference temperature and humidity.
!------------------------------------------------------------
! sfctype
!
!
IF
(sfctype(1
:3
) .EQ.'ice'
)THEN
!------------------------------------------------------------
! identify grid cells with ice
!------------------------------------------------------------
!
icells =0
DO
j=jlo,jhiDO
i=ilo,ihiIF
(aicen(i, j, n) .GT. puny)THEN
icells = icells +1
indxi(icells) = i indxj(icells) = jEND IF
END DO
END DO
!------------------------------------------------------------
! define some needed variables
!------------------------------------------------------------
!
! for qsat
qqq = qqqice! for qsat
ttt = tttice! ice to vapor
lheat = lsubDO
ij=1
,icells i = indxi(ij) j = indxj(ij)IF
(umin .LT. wind(i, j))THEN
vmag(ij) = wind(i, j)ELSE
vmag(ij) = uminEND IF
! neutral coefficient
arg1 = zref/iceruf rdn(ij) = vonkar/LOG
(arg1)END DO
ELSE
IF
(sfctype(1
:3
) .EQ.'ocn'
)THEN
! ij
!
!------------------------------------------------------------
! identify ocean cells
!------------------------------------------------------------
!
icells =0
DO
j=jlo,jhiDO
i=ilo,ihiIF
(tmask(i, j))THEN
icells = icells +1
indxi(icells) = i indxj(icells) = jEND IF
END DO
END DO
!------------------------------------------------------------
! define some needed variables
!------------------------------------------------------------
!
qqq = qqqocn ttt = tttocn! liquid to vapor
lheat = lvapDO
ij=1
,icells i = indxi(ij) j = indxj(ij)IF
(umin .LT. wind(i, j))THEN
vmag(ij) = wind(i, j)ELSE
vmag(ij) = uminEND IF
arg1 =0.0027_dbl_kind
/vmag(ij) +.000142_dbl_kind
+&
&
.0000764_dbl_kind
*vmag(ij) rdn(ij) =SQRT
(arg1)END DO
END IF
! ij
!
!
!
DO
ij=1
,icells i = indxi(ij)!------------------------------------------------------------
! define some more needed variables
!------------------------------------------------------------
j = indxj(ij)!
!
! surface temp (K)
tsfk = tsf(i, j) + tffresh! saturation humidity (kg/m^3)
arg1 = -(ttt/tsfk) qsat = qqq*EXP
(arg1)! sat surf hum (kg/kg)
ssq = qsat/rhoa(i, j)!
! virtual pot temp (K)
thva(ij) = pott(i, j)*(c1+zvir*qa(i, j))! pot temp diff (K)
delt(i, j) = pott(i, j) - tsfk! spec hum dif (kg/kg)
delq(i, j) = qa(i, j) - ssq arg1 = zlvl(i, j)/zref alz(ij) =LOG
(arg1)!------------------------------------------------------------
! first estimate of Z/L and ustar, tstar and qstar
!------------------------------------------------------------
! neutral coefficients, z/L = 0.0
cp(ij) = cp_air*(c1+cpvir*ssq)!
!
rhn(ij) = rdn(ij)! ustar,tstar,qstar
ren(ij) = rdn(ij)!
ustar(ij) = rdn(ij)*vmag(ij) tstar(ij) = rhn(ij)*delt(i, j) qstar(ij) = ren(ij)*delq(i, j)END DO
!
! ij
!------------------------------------------------------------
! iterate to converge on Z/L, ustar, tstar and qstar
!------------------------------------------------------------
!
!
! do k=1,2 !!! ech thinks this should be at least 5
DO
k=1
,5
!
DO
ij=1
,icells i = indxi(ij)! compute stability & evaluate all stability functions
j = indxj(ij)!
hol(ij) = vonkar*gravit*zlvl(i, j)*(tstar(ij)/thva(ij)+qstar(ij)&
&
/(c1/zvir+qa(i, j)))/ustar(ij)**2
IF
(hol(ij) .GE.0.
)THEN
x1 = hol(ij)ELSE
x1 = -hol(ij)END IF
IF
(x1 .GT. c10)THEN
min1 = c10ELSE
min1 = x1END IF
hol(ij) =SIGN
(min1, hol(ij)) stable(ij) = p5 +SIGN
(p5, hol(ij))IF
(c1 - c16*hol(ij) .GE.0.
)THEN
abs1 = c1 - c16*hol(ij)ELSE
abs1 = -(c1-c16*hol(ij))END IF
x2 =SQRT
(abs1)IF
(x2 .LT. c1)THEN
xqq = c1ELSE
xqq = x2END IF
xqq =SQRT
(xqq)!ccsm psimh = -c5*hol(ij)*stable(ij) + (c1-stable(ij))*psimhu(xqq)
!ccsm psixh(ij) = -c5*hol(ij)*stable(ij) + (c1-stable(ij))*psixhu(xqq)
!
! Jordan et al 1999
arg1 = -(0.35_dbl_kind
*hol(ij)) psimhs = -(0.7_dbl_kind
*hol(ij)+0.75_dbl_kind
*(hol(ij)-&
&
14.3_dbl_kind
)*EXP
(arg1)+10.7_dbl_kind
)!
arg1 = (c1+xqq*(c2+xqq))*(c1+xqq*xqq)/c8 psimh = psimhs*stable(ij) + (c1-stable(ij))*(LOG
(arg1)-c2*ATAN
(&
&
xqq)+pi*p5)! shift all coeffs to measurement height and stability
arg1 = (c1+xqq*xqq)/c2 psixh(ij) = psimhs*stable(ij) + (c1-stable(ij))*(c2*LOG
(arg1))!
rd(ij) = rdn(ij)/(c1+rdn(ij)/vonkar*(alz(ij)-psimh)) rh(ij) = rhn(ij)/(c1+rhn(ij)/vonkar*(alz(ij)-psixh(ij)))! update ustar, tstar, qstar using updated, shifted coeffs
re(ij) = ren(ij)/(c1+ren(ij)/vonkar*(alz(ij)-psixh(ij)))!
ustar(ij) = rd(ij)*vmag(ij) tstar(ij) = rh(ij)*delt(i, j) qstar(ij) = re(ij)*delq(i, j)END DO
END DO
!
! ij
! end iteration
!
!
DO
ij=1
,icells i = indxi(ij)!------------------------------------------------------------
! coefficients for turbulent flux calculation
!------------------------------------------------------------
! add windless coefficient for sensible heat flux
! as in Jordan et al (JGR, 1999)
j = indxj(ij)!
!
!ccsm shcoef(i,j) = rhoa(i,j) * ustar(ij) * cp(ij) * rh(ij)
shcoef(i, j) = rhoa(i, j)*ustar(ij)*cp(ij)*rh(ij) + c1!------------------------------------------------------------
! momentum flux
!------------------------------------------------------------
! tau = rhoa(i,j) * ustar * ustar
! strx = tau * uatm(i,j) / vmag
! stry = tau * vatm(i,j) / vmag
!------------------------------------------------------------
lhcoef(i, j) = rhoa(i, j)*ustar(ij)*lheat*re(ij)!
!
! not the stress at zlvl(i,j)
tau = rhoa(i, j)*ustar(ij)*rd(ij) strx(i, j) = tau*uatm(i, j)!------------------------------------------------------------
! Compute diagnostics: 2m ref T & Q
!------------------------------------------------------------
stry(i, j) = tau*vatm(i, j)!
hol(ij) = hol(ij)*ztrf/zlvl(i, j)IF
(c1 - c16*hol(ij) .GE.0.
)THEN
abs2 = c1 - c16*hol(ij)ELSE
abs2 = -(c1-c16*hol(ij))END IF
y1 =SQRT
(abs2)IF
(c1 .LT. y1)THEN
xqq = y1ELSE
xqq = c1END IF
xqq =SQRT
(xqq) arg1 = (c1+xqq*xqq)/c2 psix2 = -(c5*hol(ij)*stable(ij)) + (c1-stable(ij))*(c2*LOG
(arg1)) fac = rh(ij)/vonkar*(alz(ij)+al2-psixh(ij)+psix2) trf(i, j) = pott(i, j) - delt(i, j)*fac! pot temp to temp correction
trf(i, j) = trf(i, j) -0.01
*ztrf fac = re(ij)/vonkar*(alz(ij)+al2-psixh(ij)+psix2) qrf(i, j) = qa(i, j) - delq(i, j)*facEND DO
END
SUBROUTINE
ATMO_BOUNDARY_LAYER
END MODULE
ICE_ATMO