phys/module_sf_temfsfclay.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2 !
   3 MODULE module_sf_temfsfclay
   4
   5 CONTAINS
   6
   7 !-------------------------------------------------------------------
   8    SUBROUTINE temfsfclay(u3d,v3d,th3d,qv3d,p3d,pi3d,rho,z,ht,         &
   9                      cp,g,rovcp,r,xlv,psfc,chs,chs2,cqs2,cpm,      &
  10                      znt,ust,mavail,xland,                         &
  11                      hfx,qfx,lh,tsk,flhc,flqc,qgh,qsfc,            &
  12                      u10,v10,th2,t2,q2,                            &
  13                      svp1,svp2,svp3,svpt0,ep1,ep2,                 &
  14                      karman,fCor,te_temf,                          &
  15                      hd_temf,exch_temf,wm_temf,                           &
  16                      ids,ide, jds,jde, kds,kde,                    &
  17                      ims,ime, jms,jme, kms,kme,                    &
  18                      its,ite, jts,jte, kts,kte                    &
  19                      )
  20 !-------------------------------------------------------------------
  21       IMPLICIT NONE
  22 !-------------------------------------------------------------------
  23 !
  24 ! This is the Total Energy - Mass Flux (TEMF) surface layer scheme.
  25 ! Initial implementation 2010 by Wayne Angevine, CIRES/NOAA ESRL.
  26 ! References:
  27 !    Angevine et al., 2010, MWR
  28 !    Angevine, 2005, JAM
  29 !    Mauritsen et al., 2007, JAS
  30 !
  31 !-------------------------------------------------------------------
  32 !-------------------------------------------------------------------
  33 !-- u3d         3D u-velocity interpolated to theta points (m/s)
  34 !-- v3d         3D v-velocity interpolated to theta points (m/s)
  35 !-- th3d        potential temperature (K)
  36 !-- qv3d        3D water vapor mixing ratio (Kg/Kg)
  37 !-- p3d         3D pressure (Pa)
  38 !-- cp          heat capacity at constant pressure for dry air (J/kg/K)
  39 !-- g           acceleration due to gravity (m/s^2)
  40 !-- rovcp       R/CP
  41 !-- r           gas constant for dry air (J/kg/K)
  42 !-- xlv         latent heat of vaporization for water (J/kg)
  43 !-- psfc        surface pressure (Pa)
  44 !-- chs         heat/moisture exchange coefficient for LSM (m/s)
  45 !-- chs2
  46 !-- cqs2
  47 !-- cpm
  48 !-- znt         roughness length (m)
  49 !-- ust         u* in similarity theory (m/s)
  50 !-- mavail      surface moisture availability (between 0 and 1)
  51 !-- xland       land mask (1 for land, 2 for water)
  52 !-- hfx         upward heat flux at the surface (W/m^2)
  53 !-- qfx         upward moisture flux at the surface (kg/m^2/s)
  54 !-- lh          net upward latent heat flux at surface (W/m^2)
  55 !-- tsk         surface temperature (K)
  56 !-- flhc        exchange coefficient for heat (W/m^2/K)
  57 !-- flqc        exchange coefficient for moisture (kg/m^2/s)
  58 !-- qgh         lowest-level saturated mixing ratio
  59 !-- qsfc        ground saturated mixing ratio
  60 !-- u10         diagnostic 10m u wind
  61 !-- v10         diagnostic 10m v wind
  62 !-- th2         diagnostic 2m theta (K)
  63 !-- t2          diagnostic 2m temperature (K)
  64 !-- q2          diagnostic 2m mixing ratio (kg/kg)
  65 !-- svp1        constant for saturation vapor pressure (kPa)
  66 !-- svp2        constant for saturation vapor pressure (dimensionless)
  67 !-- svp3        constant for saturation vapor pressure (K)
  68 !-- svpt0       constant for saturation vapor pressure (K)
  69 !-- ep1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
  70 !-- ep2         constant for specific humidity calculation
  71 !               (R_d/R_v) (dimensionless)
  72 !-- karman      Von Karman constant
  73 !-- fCor        Coriolis parameter
  74 !-- ids         start index for i in domain
  75 !-- ide         end index for i in domain
  76 !-- jds         start index for j in domain
  77 !-- jde         end index for j in domain
  78 !-- kds         start index for k in domain
  79 !-- kde         end index for k in domain
  80 !-- ims         start index for i in memory
  81 !-- ime         end index for i in memory
  82 !-- jms         start index for j in memory
  83 !-- jme         end index for j in memory
  84 !-- kms         start index for k in memory
  85 !-- kme         end index for k in memory
  86 !-- its         start index for i in tile
  87 !-- ite         end index for i in tile
  88 !-- jts         start index for j in tile
  89 !-- jte         end index for j in tile
  90 !-- kts         start index for k in tile
  91 !-- kte         end index for k in tile
  92 !-------------------------------------------------------------------
  93       INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
  94                                         ims,ime, jms,jme, kms,kme, &
  95                                         its,ite, jts,jte, kts,kte
  96 !
  97       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
  98                 INTENT(IN   ) :: u3d, v3d, th3d, qv3d, p3d, pi3d, rho, z
  99       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 100                 INTENT(IN   ) :: mavail, xland, tsk, fCor, ht, psfc, znt
 101       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 102                 INTENT(INOUT) :: hfx, qfx, lh, flhc, flqc
 103       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 104                 INTENT(INOUT) :: ust, chs2, cqs2, chs, cpm, qgh, qsfc
 105       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 106                 INTENT(OUT  ) :: u10, v10, th2, t2, q2
 107       REAL,     DIMENSION( ims:ime, jms:jme )           , &
 108                 INTENT(IN   ) :: hd_temf
 109       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 110                 INTENT(INOUT) :: te_temf
 111       REAL,     DIMENSION( ims:ime, jms:jme )           , &
 112                 INTENT(  OUT) :: exch_temf
 113       REAL,     DIMENSION( ims:ime, jms:jme )           , &
 114                 INTENT(INOUT) :: wm_temf
 115
 116
 117       REAL,     INTENT(IN   ) :: cp,g,rovcp,r,xlv
 118       REAL,     INTENT(IN   ) :: svp1,svp2,svp3,svpt0
 119       REAL,     INTENT(IN   ) :: ep1,ep2,karman
 120 !
 121 ! LOCAL VARS
 122
 123       INTEGER ::  J
 124 !
 125
 126       DO J=jts,jte
 127
 128         CALL temfsfclay1d(j,u1d=u3d(ims,kms,j),v1d=v3d(ims,kms,j),     &
 129                 th1d=th3d(ims,kms,j),qv1d=qv3d(ims,kms,j),p1d=p3d(ims,kms,j), &
 130                 pi1d=pi3d(ims,kms,j),rho=rho(ims,kms,j),z=z(ims,kms,j),&
 131                 zsrf=ht(ims,j),      &
 132                 cp=cp,g=g,rovcp=rovcp,r=r,xlv=xlv,psfc=psfc(ims,j),    &
 133                 chs=chs(ims,j),chs2=chs2(ims,j),cqs2=cqs2(ims,j),      &
 134                 cpm=cpm(ims,j),znt=znt(ims,j),ust=ust(ims,j),          &
 135                 mavail=mavail(ims,j),xland=xland(ims,j),    &
 136                 hfx=hfx(ims,j),qfx=qfx(ims,j),lh=lh(ims,j),tsk=tsk(ims,j), &
 137                 flhc=flhc(ims,j),flqc=flqc(ims,j),qgh=qgh(ims,j),      &
 138                 qsfc=qsfc(ims,j),u10=u10(ims,j),v10=v10(ims,j),        &
 139                 th2=th2(ims,j),t2=t2(ims,j),q2=q2(ims,j),        &
 140                 svp1=svp1,svp2=svp2,svp3=svp3,svpt0=svpt0,             &
 141                 ep1=ep1,ep2=ep2,karman=karman,fCor=fCor(ims,j),  &
 142                 te_temfx=te_temf(ims,kms,j),hd_temfx=hd_temf(ims,j), &
 143                 exch_temfx=exch_temf(ims,j),wm_temfx=wm_temf(ims,j), &
 144                 ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde,     &
 145                 ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme,     &
 146                 its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte      &
 147                                                                    )
 148       ENDDO
 149
 150    END SUBROUTINE temfsfclay
 151
 152
 153 !-------------------------------------------------------------------
 154    SUBROUTINE temfsfclay1d(j,u1d,v1d,th1d,qv1d,p1d, &
 155                 pi1d,rho,z,zsrf,cp,g,rovcp,r,xlv,psfc,    &
 156                 chs,chs2,cqs2,cpm,znt,ust,          &
 157                 mavail,xland,hfx,qfx,lh,tsk, &
 158                 flhc,flqc,qgh,qsfc,u10,v10,        &
 159                 th2,t2,q2,svp1,svp2,svp3,svpt0,             &
 160                 ep1,ep2,karman,fCor,  &
 161                 te_temfx,hd_temfx,exch_temfx,wm_temfx,       &
 162                 ids,ide, jds,jde, kds,kde,                    &
 163                 ims,ime, jms,jme, kms,kme,                    &
 164                 its,ite, jts,jte, kts,kte                    &
 165                      )
 166 !!-------------------------------------------------------------------
 167       IMPLICIT NONE
 168 !!-------------------------------------------------------------------
 169       INTEGER,  INTENT(IN   ) ::        ids,ide, jds,jde, kds,kde, &
 170                                         ims,ime, jms,jme, kms,kme, &
 171                                         its,ite, jts,jte, kts,kte, &
 172                                         j
 173
 174       REAL,     DIMENSION( ims:ime ), INTENT(IN   ) ::             &
 175                                         u1d,v1d,qv1d,p1d,th1d,pi1d,rho,z,zsrf
 176       REAL,     INTENT(IN   ) ::        cp,g,rovcp,r,xlv
 177       REAL,     DIMENSION( ims:ime ), INTENT(IN   ) :: psfc,znt
 178       REAL,     DIMENSION( ims:ime ), INTENT(INOUT) ::             &
 179                                         chs,chs2,cqs2,cpm,ust
 180       REAL,     DIMENSION( ims:ime ), INTENT(IN   ) :: mavail,xland
 181       REAL,     DIMENSION( ims:ime ), INTENT(INOUT) ::             &
 182                                         hfx,qfx,lh
 183       REAL,     DIMENSION( ims:ime ), INTENT(IN   ) :: tsk
 184       REAL,     DIMENSION( ims:ime ), INTENT(  OUT) ::             &
 185                                         flhc,flqc
 186       REAL,     DIMENSION( ims:ime ), INTENT(INOUT) ::             &
 187                                         qgh,qsfc
 188       REAL,     DIMENSION( ims:ime ), INTENT(  OUT) ::             &
 189                                         u10,v10,th2,t2,q2
 190       REAL,     INTENT(IN   ) ::        svp1,svp2,svp3,svpt0
 191       REAL,     INTENT(IN   ) ::        ep1,ep2,karman
 192       REAL,     DIMENSION( ims:ime ), INTENT(IN   ) :: fCor,hd_temfx
 193       REAL,     DIMENSION( ims:ime ), INTENT(INOUT) :: te_temfx
 194       REAL,     DIMENSION( ims:ime ), INTENT(  OUT) :: exch_temfx, wm_temfx
 195 !
 196 !! LOCAL VARS
 197 ! TE model constants
 198    real, parameter :: visc_temf = 1.57e-5
 199    real, parameter :: conduc_temf = 1.57e-5 / 0.733
 200    logical, parameter :: MFopt = .true.  ! Use mass flux or not
 201    real, parameter :: TEmin = 1e-3
 202    real, parameter :: ftau0 = 0.17
 203    real, parameter :: fth0 = 0.145
 204 !   real, parameter :: fth0 = 0.12    ! WA 10/13/10 to make PrT0 ~= 1
 205    real, parameter :: Cf = 0.185
 206    real, parameter :: CN = 2.0
 207 !   real, parameter :: Ceps = ftau0**1.5
 208    real, parameter :: Ceps = 0.070
 209    real, parameter :: Cgamma = Ceps
 210    real, parameter :: Cphi = Ceps
 211 !   real, parameter :: PrT0 = Cphi/Ceps * ftau0**2. / 2 / fth0**2.
 212    real, parameter :: PrT0 = Cphi/Ceps * ftau0**2 / 2. / fth0**2
 213 !
 214    integer :: i
 215    real :: e1
 216    real, dimension( its:ite)    ::  wstr, ang, wm
 217    real, dimension( its:ite)    ::  z0t
 218    real, dimension( its:ite) :: dthdz, dqtdz, dudz, dvdz
 219    real, dimension( its:ite) :: lepsmin
 220    real, dimension( its:ite) :: thetav
 221    real, dimension( its:ite) :: zt,zm
 222    real, dimension( its:ite) :: N2, S, Ri, beta, ftau, fth, ratio
 223    real, dimension( its:ite) :: TKE, TE2
 224    real, dimension( its:ite) :: ustrtilde, linv, leps
 225    real, dimension( its:ite) :: km, kh
 226    real, dimension( its:ite) :: qsfc_air
 227 !!-------------------------------------------------------------------
 228
 229 !!!!!!! ******
 230 ! WA Known outages:  None
 231
 232    do i = its,ite      ! Main loop
 233
 234       ! Calculate surface saturated q and q in air at surface
 235       e1=svp1*exp(svp2*(tsk(i)-svpt0)/(tsk(i)-svp3))
 236       qsfc(i)=ep2*e1/((psfc(i)/1000.)-e1)
 237       qsfc_air(i) = qsfc(i) * mavail(i)
 238       thetav(i) = (tsk(i)/pi1d(i)) * (1. + 0.608*qsfc_air(i))  ! WA Assumes ql(env)=0, what if it isn't?
 239       ! WA TEST (R5) set z0t = z0
 240       ! z0t(i) = znt(i) / 10.0   ! WA this is hard coded in Matlab version
 241       z0t(i) = znt(i)
 242
 243       ! Get height and delta at turbulence levels and mass levels
 244       zt(i) = (z(i) - zsrf(i) - znt(i)) / 2.
 245       zm(i) = z(i) - zsrf(i)
 246
 247       ! Gradients at first level
 248       dthdz(i) = (th1d(i)-(tsk(i)/pi1d(i))) / (zt(i) * log10(zm(i)/z0t(i)))
 249       dqtdz(i) = (qv1d(i)-qsfc_air(i)) / (zt(i) * log10(zm(i)/z0t(i)))
 250       dudz(i) = u1d(i) / (zt(i) * log10(zm(i)/znt(i)))
 251       dvdz(i) = v1d(i) / (zt(i) * log10(zm(i)/znt(i)))
 252
 253       ! WA doing this because te_temf may not be initialized,
 254       ! would be better to do it in initialization routine but it's
 255       ! not available in module_physics_init.
 256       if (te_temfx(i) < TEmin) te_temfx(i) = TEmin
 257
 258       if ( hfx(i) > 0.) then
 259          wstr(i) = (g * hd_temfx(i) / thetav(i) * (hfx(i)/(rho(i)*cp))) ** (1./3.)
 260       else
 261          wstr(i) = 0.
 262       end if
 263
 264       ! Find stability parameters and length scale
 265       ! WA Calculation of N should really use d(thetaV)/dz not dthdz
 266       ! WA 7/1/09 allow N to be negative
 267       ! if ( dthdz(i) >= 0.) then
 268       ! N(i) = csqrt(g / thetav(i) * dthdz(i))
 269       ! else
 270       !    N(i) = 0.
 271       ! end if
 272       N2(i) = g / thetav(i) * dthdz(i)
 273       S(i) = sqrt(dudz(i)**2. + dvdz(i)**2.)
 274       ! Ri(i) = N(i)**2. / S(i)**2.
 275       Ri(i) = N2(i) / S(i)**2.
 276       ! if (S(i) < 1e-15) Ri(i) = 1./1e-15
 277       if (S(i) < 1e-15) then
 278          print *,'In TEMF SFC Limiting Ri,S,N2,Ri,u,v = ',S(i),N2(i),Ri(i),u1d(i),v1d(i)
 279          if (N2(i) >= 0) then
 280             Ri(i) = 0.2
 281          else
 282             Ri(i) = -1.
 283          end if
 284       end if
 285       if (Ri(i) > 0.2) then  ! WA TEST to prevent runaway
 286          Ri(i) = 0.2
 287       end if
 288       beta(i) = g / thetav(i)
 289       ! WA 7/1/09 adjust ratio, ftau, fth for Ri>0
 290       if (Ri(i) > 0) then
 291          ratio(i) = Ri(i)/(Cphi**2.*ftau0**2./(2.*Ceps**2.*fth0**2.)+3.*Ri(i))
 292          ftau(i) = ftau0 * ((3./4.) / (1.+4.*Ri(i)) + 1./4.)
 293          fth(i) = fth0 / (1.+4.*Ri(i))
 294          ! TE2(i) = 2. * te_temfx(i) * ratio(i) * N(i)**2. / beta(i)**2.
 295          TE2(i) = 2. * te_temfx(i) * ratio(i) * N2(i) / beta(i)**2.
 296       else
 297          ratio(i) = Ri(i)/(Cphi**2.*ftau0**2./(-2.*Ceps**2.*fth0**2.)+2.*Ri(i))
 298          ftau(i) = ftau0
 299          fth(i) = fth0
 300          TE2(i) = 0.
 301       end if
 302       TKE(i) = te_temfx(i) * (1. - ratio(i))
 303       ustrtilde(i) = sqrt(ftau(i) * TKE(i))
 304       ! linv(i) = 1./karman / zt(i) + abs(fCor(i)) / (Cf*ustrtilde(i)) + N(i)/(CN*ustrtilde(i))
 305       if (N2(i) > 0.) then
 306          linv(i) = 1./karman / zt(i) + abs(fCor(i)) / (Cf*ustrtilde(i)) + sqrt(N2(i))/(CN*ustrtilde(i))
 307       else
 308          linv(i) = 1./karman / zt(i) + abs(fCor(i)) / (Cf*ustrtilde(i))
 309       end if
 310       leps(i) = 1./linv(i)
 311       ! WA TEST (R4) remove lower limit on leps
 312       ! lepsmin(i) = min(0.4*zt(i), 5.)
 313       lepsmin(i) = 0.
 314       leps(i) = max(leps(i),lepsmin(i))
 315
 316
 317       ! Find diffusion coefficients
 318       ! First use basic formulae for stable and neutral cases,
 319       ! then for convective conditions, and finally choose the larger
 320       km(i) = TKE(i)**1.5 * ftau(i)**2. / (-beta(i) * fth(i) * sqrt(TE2(i)) + Ceps * sqrt(TKE(i)*te_temfx(i)) / leps(i))
 321       kh(i) = 2. * leps(i) * fth(i)**2. * TKE(i) / sqrt(te_temfx(i)) / Cphi
 322       km(i) = max(km(i),visc_temf)
 323       kh(i) = max(kh(i),conduc_temf)
 324
 325       ! Surface fluxes
 326       ! WA TEST 11/7/13 use w* as a component of the mean wind inside the
 327       ! u* calculation instead of in the velocity scale below (Felix)
 328       ! ust(i) = sqrt(ftau(i)/ftau0) * sqrt(u1d(i)**2. + v1d(i)**2.) * leps(i) / log(zm(i)/znt(i)) / zt(i)
 329       ust(i) = sqrt(ftau(i)/ftau0) * sqrt(u1d(i)**2. + v1d(i)**2. + (0.5*wstr(i))**2.) * leps(i) / log(zm(i)/znt(i)) / zt(i)
 330       ang(i) = atan2(v1d(i),u1d(i))
 331
 332       ! WA TEST 11/7/13 back to wm = u* but with "whole" wind in u* above
 333       wm(i) = ust(i)
 334       ! Calculate mixed scaling velocity (Moeng & Sullivan 1994 JAS p.1021)
 335       ! wm(i) = 0.5 * (1./5. * (wstr(i)**3. + 5. * ust(i)**3.)) ** (1./3.)
 336
 337       ! WA TEST 2/22/11 average with previous value to reduce instability
 338       wm(i) = (wm(i) + wm_temfx(i)) / 2.0
 339       ! WA TEST 11/26/13 set min value
 340       wm_temfx(i) = max(wm(i),1e-2)
 341
 342       ! Populate surface exchange coefficient variables to go back out
 343       ! for next time step of surface scheme
 344       ! Unit specifications in SLAB and sfclay are conflicting and probably
 345       ! incorrect.  This will give a dynamic heat flux (W/m^2) or moisture
 346       ! flux (kg(water)/(m^2*s)) when multiplied by a difference.
 347       ! These formulae are the same as what's used above to get surface
 348       ! flux from surface temperature and specific humidity.
 349       flhc(i) = rho(i) * cp * fth(i)/fth0 * wm(i) * leps(i) / PrT0 / log(zm(i)/z0t(i)) / zt(i)
 350       flqc(i)  = rho(i) * fth(i)/fth0 * wm(i) * leps(i) / PrT0 / log(zm(i)/z0t(i)) / zt(i) * mavail(i)
 351       exch_temfx(i)  = flqc(i) / mavail(i)
 352       chs(i) = flqc(i) / rho(i) / mavail(i)
 353       ! WA Must exchange coeffs be limited to avoid runaway in some
 354       ! (convective?) conditions?  Something like this is done in sfclay.
 355       ! Doing nothing for now.
 356
 357       ! Populate surface heat and moisture fluxes
 358       hfx(i) = flhc(i) * (tsk(i) - th1d(i)*pi1d(i))
 359       ! qfx(i) = flqc(i) * (qsfc_air(i) - qv1d(i))  ! WA 2/16/11
 360       qfx(i) = flqc(i) * (qsfc(i) - qv1d(i))
 361       qfx(i) = max(qfx(i),0.)  ! WA this is done in sfclay, is it right?
 362       lh(i)=xlv*qfx(i)
 363
 364
 365       ! Populate 10 m winds and 2 m temp and 2 m exchange coeffs
 366       ! WA Note this only works if first mass level is above 10 m
 367       u10(i) = u1d(i) * log(10.0/znt(i)) / log(zm(i)/znt(i))
 368       v10(i) = v1d(i) * log(10.0/znt(i)) / log(zm(i)/znt(i))
 369       t2(i) = (tsk(i)/pi1d(i) + (th1d(i) - tsk(i)/pi1d(i)) * log(2.0/z0t(i)) / log(zm(i)/z0t(i))) * pi1d(i)  ! WA this should also use pi at z0
 370       th2(i) = t2(i) / pi1d(i)
 371       q2(i) = (qsfc_air(i) + (qv1d(i) - qsfc_air(i)) * log(2.0/znt(i)) / log(zm(i)/znt(i)))
 372       ! WA are these correct?  Difference between chs2 and cqs2 is unclear
 373       ! At the moment the only difference is z0t vs. znt
 374       chs2(i)  = fth(i)/fth0 * wm(i) * leps(i) / PrT0 / log(2.0/z0t(i)) / zt(i)
 375       cqs2(i)  = fth(i)/fth0 * wm(i) * leps(i) / PrT0 / log(2.0/znt(i)) / zt(i)
 376
 377       ! Calculate qgh (saturated at first-level temp) and cpm
 378       e1=svp1*exp(svp2*((th1d(i)*pi1d(i))-svpt0)/((th1d(i)*pi1d(i))-svp3))
 379       qgh(i)=ep2*e1/((p1d(i)/1000.)-e1)
 380       cpm(i)=cp*(1.+0.8*qv1d(i))
 381
 382    end do  ! Main loop
 383
 384    END SUBROUTINE temfsfclay1d
 385
 386 !====================================================================
 387    SUBROUTINE temfsfclayinit( restart, allowed_to_read,                &
 388                               wm_temf,                                 &
 389                               ids, ide, jds, jde, kds, kde,            &
 390                               ims, ime, jms, jme, kms, kme,            &
 391                               its, ite, jts, jte, kts, kte                 )
 392
 393    logical , intent(in)          :: restart, allowed_to_read
 394    REAL,     DIMENSION( ims:ime, jms:jme )           , &
 395                 INTENT(  OUT) :: wm_temf
 396    integer , intent(in)          ::  ids, ide, jds, jde, kds, kde,             &
 397                                      ims, ime, jms, jme, kms, kme,             &
 398                                      its, ite, jts, jte, kts, kte
 399
 400 ! Local variables
 401    integer :: i, j, itf, jtf
 402 !
 403    CALL wrf_debug( 100, 'in temfsfclayinit' )
 404    jtf = min0(jte,jde-1)
 405    itf = min0(ite,ide-1)
 406 !
 407    if(.not.restart)then
 408      do j = jts,jtf
 409      do i = its,itf
 410      ! do j = jms,jme
 411      ! do i = ims,ime
 412         wm_temf(i,j) = 0.0
 413      enddo
 414      enddo
 415    endif
 416
 417    END SUBROUTINE temfsfclayinit
 418
 419 !-------------------------------------------------------------------
 420
 421 END MODULE module_sf_temfsfclay