PageRenderTime 54ms CodeModel.GetById 16ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_sf_mynn.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1595 lines | 759 code | 294 blank | 542 comment | 5 complexity | 517833b1cd11f840640a0c9b4bc8f439 MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. MODULE module_sf_mynn
    2. 

  2. 

  3. !-------------------------------------------------------------------
    4. 

  4. !Modifications implemented by Joseph Olson NOAA/GSD/AMB - CU/CIRES
    5. 

  5. !for WRFv3.4:
    6. 

  6. !
    7. 

  7. !   BOTH LAND AND WATER:
    8. 

  8. !1) Calculation of stability parameter (z/L) taken from Li et al. (2010 BLM)
    9. 

  9. !2) Fixed isfflx=0 option to turn off scalar fluxes, but keep momentum
    10. 

  10. !   fluxes for idealized studies (credit: Anna Fitch).
    11. 

  11. !3) Kinematic viscosity now varies with temperature
    12. 

  12. !4) Uses Monin-Obukhov flux-profile relationships more consistent with
    13. 

  13. !   those used in the MYNN PBL code.
    14. 

  14. !5) Allows negative QFX, similar to MYJ scheme
    15. 

  15. !
    16. 

  16. !   LAND only:
    17. 

  17. !1) Scalar roughness length parameterization of Yang et al. (2002 QJRMS,
    18. 

  18. !   2008 JAMC) and Chen et al. (2010 J. of Hydromet), with modifications.
    19. 

  19. !2) Relaxed u* minimum from 0.1 to 0.01
    20. 

  20. !
    21. 

  21. !   WATER only:
    22. 

  22. !1) Charnock relationship with varying Charnock parameter is taken from
    23. 

  23. !   the COARE3.0 bulk algorithm (Fairall et al 2003). A formulation for
    24. 

  24. !   the aerodynamic roughness length was also taken from Davis et al (2008)
    25. 

  25. !   and Donelan et al. (2004), but this formulation seems to result in a 
    26. 

  26. !   hign surface wind speed bias for low/moderate wind speeds, so the
    27. 

  27. !   varying Charnock relationship is default.
    28. 

  28. !2) Scalar roughness length is taken from the COARE3.0 bulk algorithm 
    29. 

  29. !   (Fairall et al 2003). A formulation for the scalar roughness length 
    30. 

  30. !   was also taken from Garrat (1992), but since this formula makes zt 
    31. 

  31. !   proportional to zo, heat and moisture fluxes are dependent
    32. 

  32. !   upon which formulation for aerodynamic roughness length is used. 
    33. 

  33. !   Therefore, the COARE3.0 formulation (with slight modifications) is 
    34. 

  34. !   used as default instead.
    35. 

  35. !-------------------------------------------------------------------
    36. 

  36.   USE module_model_constants, only: &
    37. 

  37.        &p1000mb, cp, xlv, ep_2
    38. 

  38. 

  39.   USE module_sf_sfclay, ONLY: sfclayinit
    40. 

  40.   USE module_bl_mynn,   only: tv0, mym_condensation
    41. 

  41.   
    42. 

  42. !-------------------------------------------------------------------
    43. 

  43.   IMPLICIT NONE
    44. 

  44. !-------------------------------------------------------------------
    45. 

  45. 

  46.   REAL, PARAMETER :: xlvcp=xlv/cp, ep_3=1.-ep_2
    47. 

  47.  
    48. 

  48.   REAL, PARAMETER :: wmin=0.1    ! Minimum wind speed
    49. 

  49.   REAL, PARAMETER :: zm2h=7.4    ! = z_0m/z_0h
    50. 

  50. 

  51.   REAL, PARAMETER :: charnock=0.016, bvisc=1.5e-5, z0hsea=5.e-5
    52. 

  52. 

  53.   REAL, PARAMETER :: VCONVC=1.0
    54. 

  54.   REAL, PARAMETER :: CZ0=charnock
    55. 

  55.   
    56. 

  56.   REAL, DIMENSION(0:1000 ),SAVE          :: PSIMTB,PSIHTB
    57. 

  57. 

  58. 

  59. CONTAINS
    60. 

  60. 

  61. !-------------------------------------------------------------------
    62. 

  62.   SUBROUTINE mynn_sf_init_driver(allowed_to_read)
    63. 

  63. 

  64.     LOGICAL, INTENT(in) :: allowed_to_read
    65. 

  65. 

  66.     !Fill the PSIM and PSIH tables. The subroutine "sfclayinit" 
    67. 

  67.     !can be found in module_sf_sfclay.F. This subroutine returns
    68. 

  68.     !the forms from Dyer and Hicks (1974) or Paulson (1970)???.
    69. 

  69.        
    70. 

  70.     CALL sfclayinit(allowed_to_read)
    71. 

  71. 

  72.   END SUBROUTINE mynn_sf_init_driver
    73. 

  73. 

  74. 

  75. !-------------------------------------------------------------------
    76. 

  76.    SUBROUTINE SFCLAY_mynn(U3D,V3D,T3D,QV3D,P3D,dz8w,               &
    77. 

  77.                      CP,G,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM,      &
    78. 

  78.                      ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
    79. 

  79.                      XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, &
    80. 

  80.                      U10,V10,TH2,T2,Q2,                            &
    81. 

  81.                      GZ1OZ0,WSPD,BR,ISFFLX,DX,                     &
    82. 

  82.                      SVP1,SVP2,SVP3,SVPT0,EP1,EP2,                 &
    83. 

  83.                      KARMAN,EOMEG,STBOLT,                          &
    84. 

  84.                      itimestep,ch,th3d,pi3d,qc3d,                  &
    85. 

  85.                      tsq,qsq,cov,qcg,                              &
    86. 

  86. !JOE-add zo/zt
    87. 

  87. !                     z0zt_ratio,BulkRi,wstar,&
    88. 

  88. !JOE-end z0/zt
    89. 

  89.                      ids,ide, jds,jde, kds,kde,                    &
    90. 

  90.                      ims,ime, jms,jme, kms,kme,                    &
    91. 

  91.                      its,ite, jts,jte, kts,kte                     )
    92. 

  92. !-------------------------------------------------------------------
    93. 

  93.       IMPLICIT NONE
    94. 

  94. !-------------------------------------------------------------------
    95. 

  95. !-- U3D         3D u-velocity interpolated to theta points (m/s)
    96. 

  96. !-- V3D         3D v-velocity interpolated to theta points (m/s)
    97. 

  97. !-- T3D         temperature (K)
    98. 

  98. !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
    99. 

  99. !-- P3D         3D pressure (Pa)
    100. 

  100. !-- dz8w        dz between full levels (m)
    101. 

  101. !-- CP          heat capacity at constant pressure for dry air (J/kg/K)
    102. 

  102. !-- G           acceleration due to gravity (m/s^2)
    103. 

  103. !-- ROVCP       R/CP
    104. 

  104. !-- R           gas constant for dry air (J/kg/K)
    105. 

  105. !-- XLV         latent heat of vaporization for water (J/kg)
    106. 

  106. !-- PSFC        surface pressure (Pa)
    107. 

  107. !-- ZNT         roughness length (m)
    108. 

  108. !-- UST         u* in similarity theory (m/s)
    109. 

  109. !-- PBLH        PBL height from previous time (m)
    110. 

  110. !-- MAVAIL      surface moisture availability (between 0 and 1)
    111. 

  111. !-- ZOL         z/L height over Monin-Obukhov length
    112. 

  112. !-- MOL         T* (similarity theory) (K)
    113. 

  113. !-- REGIME      flag indicating PBL regime (stable, unstable, etc.)
    114. 

  114. !-- PSIM        similarity stability function for momentum
    115. 

  115. !-- PSIH        similarity stability function for heat
    116. 

  116. !-- XLAND       land mask (1 for land, 2 for water)
    117. 

  117. !-- HFX         upward heat flux at the surface (W/m^2)
    118. 

  118. !-- QFX         upward moisture flux at the surface (kg/m^2/s)
    119. 

  119. !-- LH          net upward latent heat flux at surface (W/m^2)
    120. 

  120. !-- TSK         surface temperature (K)
    121. 

  121. !-- FLHC        exchange coefficient for heat (W/m^2/K)
    122. 

  122. !-- FLQC        exchange coefficient for moisture (kg/m^2/s)
    123. 

  123. !-- CHS         heat/moisture exchange coefficient for LSM (m/s)
    124. 

  124. !-- QGH         lowest-level saturated mixing ratio
    125. 

  125. !-- U10         diagnostic 10m u wind
    126. 

  126. !-- V10         diagnostic 10m v wind
    127. 

  127. !-- TH2         diagnostic 2m theta (K)
    128. 

  128. !-- T2          diagnostic 2m temperature (K)
    129. 

  129. !-- Q2          diagnostic 2m mixing ratio (kg/kg)
    130. 

  130. !-- GZ1OZ0      log(z/z0) where z0 is roughness length
    131. 

  131. !-- WSPD        wind speed at lowest model level (m/s)
    132. 

  132. !-- BR          bulk Richardson number in surface layer
    133. 

  133. !-- ISFFLX      isfflx=1 for surface heat and moisture fluxes
    134. 

  134. !-- DX          horizontal grid size (m)
    135. 

  135. !-- SVP1        constant for saturation vapor pressure (kPa)
    136. 

  136. !-- SVP2        constant for saturation vapor pressure (dimensionless)
    137. 

  137. !-- SVP3        constant for saturation vapor pressure (K)
    138. 

  138. !-- SVPT0       constant for saturation vapor pressure (K)
    139. 

  139. !-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
    140. 

  140. !-- EP2         constant for specific humidity calculation 
    141. 

  141. !               (R_d/R_v) (dimensionless)
    142. 

  142. !-- KARMAN      Von Karman constant
    143. 

  143. !-- EOMEG       angular velocity of earth's rotation (rad/s)
    144. 

  144. !-- STBOLT      Stefan-Boltzmann constant (W/m^2/K^4)
    145. 

  145. !-- ids         start index for i in domain
    146. 

  146. !-- ide         end index for i in domain
    147. 

  147. !-- jds         start index for j in domain
    148. 

  148. !-- jde         end index for j in domain
    149. 

  149. !-- kds         start index for k in domain
    150. 

  150. !-- kde         end index for k in domain
    151. 

  151. !-- ims         start index for i in memory
    152. 

  152. !-- ime         end index for i in memory
    153. 

  153. !-- jms         start index for j in memory
    154. 

  154. !-- jme         end index for j in memory
    155. 

  155. !-- kms         start index for k in memory
    156. 

  156. !-- kme         end index for k in memory
    157. 

  157. !-- its         start index for i in tile
    158. 

  158. !-- ite         end index for i in tile
    159. 

  159. !-- jts         start index for j in tile
    160. 

  160. !-- jte         end index for j in tile
    161. 

  161. !-- kts         start index for k in tile
    162. 

  162. !-- kte         end index for k in tile
    163. 

  163. !-------------------------------------------------------------------
    164. 

  164.       INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
    165. 

  165.                                         ims,ime, jms,jme, kms,kme, &
    166. 

  166.                                         its,ite, jts,jte, kts,kte
    167. 

  167. !                                                               
    168. 

  168.       INTEGER,  INTENT(IN )   ::        ISFFLX
    169. 

  169.       REAL,     INTENT(IN )   ::        SVP1,SVP2,SVP3,SVPT0
    170. 

  170.       REAL,     INTENT(IN )   ::        EP1,EP2,KARMAN,EOMEG,STBOLT
    171. 

  171. !
    172. 

  172.       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
    173. 

  173.                 INTENT(IN   )   ::                           dz8w
    174. 

  174.                                         
    175. 

  175.       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
    176. 

  176.                 INTENT(IN   )   ::                           QV3D, &
    177. 

  177.                                                               P3D, &
    178. 

  178.                                                               T3D, &
    179. 

  179.                                                              QC3D,&
    180. 

  180.                                             th3d,pi3d,tsq,qsq,cov
    181. 

  181. 

  182.       INTEGER, INTENT(in) :: itimestep
    183. 

  183. 

  184.       REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(IN) ::&
    185. 

  185.            &    qcg
    186. 

  186.       
    187. 

  187.       REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) ::&
    188. 

  188.            & ch
    189. 

  189. 

  190.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    191. 

  191.                 INTENT(IN   )               ::             MAVAIL, &
    192. 

  192.                                                              PBLH, &
    193. 

  193.                                                             XLAND, &
    194. 

  194.                                                               TSK
    195. 

  195.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    196. 

  196.                 INTENT(OUT  )               ::                U10, &
    197. 

  197.                                                               V10, &
    198. 

  198.                                                               TH2, &
    199. 

  199.                                                                T2, &
    200. 

  200. !JOE-use value from LSM                                        Q2, &
    201. 

  201.                                                                Q2
    202. 

  202. !JOE-moved down below                                          QSFC
    203. 

  203. 

  204. !
    205. 

  205.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    206. 

  206.                 INTENT(INOUT)               ::             REGIME, &
    207. 

  207.                                                               HFX, &
    208. 

  208.                                                               QFX, &
    209. 

  209.                                                                LH, &
    210. 

  210.                                                     MOL,RMOL,QSFC
    211. 

  211. !m the following 5 are change to memory size
    212. 

  212. !
    213. 

  213.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    214. 

  214.                 INTENT(INOUT)   ::                 GZ1OZ0,WSPD,BR, &
    215. 

  215.                                                         PSIM,PSIH
    216. 

  216. 

  217.       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
    218. 

  218.                 INTENT(IN   )   ::                            U3D, &
    219. 

  219.                                                               V3D
    220. 

  220.                                         
    221. 

  221.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    222. 

  222.                 INTENT(IN   )               ::               PSFC
    223. 

  223. 

  224.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    225. 

  225.                 INTENT(INOUT)   ::                            ZNT, &
    226. 

  226.                                                               ZOL, &
    227. 

  227.                                                               UST, &
    228. 

  228.                                                               CPM, &
    229. 

  229.                                                              CHS2, &
    230. 

  230.                                                              CQS2, &
    231. 

  231.                                                               CHS
    232. 

  232. 

  233.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    234. 

  234.                 INTENT(INOUT)   ::                      FLHC,FLQC
    235. 

  235. 

  236.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    237. 

  237.                 INTENT(INOUT)   ::                            QGH
    238. 

  238. 

  239. !JOE-begin
    240. 

  240. !      REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    241. 

  241. !                INTENT(OUT)               ::         z0zt_ratio, &
    242. 

  242. !                                                      BulkRi,wstar
    243. 

  243. !JOE-end 
    244. 

  244.                                     
    245. 

  245.       REAL,     INTENT(IN   )               ::   CP,G,ROVCP,R,XLV,DX
    246. 

  246. 

  247. 

  248. !----------- LOCAL VARS -----------------------------------
    249. 

  249. 

  250.       REAL,     DIMENSION( its:ite ) ::                       U1D, &
    251. 

  251.                                                               V1D, &
    252. 

  252.                                                              QV1D, &
    253. 

  253.                                                               P1D, &
    254. 

  254.                                                          T1D,qc1d
    255. 

  255. 

  256.       REAL,     DIMENSION( its:ite ) ::                    dz8w1d
    257. 

  257. 

  258.       REAL,     DIMENSION( its:ite ) ::  vt1,vq1
    259. 

  259.       REAL,     DIMENSION(kts:kts+1) ::  thl, qw, vt, vq
    260. 

  260.       REAL                           ::  ql
    261. 

  261. 

  262.       INTEGER ::  I,J,K
    263. 

  263. !-----------------------------------------------------------
    264. 

  264. 

  265.       DO J=jts,jte
    266. 

  266.         DO i=its,ite
    267. 

  267.           dz8w1d(I) = dz8w(i,kts,j)
    268. 

  268.         ENDDO
    269. 

  269.    
    270. 

  270.         DO i=its,ite
    271. 

  271.            U1D(i) =U3D(i,kts,j)
    272. 

  272.            V1D(i) =V3D(i,kts,j)
    273. 

  273.            QV1D(i)=QV3D(i,kts,j)
    274. 

  274.            QC1D(i)=QC3D(i,kts,j)
    275. 

  275.            P1D(i) =P3D(i,kts,j)
    276. 

  276.            T1D(i) =T3D(i,kts,j)
    277. 

  277.         ENDDO
    278. 

  278. 

  279.         IF (itimestep==1) THEN
    280. 

  280.            DO i=its,ite
    281. 

  281.               vt1(i)=0.
    282. 

  282.               vq1(i)=0.
    283. 

  283.               UST(i,j)=0.02*SQRT(U1D(i)*U1D(i) + V1D(i)*V1D(i))
    284. 

  284.               MOL(i,j)=0.     ! Tstar
    285. 

  285.            ENDDO
    286. 

  286.         ELSE
    287. 

  287.            DO i=its,ite
    288. 

  288.               do k = kts,kts+1
    289. 

  289.                 ql = qc3d(i,k,j)/(1.+qc3d(i,k,j))
    290. 

  290.                 qw(k) = qv3d(i,k,j)/(1.+qv3d(i,k,j)) + ql
    291. 

  291.                 thl(k) = th3d(i,k,j)-xlvcp*ql/pi3d(i,k,j)
    292. 

  292.               end do
    293. 

  293. 

  294. ! NOTE: The last grid number is kts+1 instead of kte.
    295. 

  295.               CALL mym_condensation (kts,kts+1, &
    296. 

  296.                    &            dz8w(i,kts:kts+1,j), &
    297. 

  297.                    &            thl(kts:kts+1), qw(kts:kts+1), &
    298. 

  298.                    &            p3d(i,kts:kts+1,j),&
    299. 

  299.                    &            pi3d(i,kts:kts+1,j), &
    300. 

  300.                    &            tsq(i,kts:kts+1,j), &
    301. 

  301.                    &            qsq(i,kts:kts+1,j), &
    302. 

  302.                    &            cov(i,kts:kts+1,j), &
    303. 

  303.                    &            vt(kts:kts+1), vq(kts:kts+1))
    304. 

  304. 

  305.               vt1(i) = vt(kts)
    306. 

  306.               vq1(i) = vq(kts)
    307. 

  307.            ENDDO
    308. 

  308.         ENDIF
    309. 

  309. 

  310.         CALL SFCLAY1D_mynn(J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,          &
    311. 

  311.                 CP,G,ROVCP,R,XLV,PSFC(ims,j),CHS(ims,j),CHS2(ims,j),&
    312. 

  312.                 CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j),   &
    313. 

  313.                 ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j),    &
    314. 

  314.                 MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j),  &
    315. 

  315.                 XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j),     &
    316. 

  316.                 U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j),        &
    317. 

  317.                 Q2(ims,j),FLHC(ims,j),FLQC(ims,j),QGH(ims,j),      &
    318. 

  318.                 QSFC(ims,j),LH(ims,j),                             &
    319. 

  319.                 GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX,     &
    320. 

  320.                 SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,  &
    321. 

  321.                 ch(ims,j),vt1,vq1,qc1d,qcg(ims,j),&
    322. 

  322. !JOE-begin
    323. 

  323. !                z0zt_ratio(ims,j),BulkRi(ims,j),wstar(ims,j),      &
    324. 

  324.                 itimestep,&
    325. 

  325. !JOE-end
    326. 

  326.                 ids,ide, jds,jde, kds,kde,                         &
    327. 

  327.                 ims,ime, jms,jme, kms,kme,                         &
    328. 

  328.                 its,ite, jts,jte, kts,kte                          )
    329. 

  329. 

  330.       ENDDO
    331. 

  331. 

  332. 

  333.     END SUBROUTINE SFCLAY_MYNN
    334. 

  334. 

  335. !-------------------------------------------------------------------
    336. 

  336.    SUBROUTINE SFCLAY1D_mynn(J,UX,VX,T1D,QV1D,P1D,dz8w1d,                &
    337. 

  337.                      CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, &
    338. 

  338.                      ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH,      &
    339. 

  339.                      XLAND,HFX,QFX,TSK,                            &
    340. 

  340.                      U10,V10,TH2,T2,Q2,FLHC,FLQC,QGH,              &
    341. 

  341.                      QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX,             &
    342. 

  342.                      SVP1,SVP2,SVP3,SVPT0,EP1,EP2,                 &
    343. 

  343.                      KARMAN,EOMEG,STBOLT,                          &
    344. 

  344.                      ch,vt1,vq1,qc1d,qcg,&
    345. 

  345. !JOE-add
    346. 

  346. !                     zratio,BRi,wstar,itimestep,                   &
    347. 

  347.                      itimestep,                                    &
    348. 

  348. !JOE-end
    349. 

  349.                      ids,ide, jds,jde, kds,kde,                    &
    350. 

  350.                      ims,ime, jms,jme, kms,kme,                    &
    351. 

  351.                      its,ite, jts,jte, kts,kte                     )
    352. 

  352. !-------------------------------------------------------------------
    353. 

  353.       IMPLICIT NONE
    354. 

  354. !-------------------------------------------------------------------
    355. 

  355.       REAL,     PARAMETER  :: XKA=2.4E-5   !molecular diffusivity
    356. 

  356.       REAL,     PARAMETER  :: PRT=1.       !prandlt number
    357. 

  357. 

  358.       INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
    359. 

  359.                                         ims,ime, jms,jme, kms,kme, &
    360. 

  360.                                         its,ite, jts,jte, kts,kte, &
    361. 

  361.                                         J
    362. 

  362. !                                                               
    363. 

  363.       INTEGER,  INTENT(in)    :: itimestep
    364. 

  364.       INTEGER,  INTENT(IN )   ::        ISFFLX
    365. 

  365.       REAL,     INTENT(IN )   ::        SVP1,SVP2,SVP3,SVPT0
    366. 

  366.       REAL,     INTENT(IN )   ::        EP1,EP2,KARMAN,EOMEG,STBOLT
    367. 

  367. 

  368. !
    369. 

  369.       REAL,     DIMENSION( ims:ime )                             , &
    370. 

  370.                 INTENT(IN   )               ::             MAVAIL, &
    371. 

  371.                                                              PBLH, &
    372. 

  372.                                                             XLAND, &
    373. 

  373.                                                               TSK
    374. 

  374. !
    375. 

  375.       REAL,     DIMENSION( ims:ime )                             , &
    376. 

  376.                 INTENT(IN   )               ::             PSFCPA
    377. 

  377. 

  378.       REAL,     DIMENSION( ims:ime )                             , &
    379. 

  379.                 INTENT(INOUT)               ::             REGIME, &
    380. 

  380.                                                               HFX, &
    381. 

  381.                                                               QFX, &
    382. 

  382.                                                          MOL,RMOL
    383. 

  383. !m the following 5 are changed to memory size---
    384. 

  384. !
    385. 

  385.       REAL,     DIMENSION( ims:ime )                             , &
    386. 

  386.                 INTENT(INOUT)   ::                 GZ1OZ0,WSPD,BR, &
    387. 

  387.                                                         PSIM,PSIH
    388. 

  388. 

  389.       REAL,     DIMENSION( ims:ime )                             , &
    390. 

  390.                 INTENT(INOUT)   ::                            ZNT, &
    391. 

  391.                                                               ZOL, &
    392. 

  392.                                                               UST, &
    393. 

  393.                                                               CPM, &
    394. 

  394.                                                              CHS2, &
    395. 

  395.                                                              CQS2, &
    396. 

  396.                                                               CHS
    397. 

  397. !JOE-add
    398. 

  398. !      REAL,     DIMENSION( ims:ime )                             , &
    399. 

  399. !                INTENT(OUT)     ::                zratio,BRi,wstar
    400. 

  400.       REAL,     DIMENSION( ims:ime )  ::           zratio,BRi,wstar
    401. 

  401. !JOE-end
    402. 

  402. 

  403.       REAL,     DIMENSION( ims:ime )                             , &
    404. 

  404.                 INTENT(INOUT)   ::                      FLHC,FLQC
    405. 

  405. 

  406.       REAL,     DIMENSION( ims:ime )                             , &
    407. 

  407.                 INTENT(INOUT)   ::                       QGH,QSFC
    408. 

  408. 

  409.       REAL,     DIMENSION( ims:ime )                             , &
    410. 

  410.                 INTENT(OUT)     ::                        U10,V10, &
    411. 

  411. !JOE-make qsfc inout (moved up)                   TH2,T2,Q2,QSFC,LH
    412. 

  412.                                                 TH2,T2,Q2,LH
    413. 

  413.                                     
    414. 

  414.       REAL,     INTENT(IN)               ::   CP,G,ROVCP,R,XLV,DX
    415. 

  415. 

  416. ! MODULE-LOCAL VARIABLES, DEFINED IN SUBROUTINE SFCLAY
    417. 

  417.       REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::  dz8w1d
    418. 

  418. 

  419.       REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::      UX, &
    420. 

  420.                                                                VX, &
    421. 

  421.                                                              QV1D, &
    422. 

  422.                                                               P1D, &
    423. 

  423.                                                          T1D,qc1d
    424. 

  424.  
    425. 

  425.       REAL,     DIMENSION( ims:ime ), INTENT(IN)    :: qcg
    426. 

  426.       REAL,     DIMENSION( ims:ime ), INTENT(INOUT) :: ch
    427. 

  427. 

  428.       REAL,     DIMENSION( its:ite ), INTENT(IN)    :: vt1,vq1
    429. 

  429. 

  430. ! LOCAL VARS
    431. 

  431. 

  432.       REAL,     DIMENSION( its:ite ) :: z_t,z_q
    433. 

  433. 

  434.       REAL :: thl1,sqv1,sqc1,exner1,sqvg,sqcg,vv,ww
    435. 

  435. 

  436.       REAL,     DIMENSION( its:ite )        ::                 ZA, &
    437. 

  437.                                                         THVX,ZQKL, &
    438. 

  438.                                                            THX,QX, &
    439. 

  439.                                                             PSIH2, &
    440. 

  440.                                                             PSIM2, &
    441. 

  441.                                                            PSIH10, &
    442. 

  442.                                                            PSIM10, &
    443. 

  443.                                                            GZ2OZ0, &
    444. 

  444.                                                           GZ10OZ0
    445. 

  445. !
    446. 

  446.       REAL,     DIMENSION( its:ite )        ::        RHOX,GOVRTH
    447. 

  447. !
    448. 

  448.       REAL,     DIMENSION( its:ite)         ::          SCR4
    449. 

  449.       REAL,     DIMENSION( its:ite )        ::         THGB, PSFC
    450. 

  450. 

  451.       REAL,     DIMENSION( its:ite )        :: GZ2OZt,GZ10OZt,GZ1OZt
    452. 

  452. 

  453. !
    454. 

  454.       INTEGER ::  N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10
    455. 

  455. 

  456.       REAL    ::  PL,THCON,TVCON,E1
    457. 

  457.       REAL    ::  ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10
    458. 

  458.       REAL    ::  DTG,PSIX,USTM,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2
    459. 

  459.       REAL    ::  FLUXC,VSGD
    460. 

  460.       real    ::  restar,VISC,psilim
    461. 

  461. !-------------------------------------------------------------------
    462. 

  462. 

  463. !----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:  
    464. 

  464.       DO I=its,ite
    465. 

  465.          ! PSFC cmb
    466. 

  466.          PSFC(I)=PSFCPA(I)/1000.
    467. 

  467.          THGB(I)=TSK(I)*(100./PSFC(I))**ROVCP                
    468. 

  468.       ENDDO                                              
    469. 

  469. !                                                                                   
    470. 

  470. !     SCR4(I,K) STORES EITHER TEMPERATURE OR VIRTUAL TEMPERATURE,
    471. 

  471. !     DEPENDING ON IFDRY (CURRENTLY NOT USED, SO SCR4 == TVX).                                       
    472. 

  472.                                                                                  
    473. 

  473.       DO 30 I=its,ite
    474. 

  474.          ! PL cmb
    475. 

  475.          PL=P1D(I)/1000.                                                   
    476. 

  476.          THCON=(100./PL)**ROVCP                                                 
    477. 

  477.          THX(I)=T1D(I)*THCON                                               
    478. 

  478.          SCR4(I)=T1D(I)                                                    
    479. 

  479.          THVX(I)=THX(I)                                                     
    480. 

  480.          QX(I)=0.                                                             
    481. 

  481.    30 CONTINUE                                                                 
    482. 

  482. 

  483.       ! INITIALIZE SOME VARIABLES HERE:
    484. 

  484.       DO I=its,ite
    485. 

  485.          QGH(I)=0.                                                                
    486. 

  486.          FLHC(I)=0.                                                               
    487. 

  487.          FLQC(I)=0.  
    488. 

  488.          CPM(I)=CP                                                             
    489. 

  489.       ENDDO
    490. 

  490.                                                                                 
    491. 

  491. !     IF(IDRY.EQ.1)GOTO 80                                                   
    492. 

  492.       DO 50 I=its,ite
    493. 

  493.          QX(I)=QV1D(I)/(1.+QV1D(I))        !CONVERT TO SPEC HUM
    494. 

  494.          TVCON=(1.+EP1*QX(I))                                      
    495. 

  495.          THVX(I)=THX(I)*TVCON                                               
    496. 

  496.          SCR4(I)=T1D(I)*TVCON                                              
    497. 

  497.    50 CONTINUE                                                                 
    498. 

  498. !                                                                                
    499. 

  499.       DO 60 I=its,ite
    500. 

  500.         IF (TSK(I) .LT. 273.15) THEN
    501. 

  501.            !SATURATION VAPOR PRESSURE WRT ICE
    502. 

  502.            E1=SVP1*EXP(4648*(1./273.15 - 1./TSK(I)) - &
    503. 

  503.                   11.64*LOG(273.15/TSK(I)) + 0.02265*(273.15 - TSK(I)))
    504. 

  504.         ELSE
    505. 

  505.            !SATURATION VAPOR PRESSURE WRT WATER
    506. 

  506.            E1=SVP1*EXP(SVP2*(TSK(I)-SVPT0)/(TSK(I)-SVP3))
    507. 

  507.         ENDIF
    508. 

  508.         QSFC(I)=EP2*E1/(PSFC(I)-ep_3*E1)   !specific humidity
    509. 

  509.         !QSFC(I)=EP2*E1/(PSFC(I)-E1)       !mixing ratio                                  
    510. 

  510.         !FOR LAND POINTS, QSFC can come from previous time step (in LSM)
    511. 

  511.         !if(xland(i).gt.1.5 .or. QSFC(i).le.0.0) QSFC(I)=EP2*E1/(PSFC(I)-ep_3*E1)     
    512. 

  512.   
    513. 

  513.         ! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP CONSISTENT WITH MYJSFC CHANGE
    514. 

  514.         ! Q2SAT = QGH IN LSM
    515. 

  515.         E1=SVP1*EXP(SVP2*(T1D(I)-SVPT0)/(T1D(I)-SVP3))
    516. 

  516.         PL=P1D(I)/1000.
    517. 

  517.         QGH(I)=EP2*E1/(PL-ep_3*E1)     !specific humidity
    518. 

  518.         !QGH(I)=EP2*E1/(PL-E1)         !mixing ratio    
    519. 

  519.         CPM(I)=CP*(1.+0.84*QX(I)/(1.-qx(i)))
    520. 

  520. 

  521.    60 CONTINUE                                                                   
    522. 

  522.    80 CONTINUE
    523. 

  523.                                                                                  
    524. 

  524. !-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND             
    525. 

  525. !     LEVEL, AND THE LAYER THICKNESSES.                                          
    526. 

  526.                                                                                  
    527. 

  527.       DO I=its,ite
    528. 

  528.          RHOX(I)=PSFC(I)*1000./(R*SCR4(I))                                       
    529. 

  529.          ZQKL(I)=dz8w1d(I)            !first full-sigma level
    530. 

  530.          ZA(I)=0.5*ZQKL(I)            !first half-sigma level                            
    531. 

  531.          GOVRTH(I)=G/THX(I)                                                    
    532. 

  532.       ENDDO
    533. 

  533.                                                                                  
    534. 

  534. !-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO               
    535. 

  535. !     AKB(1976), EQ(12).                                                         
    536. 

  536.                    
    537. 

  537.       DO 260 I=its,ite
    538. 

  538. 

  539.          WSPD(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))                        
    540. 

  540. 

  541.          IF((XLAND(I)-1.5).GE.0)THEN
    542. 

  542. 

  543.             !COMPUTE KINEMATIC VISCOSITY
    544. 

  544.             VISC=(1.32+0.009*(T1D(I)-273.15))*1.E-5
    545. 

  545. 

  546.             !--------------------------------------
    547. 

  547.             ! WATER
    548. 

  548.             !--------------------------------------
    549. 

  549.             !GET ZNT
    550. 

  550.             !--------------------------------------
    551. 

  551.             !CALL davis_etal_2008(ZNT(i),UST(i))
    552. 

  552.             !CALL Taylor_Yelland_2001(ZNT(i),UST(i),WSPD(i))
    553. 

  553.             CALL charnock_1955(ZNT(i),UST(i),WSPD(i),visc)
    554. 

  554. 

  555.             !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar) USING NEW ZNT
    556. 

  556.             ! AHW: Garratt formula: Calculate roughness Reynolds number
    557. 

  557.             !      Kinematic viscosity of air (linear approx to
    558. 

  558.             !      temp dependence at sea level)
    559. 

  559.             restar=MAX(ust(i)*ZNT(i)/visc, 0.1)
    560. 

  560. 

  561.             !--------------------------------------
    562. 

  562.             !GET z_t and z_q
    563. 

  563.             !--------------------------------------
    564. 

  564.             !CALL zilitinkevich_1995(ZNT(i),z_t(i),z_q(i),restar,UST(I),KARMAN,XLAND(I))
    565. 

  565.             !CALL garrat_1992(z_t(i),z_q(i),ZNT(i),restar,XLAND(I))
    566. 

  566.             CALL fairall_2001(z_t(i),z_q(i),restar,UST(i),visc)
    567. 

  567. 

  568.             zratio(i)=znt(i)/z_t(i)
    569. 

  569. 

  570.          ELSE
    571. 

  571. 

  572.             !--------------------------------------
    573. 

  573.             ! LAND
    574. 

  574.             !--------------------------------------
    575. 

  575. 

  576.             !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar) USING DEFAULT ZNT
    577. 

  577.             !restar=MAX(ust(i)*ZNT(i)/bvisc, 1.0)
    578. 

  578.             VISC=(1.32+0.009*(T1D(I)-273.15))*1.E-5
    579. 

  579.             restar=MAX(ust(i)*ZNT(i)/visc, 0.1)
    580. 

  580. 

  581.             !--------------------------------------
    582. 

  582.             !GET z_t and z_q
    583. 

  583.             !--------------------------------------
    584. 

  584.             !CALL zilitinkevich_1995(ZNT(i),z_t(i),z_q(i),restar,UST(I),KARMAN,XLAND(I))
    585. 

  585.             !CALL garrat_1992(z_t(i),z_q(i),ZNT(i),restar,XLAND(I))
    586. 

  586.             CALL Yang_2008(ZNT(i),z_t(i),z_q(i),UST(i),MOL(I),restar,visc,XLAND(I))
    587. 

  587. 

  588.             zratio(i)=znt(i)/z_t(i)
    589. 

  589. 

  590.          ENDIF
    591. 

  591. 

  592. 

  593.          GZ1OZ0(I)= LOG(ZA(I)/ZNT(I))
    594. 

  594.          GZ1OZt(I)= LOG(ZA(I)/z_t(i))           
    595. 

  595.          GZ2OZ0(I)= LOG(2./ZNT(I))                                        
    596. 

  596.          GZ2OZt(I)= LOG(2./z_t(i))                                        
    597. 

  597.          GZ10OZ0(I)=LOG(10./ZNT(I)) 
    598. 

  598.          GZ10OZt(I)=LOG(10./z_t(i)) 
    599. 

  599. 

  600.          !account for partial condensation
    601. 

  601.          exner1=(p1d(i)/p1000mb)**ROVCP
    602. 

  602.          sqc1=qc1d(i)/(1.+qc1d(i))         !convert to spec hum.
    603. 

  603.          sqv1=qx(i)
    604. 

  604.          thl1=THX(I)-xlvcp/exner1*sqc1
    605. 

  605.          sqvg=qsfc(i)
    606. 

  606.          sqcg=qcg(i)/(1.+qcg(i))  !convert to specific humidity
    607. 

  607. 

  608.          vv = thl1-THGB(I)
    609. 

  609.          ww = mavail(i)*(sqv1-sqvg) + (sqc1-sqcg)
    610. 

  610. 

  611.          TSKV=THGB(I)*(1.+EP1*QSFC(I)*MAVAIL(I)) 
    612. 

  612.          !TSKV=THGB(I)*(1.+EP1*QSFC(I))                     
    613. 

  613. 

  614.          !DTHVDZ=(THVX(I)-TSKV)                     
    615. 

  615.          DTHVDZ= (vt1(i) + 1.0)*vv + (vq1(i) + tv0)*ww
    616. 

  616. 

  617.          !--------------------------------------------------------
    618. 

  618.          !  Calculate the convective velocity scale (WSTAR) and 
    619. 

  619.          !  subgrid-scale velocity (VSGD) following Beljaars (1995, QJRMS) 
    620. 

  620.          !  and Mahrt and Sun (1995, MWR), respectively
    621. 

  621.          !-------------------------------------------------------
    622. 

  622.          !       VCONV = 0.25*sqrt(g/tskv*pblh(i)*dthvm)
    623. 

  623.          !  Use Beljaars over land, old MM5 (Wyngaard) formula over water
    624. 

  624.          if (xland(i).lt.1.5) then     !LAND (xland == 1)
    625. 

  625. 

  626.             fluxc = max(hfx(i)/rhox(i)/cp                    &
    627. 

  627.                  + ep1*tskv*qfx(i)/rhox(i),0.)
    628. 

  628.             !JOE:VCONV = vconvc*(g/TSK(i)*pblh(i)*fluxc)**.33
    629. 

  629.             WSTAR(I) = vconvc*(g/TSK(i)*pblh(i)*fluxc)**.33
    630. 

  630. 

  631.          else                          !WATER (xland == 2)
    632. 

  632. 

  633.             IF(-DTHVDZ.GE.0)THEN
    634. 

  634.                DTHVM=-DTHVDZ
    635. 

  635.             ELSE
    636. 

  636.                DTHVM=0.
    637. 

  637.             ENDIF
    638. 

  638.             !JOE-the Wyngaard formula is 2-3 times larger than the Beljaars 
    639. 

  639.             !formula, so switch to Beljaars for water, but use VCONVC = 1.25,
    640. 

  640.             !as in the COARE3.0 bulk parameterizations.
    641. 

  641.             !WSTAR(I) = 2.*SQRT(DTHVM)
    642. 

  642.             fluxc = max(hfx(i)/rhox(i)/cp                    &
    643. 

  643.                 + ep1*tskv*qfx(i)/rhox(i),0.)
    644. 

  644.             WSTAR(I) = 1.25*(g/TSK(i)*pblh(i)*fluxc)**.33
    645. 

  645. 

  646.          endif
    647. 

  647. 

  648.          !--------------------------------------------------------
    649. 

  649.          ! Mahrt and Sun low-res correction
    650. 

  650.          ! (for 13 km ~ 0.37 m/s; for 3 km == 0 m/s)
    651. 

  651.          !--------------------------------------------------------
    652. 

  652.          VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33
    653. 

  653. 

  654.          WSPD(I)=SQRT(WSPD(I)*WSPD(I)+WSTAR(I)*WSTAR(I)+vsgd*vsgd)
    655. 

  655.          WSPD(I)=AMAX1(WSPD(I),wmin)
    656. 

  656. 

  657.          !--------------------------------------------------------
    658. 

  658.          ! CALCULATE THE BULK RICHARDSON NUMBER OF SURFACE LAYER, 
    659. 

  659.          ! ACCORDING TO AKB(1976), EQ(12). 
    660. 

  660.          !--------------------------------------------------------
    661. 

  661.          BR(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD(I)*WSPD(I))
    662. 

  662.          !SET LIMITS ACCORDING TO Li et al. (2010) Boundary-Layer Meteorol (p.158)
    663. 

  663.          BR(I)=MAX(BR(I),-2.0)
    664. 

  664.          BR(I)=MIN(BR(I),1.0)
    665. 

  665.          BRi(I)=BR(I)  !new variable for output - BR is not a "state" variable.
    666. 

  666.                
    667. 

  667.          !  IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2 (STABLE)
    668. 

  668.          if (itimestep .GT. 1) THEN
    669. 

  669.             IF(MOL(I).LT.0.)BR(I)=AMIN1(BR(I),0.0)
    670. 

  670.          ENDIF
    671. 

  671. !jdf
    672. 

  672.          RMOL(I)=-GOVRTH(I)*DTHVDZ*ZA(I)*KARMAN
    673. 

  673. !jdf
    674. 

  674. !JOE-begin
    675. 

  675. !        IF(I .eq. 2)THEN
    676. 

  676. !          write(*,1006)"BR:",BR(I)," fluxc:",fluxc," vt1:",vt1(i)," vq1:",vq1(i)
    677. 

  677. !          write(*,1007)"XLAND:",XLAND(I)," WSPD:",WSPD(I)," DTHVDZ:",DTHVDZ," WSTAR:",WSTAR(I)
    678. 

  678. !        ENDIF
    679. 

  679. !JOE-end
    680. 

  680.                          
    681. 

  681.   260 CONTINUE                                                                   
    682. 

  682. 

  683.  1006   format(A,F7.3,A,f9.4,A,f9.5,A,f9.4)
    684. 

  684.  1007   format(A,F2.0,A,f6.2,A,f7.3,A,f7.2)
    685. 

  685. 

  686. !--------------------------------------------------------------------      
    687. 

  687. !--- DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:            
    688. 

  688. ! 
    689. 

  689. !                                                                                
    690. 

  690. !     THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)           
    691. 

  691. !     AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).                              
    692. 

  692. !                                                                                
    693. 

  693. !     CRITERIA FOR THE CLASSES ARE AS FOLLOWS:                                   
    694. 

  694. !                                                                                
    695. 

  695. !        1. BR .GE. 0.2;                                                         
    696. 

  696. !               REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),               
    697. 

  697. !                                                                                
    698. 

  698. !        2. BR .LT. 0.2 .AND. BR .GT. 0.0;                                       
    699. 

  699. !               REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS                
    700. 

  700. !               (REGIME=2),                                                      
    701. 

  701. !                                                                                
    702. 

  702. !        3. BR .EQ. 0.0                                                          
    703. 

  703. !               REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),              
    704. 

  704. !                                                                                
    705. 

  705. !        4. BR .LT. 0.0                                                          
    706. 

  706. !               REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).                
    707. 

  707. !                                                                                
    708. 

  708. !--------------------------------------------------------------------
    709. 

  709. 

  710.   DO 320 I=its,ite
    711. 

  711. 

  712.      IF (BR(I) .GT. 0.2) THEN                                           
    713. 

  713.         !===================================================
    714. 

  714.         !---CLASS 1; STABLE (NIGHTTIME) CONDITIONS:                                    
    715. 

  715.         !===================================================
    716. 

  716.         REGIME(I)=1.
    717. 

  717. 

  718.         !COMPUTE z/L
    719. 

  719.         CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNT(I),zratio(I))
    720. 

  720. 

  721.         !COMPUTE PSIM and PSIH
    722. 

  722.         IF((XLAND(I)-1.5).GE.0)THEN                                            
    723. 

  723.            ! WATER
    724. 

  724.            !CALL PSI_Suselj_Sood_2010(PSIM(I),PSIH(I),ZOL(I)) !produces neg TKE
    725. 

  725.            !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I))
    726. 

  726.            !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
    727. 

  727.            CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I))
    728. 

  728. 

  729.            !LOWER LIMIT ON PSI IN STABLE CONDITIONS
    730. 

  730.             psilim = -20.
    731. 

  731.         ELSE           
    732. 

  732.            ! LAND  
    733. 

  733.            !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I))
    734. 

  734.            !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
    735. 

  735.            !CALL PSI_Zilitinkevich_Esau_2007(PSIM(I),PSIH(I),ZOL(I))
    736. 

  736.            CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I))
    737. 

  737. 

  738.            !LOWER LIMIT ON PSI IN STABLE CONDITIONS
    739. 

  739.             psilim = -20.   !JOE: relaxed from -10 to -20.
    740. 

  740.         ENDIF      
    741. 

  741. 

  742.         PSIM(I)=AMAX1(PSIM(I),psilim)
    743. 

  743.         PSIH(I)=AMAX1(PSIH(I),psilim)
    744. 

  744.         PSIM10(I)=10./ZA(I)*PSIM(I)
    745. 

  745.         PSIM10(I)=AMAX1(PSIM10(I),psilim)
    746. 

  746.         PSIH10(I)=PSIM10(I)
    747. 

  747.         PSIM2(I)=2./ZA(I)*PSIM(I)
    748. 

  748.         PSIM2(I)=AMAX1(PSIM2(I),psilim)
    749. 

  749.         PSIH2(I)=PSIM2(I)    
    750. 

  750.         RMOL(I) = ZOL(I)/ZA(I) !1.0/L                                     
    751. 

  751. 

  752.      ELSEIF(BR(I) .GT. 0. .AND. BR(I) .LE. 0.2) THEN                                                 
    753. 

  753.         !========================================================
    754. 

  754.         !---CLASS 2; DAMPED MECHANICAL TURBULENCE:                                     
    755. 

  755.         !========================================================
    756. 

  756.         REGIME(I)=2.    
    757. 

  757.                                                        
    758. 

  758.         !COMPUTE z/L
    759. 

  759.         CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNT(I),zratio(I))
    760. 

  760. 

  761.         !COMPUTE PSIM and PSIH
    762. 

  762.         IF((XLAND(I)-1.5).GE.0)THEN                                            
    763. 

  763.            ! WATER
    764. 

  764.            !CALL PSI_Suselj_Sood_2010(PSIM(I),PSIH(I),ZOL(I))
    765. 

  765.            !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I))
    766. 

  766.            !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
    767. 

  767.            CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I))
    768. 

  768. 

  769.            !LOWER LIMIT ON PSI IN STABLE CONDITIONS
    770. 

  770.             psilim = -10.   !JOE: This limit is never hit.
    771. 

  771.         ELSE           
    772. 

  772.            ! LAND  
    773. 

  773.            !CALL PSI_Beljaars_Holtslag_1991(PSIM(I),PSIH(I),ZOL(I))
    774. 

  774.            !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
    775. 

  775.            !CALL PSI_Zilitinkevich_Esau_2007(PSIM(I),PSIH(I),ZOL(I))
    776. 

  776.            CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I))
    777. 

  777. 

  778.            !LOWER LIMIT ON PSI IN STABLE CONDITIONS
    779. 

  779.             psilim = -10.   !JOE: This limit is never hit.
    780. 

  780.         ENDIF              
    781. 

  781. 

  782.         ! LOWER LIMIT ON PSI IN STABLE CONDITIONS                                     
    783. 

  783.         PSIM(I)=AMAX1(PSIM(I),psilim)
    784. 

  784.         PSIH(I)=AMAX1(PSIH(I),psilim)
    785. 

  785.         PSIM10(I)=10./ZA(I)*PSIM(I)
    786. 

  786.         PSIM10(I)=AMAX1(PSIM10(I),psilim)                               
    787. 

  787.         PSIH10(I)=PSIM10(I)                                          
    788. 

  788.         PSIM2(I)=2./ZA(I)*PSIM(I)
    789. 

  789.         PSIM2(I)=AMAX1(PSIM2(I),psilim)                              
    790. 

  790.         PSIH2(I)=PSIM2(I)
    791. 

  791.         ! 1.0 over Monin-Obukhov length
    792. 

  792.         RMOL(I)= ZOL(I)/ZA(I)
    793. 

  793. 

  794.      ELSEIF(BR(I) .EQ. 0.) THEN                  
    795. 

  795.         !=========================================================  
    796. 

  796.         !-----CLASS 3; FORCED CONVECTION/NEUTRAL:                                                
    797. 

  797.         !=========================================================
    798. 

  798.         REGIME(I)=3.
    799. 

  799. 

  800.         PSIM(I)=0.0                                                              
    801. 

  801.         PSIH(I)=PSIM(I)                                                          
    802. 

  802.         PSIM10(I)=0.                                                   
    803. 

  803.         PSIH10(I)=PSIM10(I)                                           
    804. 

  804.         PSIM2(I)=0.                                                  
    805. 

  805.         PSIH2(I)=PSIM2(I)                                           
    806. 

  806.                                            
    807. 

  807.         !ZOL(I)=0.                                             
    808. 

  808.         IF(UST(I) .LT. 0.01)THEN                                                 
    809. 

  809.           ZOL(I)=BR(I)*GZ1OZ0(I)                                               
    810. 

  810.         ELSE                                                                     
    811. 

  811.           ZOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I)) 
    812. 

  812.         ENDIF                                                                    
    813. 

  813.         RMOL(I) = ZOL(I)/ZA(I)  
    814. 

  814. 

  815.      ELSEIF(BR(I) .LT. 0.)THEN            
    816. 

  816.         !==========================================================
    817. 

  817.         !-----CLASS 4; FREE CONVECTION:                                                  
    818. 

  818.         !==========================================================
    819. 

  819.         REGIME(I)=4.
    820. 

  820. 

  821.         !COMPUTE z/L
    822. 

  822.         CALL Li_etal_2010(ZOL(I),BR(I),ZA(I)/ZNT(I),zratio(I))
    823. 

  823. 

  824.         ZOL10=10./ZA(I)*ZOL(I)
    825. 

  825.         ZOL2=2./ZA(I)*ZOL(I)
    826. 

  826.         ZOL(I)=AMIN1(ZOL(I),0.)
    827. 

  827.         ZOL(I)=AMAX1(ZOL(I),-9.9999)
    828. 

  828.         ZOL10=AMIN1(ZOL10,0.)
    829. 

  829.         ZOL10=AMAX1(ZOL10,-9.9999)
    830. 

  830.         ZOL2=AMIN1(ZOL2,0.)
    831. 

  831.         ZOL2=AMAX1(ZOL2,-9.9999)
    832. 

  832.         NZOL=INT(-ZOL(I)*100.)
    833. 

  833.         RZOL=-ZOL(I)*100.-NZOL
    834. 

  834.         NZOL10=INT(-ZOL10*100.)
    835. 

  835.         RZOL10=-ZOL10*100.-NZOL10
    836. 

  836.         NZOL2=INT(-ZOL2*100.)
    837. 

  837.         RZOL2=-ZOL2*100.-NZOL2
    838. 

  838. 

  839.         !COMPUTE PSIM and PSIH
    840. 

  840.         IF((XLAND(I)-1.5).GE.0)THEN                                            
    841. 

  841.            ! WATER
    842. 

  842.            !CALL PSI_Suselj_Sood_2010(PSIM(I),PSIH(I),ZOL(I))
    843. 

  843.            !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I))
    844. 

  844.            !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
    845. 

  845.            CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I))
    846. 

  846.         ELSE           
    847. 

  847.            ! LAND  
    848. 

  848.            !CALL PSI_Hogstrom_1996(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I))
    849. 

  849.            !CALL PSI_Businger_1971(PSIM(I),PSIH(I),ZOL(I))
    850. 

  850.            CALL PSI_DyerHicks(PSIM(I),PSIH(I),ZOL(I), z_t(I), ZNT(I), ZA(I))
    851. 

  851.         ENDIF              
    852. 

  852. 

  853.         PSIM10(I)=PSIMTB(NZOL10)+RZOL10*(PSIMTB(NZOL10+1)-PSIMTB(NZOL10))
    854. 

  854.         PSIH10(I)=PSIHTB(NZOL10)+RZOL10*(PSIHTB(NZOL10+1)-PSIHTB(NZOL10))
    855. 

  855.         PSIM2(I)=PSIMTB(NZOL2)+RZOL2*(PSIMTB(NZOL2+1)-PSIMTB(NZOL2))
    856. 

  856.         PSIH2(I)=PSIHTB(NZOL2)+RZOL2*(PSIHTB(NZOL2+1)-PSIHTB(NZOL2))
    857. 

  857. 

  858.         !---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND
    859. 

  859.         !---HIGH ROUGHNESS.  THIS PREVENTS DENOMINATOR IN FLUXES
    860. 

  860.         !---FROM GETTING TOO SMALL
    861. 

  861.         PSIH(I)=AMIN1(PSIH(I),0.85*GZ1OZ0(I))
    862. 

  862.         PSIM(I)=AMIN1(PSIM(I),0.85*GZ1OZ0(I))
    863. 

  863.         PSIH2(I)=AMIN1(PSIH2(I),0.85*GZ2OZ0(I))
    864. 

  864.         PSIM10(I)=AMIN1(PSIM10(I),0.85*GZ10OZ0(I))
    865. 

  865. 

  866.         RMOL(I) = ZOL(I)/ZA(I)  
    867. 

  867. 

  868.      ENDIF
    869. 

  869. 

  870.   320 CONTINUE   
    871. 

  871. 

  872. !-----------------------------------------------------------
    873. 

  873. !-----COMPUTE U*, T*, AND SURFACE DIAGNOSTICS:                       
    874. 

  874. !-----------------------------------------------------------
    875. 

  875.   DO 330 I=its,ite
    876. 

  876. 

  877.      !------------------------------------------------------------
    878. 

  878.      !-----COMPUTE THE FRICTIONAL VELOCITY:                                           
    879. 

  879.      !     ZA(1982) EQS(2.60),(2.61).                                                 
    880. 

  880.      !------------------------------------------------------------
    881. 

  881.       DTG=THX(I)-THGB(I)                                                   
    882. 

  882.       PSIX=GZ1OZ0(I)-PSIM(I)                                                   
    883. 

  883.       PSIX10=GZ10OZ0(I)-PSIM10(I)
    884. 

  884.       ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE 
    885. 

  885.       UST(I)=0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX 
    886. 

  886.       IF((XLAND(I)-1.5).LT.0.)THEN        !LAND
    887. 

  887.           !JOE: UST(I)=AMAX1(UST(I),0.1)
    888. 

  888.           UST(I)=AMAX1(UST(I),0.01)  !Relaxing this limit
    889. 

  889.       ENDIF
    890. 

  890. 

  891.      !------------------------------------------------------------
    892. 

  892.      !-----COMPUTE THE RESISTANCE (PSIQ AND PSIT):                                           
    893. 

  893.      !------------------------------------------------------------
    894. 

  894.       ! LOWER LIMIT ADDED TO PREVENT LARGE FLHC IN SOIL MODEL
    895. 

  895.       ! ACTIVATES IN UNSTABLE CONDITIONS WITH THIN LAYERS OR HIGH Z0
    896. 

  896.       !PSIT=AMAX1(GZ1OZ0(I)-PSIH(I),2.)
    897. 

  897.       PSIT=AMAX1(GZ1OZt(I)-PSIH(I),2.)
    898. 

  898.       PSIT2=MAX(GZ2OZt(I)-PSIH2(I),2.)
    899. 

  899.                                   
    900. 

  900.       IF((XLAND(I)-1.5).GE.0)THEN                                            
    901. 

  901.          PSIQ=MAX(LOG(za(i)/z_q(I))-PSIH(I),2.) 
    902. 

  902.          PSIQ2=MAX(LOG(2./z_q(I))-PSIH2(I),2.)  
    903. 

  903.       ELSE                                                                     
    904. 

  904.          !CARLSON AND BOLAND (1978)(independent of zq):
    905. 

  905.          ZL=0.01                                                                
    906. 

  906.          PSIQ=MAX(LOG(KARMAN*UST(I)*ZA(I)/XKA + ZA(I)/ZL)-PSIH(I),2.0)   
    907. 

  907.          PSIQ2=MAX(LOG(KARMAN*UST(I)*2./XKA + 2./ZL)-PSIH2(I),2.0)                                  
    908. 

  908.       ENDIF                                                                    
    909. 

  909.                                  
    910. 

  910.       !------------------------------------------------------------
    911. 

  911.       !COMPUTE THE TEMPERATURE SCALE (or FRICTION TEMPERATURE, T*)
    912. 

  912.       !------------------------------------------------------------
    913. 

  913.       MOL(I)=KARMAN*DTG/PSIT/PRT
    914. 

  914.       !t_star(I) = -HFX(I)/(UST(I)*CPM(I)*RHOX(I))
    915. 

  915.       !t_star(I) = MOL(I)
    916. 

  916. 

  917.       !-----------------------------------------------------
    918. 

  918.       !COMPUTE 10 M WNDS
    919. 

  919.       ! If the lowest model level is close to 10-m, use it 
    920. 

  920.       ! instead of the flux-based formula.
    921. 

  921.       if (ZA(i) .gt. 7.0 .and. ZA(i) .lt. 13.0) then
    922. 

  922.          U10(I)=UX(I)                                   
    923. 

  923.          V10(I)=VX(I)
    924. 

  924.       else                                 
    925. 

  925.          U10(I)=UX(I)*PSIX10/PSIX                                    
    926. 

  926.          V10(I)=VX(I)*PSIX10/PSIX     
    927. 

  927.       endif                              
    928. 

  928. 

  929.       !-----------------------------------------------------
    930. 

  930.       !GET 2 M T, TH, AND Q
    931. 

  931.       !THESE WILL BE OVERWRITTEN FOR LAND POINTS IN THE LSM 
    932. 

  932.       TH2(I)=THGB(I)+DTG*PSIT2/PSIT                                
    933. 

  933.       Q2(I)=QSFC(I)+(QX(I)-QSFC(I))*PSIQ2/PSIQ                   
    934. 

  934.       T2(I) = TH2(I)*(PSFC(I)/100.)**ROVCP                     
    935. 

  935. 

  936.   330 CONTINUE                                                                   
    937. 

  937.                                                                                   
    938. 

  938. !-----------------------------------------------------------
    939. 

  939. !-----COMPUTE THE SURFACE SENSIBLE AND LATENT HEAT FLUXES:                       
    940. 

  940. !-----------------------------------------------------------
    941. 

  941.                                                                                  
    942. 

  942.       DO i=its,ite
    943. 

  943.         QFX(i) =0.                                                              
    944. 

  944.         HFX(i) =0.                                                              
    945. 

  945.       ENDDO
    946. 

  946. 

  947.       !-----------------------------------------------------
    948. 

  948.       !-- BUT FIRST, RECALCULATE ROUGHNESS LENGTHS (ZNT, Z_T, Z_Q).
    949. 

  949.       !-- USING AN UPDATED U*.          
    950. 

  950.       !-----------------------------------------------------
    951. 

  951.                                                                                  
    952. 

  952.       DO 360 I=its,ite
    953. 

  953. 

  954.          IF((XLAND(I)-1.5).GE.0)THEN                                            
    955. 

  955. 

  956.             !COMPUTE KINEMATIC VISCOSITY
    957. 

  957.             ! AHW: Garratt formula: Calculate kinematic viscosity 
    958. 

  958.             !      of air (linear approc to temp dependence at sea lev)
    959. 

  959.             VISC=(1.32+0.009*(T1D(I)-273.15))*1.E-5
    960. 

  960. 

  961.             !--------------------------------------
    962. 

  962.             ! WATER
    963. 

  963.             !--------------------------------------
    964. 

  964.             ! GET ZNT
    965. 

  965.             !--------------------------------------
    966. 

  966.             !CALL davis_etal_2008(ZNT(i),UST(i))
    967. 

  967.             !CALL Taylor_Yelland_2001(ZNT(i),UST(i),WSPD(i))
    968. 

  968.             CALL charnock_1955(ZNT(i),UST(i),WSPD(i),visc)
    969. 

  969. 

  970.             !COMPUTE ROUGHNESS REYNOLDS NUMBER WITH NEW ZNT
    971. 

  971.              restar=MAX(ust(i)*ZNT(i)/visc, 0.1)
    972. 

  972. 

  973.             !--------------------------------------
    974. 

  974.             ! GET z_t and z_q
    975. 

  975.             !--------------------------------------
    976. 

  976.             !CALL zilitinkevich_1995(ZNT(i),z_t(i),z_q(i),restar,UST(I),KARMAN,XLAND(I))
    977. 

  977.             !CALL garrat_1992(z_t(i),z_q(i),ZNT(i),restar,XLAND(I))
    978. 

  978.             CALL fairall_2001(z_t(i),z_q(i),restar,UST(i),visc)
    979. 

  979. 

  980.             zratio(i)=znt(i)/z_t(i)
    981. 

  981. 

  982.          ELSE
    983. 

  983. 

  984.             !--------------------------------------
    985. 

  985.             ! LAND
    986. 

  986.             !--------------------------------------
    987. 

  987. 

  988.             !COMPUTE ROUGHNESS REYNOLDS NUMBER (restar)
    989. 

  989.             VISC=(1.32+0.009*(T1D(I)-273.15))*1.E-5
    990. 

  990.             restar=MAX(ust(i)*ZNT(i)/visc, 0.1)
    991. 

  991. 

  992.             !--------------------------------------
    993. 

  993.             ! GET z_t and z_q
    994. 

  994.             !--------------------------------------
    995. 

  995.             !CALL zilitinkevich_1995(ZNT(i),z_t(i),z_q(i),restar,UST(I),KARMAN,XLAND(I))
    996. 

  996.             !CALL garrat_1992(z_t(i),z_q(i),ZNT(i),restar,XLAND(I))
    997. 

  997.             CALL Yang_2008(ZNT(i),z_t(i),z_q(i),UST(i),MOL(I),restar,visc,XLAND(I))
    998. 

  998. 

  999.             zratio(i)=znt(i)/z_t(i)
    1000. 

  1000. 

  1001.          ENDIF    
    1002. 

  1002. 

  1003.         !------------------------------------------
    1004. 

  1004.         ! CALCULATE THE EXCHANGE COEFFICIENTS FOR HEAT (FLHC)
    1005. 

  1005.         ! AND MOISTURE (FLQC)
    1006. 

  1006.         !------------------------------------------
    1007. 

  1007.         IF (ISFFLX.GT.0) THEN
    1008. 

  1008.            IF((XLAND(I)-1.5).GE.0)THEN                                            
    1009. 

  1009.               PSIQ=MAX(LOG(za(i)/z_q(I))-PSIH(I),2.) 
    1010. 

  1010.            ELSE                                                                     
    1011. 

  1011.               !CARLSON AND BOLAND (1978)(independent of zq):
    1012. 

  1012.               ZL=0.01                                                                
    1013. 

  1013.               PSIQ=MAX(LOG(KARMAN*UST(I)*ZA(I)/XKA + ZA(I)/ZL)-PSIH(I),2.0)   
    1014. 

  1014.            ENDIF                                                                    
    1015. 

  1015. 

  1016.            !FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/(   &
    1017. 

  1017.            !        ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-PSIH(I))
    1018. 

  1018. 

  1019.            !FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/(   &
    1020. 

  1020.            !     &             ALOG(za(i)/z_q(i))-PSIH(I))
    1021. 

  1021.            FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/PSIQ
    1022. 

  1022. 

  1023.            DTTHX=ABS(THX(I)-THGB(I))                                            
    1024. 

  1024.            IF(DTTHX.GT.1.E-5)THEN                                                   
    1025. 

  1025.              FLHC(I)=CPM(I)*RHOX(I)*UST(I)*MOL(I)/(THX(I)-THGB(I))          
    1026. 

  1026.            ELSE                                                                     
    1027. 

  1027.              FLHC(I)=0.                                                             
    1028. 

  1028.            ENDIF   
    1029. 

  1029.         ENDIF   
    1030. 

  1030. 

  1031. !        IF (I .eq. 2) THEN
    1032. 

  1032. !          write(*,1001)"REGIME:",REGIME(I)," z/L:",ZOL(I)," UST:",UST(I)," Tstar:",MOL(I)
    1033. 

  1033. !          write(*,1002)"PSIM:",PSIM(I)," PSIH:",PSIH(I)," FLHC:",FLHC(I)," FLQC:",FLQC(I)
    1034. 

  1034. !          write(*,1003)"CPM:",CPM(I)," RHOX:",RHOX(I)," L:",ZOL(I)/ZA(I)," DTHV:",THX(I)-THGB(I)
    1035. 

  1035. !          write(*,1004)"Z0/Zt:",zratio(I)," Z0:",ZNT(I)," Zt:",z_t(I)," za:",za(I)
    1036. 

  1036. !          write(*,1005)"Re:",restar," MAVAIL:",MAVAIL(I)," QSFC(I):",QSFC(I)," QX(I):",QX(I)
    1037. 

  1037. !          print*,"VISC=",VISC," Z0:",ZNT(I)," T1D(i):",T1D(i)
    1038. 

  1038. !          write(*,*)"============================================="
    1039. 

  1039. !        ENDIF
    1040. 

  1040.                          
    1041. 

  1041.   360 CONTINUE 
    1042. 

  1042.                                                    
    1043. 

  1043.  1001   format(A,F2.0, A,f10.4,A,f5.3, A,f11.5)
    1044. 

  1044.  1002   format(A,f7.2, A,f7.2, A,f9.2, A,f10.6)
    1045. 

  1045.  1003   format(A,f7.2, A,f7.2, A,f10.3,A,f10.3)
    1046. 

  1046.  1004   format(A,f11.3,A,f9.7, A,f9.7, A,f6.2, A,f10.3)
    1047. 

  1047.  1005   format(A,f9.2,A,f6.4,A,f7.4,A,f7.4)
    1048. 

  1048. 

  1049.   IF (ISFFLX .GT. 0) THEN                                                
    1050. 

  1050.      !----------------------------------
    1051. 

  1051.      !-----COMPUTE SURFACE MOISTURE FLUX:
    1052. 

  1052.      !                                                                                
    1053. 

  1053.      !IF(IDRY .EQ. 1)GOTO 390                                                
    1054. 

  1054. 

  1055.      DO I=its,ite
    1056. 

  1056.         QFX(I)=FLQC(I)*(QSFC(I)-QX(I))                                     
    1057. 

  1057.         !JOE: QFX(I)=AMAX1(QFX(I),0.)   !does not allows neg QFX           
    1058. 

  1058.         QFX(I)=AMAX1(QFX(I),-0.02)      !allows small neg QFX, like MYJ  
    1059. 

  1059.         LH(I)=XLV*QFX(I)
    1060. 

  1060.      ENDDO                                                             
    1061. 

  1061. 

  1062.   390 CONTINUE                                                                 
    1063. 

  1063.      !----------------------------------
    1064. 

  1064.      !-----COMPUTE SURFACE HEAT FLUX:
    1065. 

  1065.      !                                                                                
    1066. 

  1066.       DO 400 I=its,ite
    1067. 

  1067.         IF(XLAND(I)-1.5.GT.0.)THEN                                           
    1068. 

  1068.           HFX(I)=FLHC(I)*(THGB(I)-THX(I))                                
    1069. 

  1069.         ELSEIF(XLAND(I)-1.5.LT.0.)THEN                                       
    1070. 

  1070.           HFX(I)=FLHC(I)*(THGB(I)-THX(I))                                
    1071. 

  1071.           HFX(I)=AMAX1(HFX(I),-250.)                                       
    1072. 

  1072.         ENDIF                                                                  
    1073. 

  1073.   400 CONTINUE                                                                 
    1074. 

  1074.          
    1075. 

  1075.       DO I=its,ite
    1076. 

  1076. 

  1077. !         CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) &
    1078. 

  1078. !                /XKA+ZA(I)/ZL)-PSIH(I))
    1079. 

  1079. 

  1080.          CHS(I)=UST(I)*KARMAN/(LOG(ZA(I)/z_t(I))- &
    1081. 

  1081.               &PSIH(I))
    1082. 

  1082. 

  1083. ! The exchange coefficient for cloud water is assumed to be the same as
    1084. 

  1084. ! that for heat. CH is multiplied by WSPD.
    1085. 

  1085. 

  1086. !!         ch(i)=chs(i)
    1087. 

  1087.          ch(i)=flhc(i)/( cpm(i)*rhox(i) )
    1088. 

  1088. 

  1089. !        GZ2OZ0(I)=ALOG(2./ZNT(I))
    1090. 

  1090. !        PSIM2(I)=-10.*GZ2OZ0(I)
    1091. 

  1091. !        PSIM2(I)=AMAX1(PSIM2(I),-10.)
    1092. 

  1092. !        PSIH2(I)=PSIM2(I)
    1093. 

  1093. !         CQS2(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*2.0  &
    1094. 

  1094. !               /XKA+2.0/ZL)-PSIH2(I))
    1095. 

  1095. !         CHS2(I)=UST(I)*KARMAN/(GZ2OZ0(I)-PSIH2(I))
    1096. 

  1096. 

  1097.          !THESE ARE USED FOR 2-M DIAGNOSTICS ONLY
    1098. 

  1098.          CQS2(I)=UST(I)*KARMAN/(LOG(2.0/z_q(I))-PSIH2(I))
    1099. 

  1099.          CHS2(I)=UST(I)*KARMAN/(GZ2OZt(I)-PSIH2(I))
    1100. 

  1100. 

  1101.       ENDDO
    1102. 

  1102. 

  1103.    ENDIF
    1104. 

  1104. 

  1105.                                                                                 
    1106. 

  1106.    END SUBROUTINE SFCLAY1D_mynn
    1107. 

  1107. 

  1108. !-------------------------------------------------------------------          
    1109. 

  1109.    SUBROUTINE zilitinkevich_1995(Z_0,Zt,Zq,restar,ustar,KARMAN,landsea)
    1110. 

  1110. 

  1111.        IMPLICIT NONE
    1112. 

  1112.        REAL, INTENT(IN) :: Z_0,restar,ustar,KARMAN,landsea
    1113. 

  1113.        REAL, INTENT(OUT) :: Zt,Zq
    1114. 

  1114.        REAL, PARAMETER :: C1=0.075
    1115. 

  1115. 

  1116.        IF (landsea-1.5 .GT. 0) THEN    !WATER
    1117. 

  1117.           !THIS IS BASED ON Zilitinkevich, Grachev, and Fairall (2001;
    1118. 

  1118.           !Their equations 15 and 16).
    1119. 

  1119.           IF (restar .LT. 0.1) THEN
    1120. 

  1120.              Zt = Z_0*EXP(KARMAN*2.0)
    1121. 

  1121.              Zt = MIN( Zt, 5.0e-5)
    1122. 

  1122.              Zt = MAX( Zt, 2.0e-9) !same lower limit as ECMWF
    1123. 

  1123.              Zq = Z_0*EXP(KARMAN*3.0)
    1124. 

  1124.              Zq = MIN( Zq, 5.0e-5)
    1125. 

  1125.              Zq = MAX( Zq, 2.0e-9)
    1126. 

  1126.           ELSE
    1127. 

  1127.              Zt = Z_0*EXP(-KARMAN*(4.0*SQRT(restar)-3.2))
    1128. 

  1128.              Zt = MIN( Zt, 5.0e-5)
    1129. 

  1129.              Zt = MAX( Zt, 2.0e-9)
    1130. 

  1130.              Zq = Z_0*EXP(-KARMAN*(4.0*SQRT(restar)-4.2))
    1131. 

  1131.              Zq = MIN( Zt, 5.0e-5)
    1132. 

  1132.              Zq = MAX( Zt, 2.0e-9)
    1133. 

  1133.           ENDIF
    1134. 

  1134.        ELSE                             !LAND
    1135. 

  1135.           Zt = Z_0*EXP(-KARMAN*C1*SQRT(restar))
    1136. 

  1136.           Zt = MAX( Zt, Z_0/200.)
    1137. 

  1137.           Zt = MIN( Zt, Z_0)
    1138. 

  1138. 

  1139.           Zq = Z_0*EXP(-KARMAN*C1*SQRT(restar))
    1140. 

  1140.           Zq = MAX( Zq, Z_0/200.)
    1141. 

  1141.           Zq = MIN( Zq, Z_0)
    1142. 

  1142. 

  1143.           !Zq = Zt
    1144. 

  1144.        ENDIF
    1145. 

  1145.                    
    1146. 

  1146.        return
    1147. 

  1147. 

  1148.    END SUBROUTINE zilitinkevich_1995
    1149. 

  1149. !--------------------------------------------------------------------
    1150. 

  1150.    SUBROUTINE davis_etal_2008(Z_0,ustar)
    1151. 

  1151. 

  1152.     !This formulation for roughness length was designed to match 
    1153. 

  1153.     !the labratory experiments of Donelan et al. (2004).
    1154. 

  1154.     !This formulation produces smaller Z_0 between 5-20 m/s than
    1155. 

  1155.     !the Taylor_Yelland_2001 or Charnock subroutines (below).
    1156. 

  1156. 

  1157.        IMPLICIT NONE
    1158. 

  1158.        REAL, INTENT(IN)  :: ustar
    1159. 

  1159.        REAL, INTENT(OUT)  :: Z_0
    1160. 

  1160. 

  1161.        Z_0 = 10.*EXP(-10./(ustar**(1./3.)))
    1162. 

  1162.        Z_0 = MAX( Z_0, 1.27e-7)  !These max/mins were suggested by
    1163. 

  1163.        Z_0 = MIN( Z_0, 2.85e-3)  !Davis et al. (2008)
    1164. 

  1164.                    
    1165. 

  1165.        return
    1166. 

  1166. 

  1167.    END SUBROUTINE davis_etal_2008
    1168. 

  1168. !--------------------------------------------------------------------
    1169. 

  1169.    SUBROUTINE Taylor_Yelland_2001(Z_0,ustar,wsp10)
    1170. 

  1170. 

  1171.     !This formulation for roughness length was designed account for 
    1172. 

  1172.     !wave steepness.
    1173. 

  1173. 

  1174.        IMPLICIT NONE
    1175. 

  1175.        REAL, INTENT(IN)  :: ustar,wsp10
    1176. 

  1176.        REAL, INTENT(OUT) :: Z_0
    1177. 

  1177.        REAL, parameter  :: g=9.81, pi=3.14159265
    1178. 

  1178.        REAL :: hs, Tp, Lp
    1179. 

  1179. 

  1180.        !hs is the significant wave height
    1181. 

  1181.         hs = 0.0248*(wsp10**2.)
    1182. 

  1182.        !Tp dominant wave period
    1183. 

  1183.         Tp = 0.729*MAX(wsp10,0.1)
    1184. 

  1184.        !Lp is the wavelength of the dominant wave
    1185. 

  1185.         Lp = g*Tp**2/(2*pi)
    1186. 

  1186. 

  1187.        Z_0 = 1200.*hs*(hs/Lp)**4.5
    1188. 

  1188.        Z_0 = MAX( Z_0, 1.27e-7)  !These max/mins were suggested by
    1189. 

  1189.        Z_0 = MIN( Z_0, 2.85e-3)  !Davis et al. (2008)
    1190. 

  1190.                    
    1191. 

  1191.        return
    1192. 

  1192. 

  1193.    END SUBROUTINE Taylor_Yelland_2001
    1194. 

  1194. !--------------------------------------------------------------------
    1195. 

  1195.    SUBROUTINE charnock_1955(Z_0,ustar,wsp10,visc)
    1196. 

  1196.  
    1197. 

  1197.     !This version of Charnock's relation employs a varying
    1198. 

  1198.     !Charnock parameter, similar to COARE3.0 [Fairall et al. (2003)].
    1199. 

  1199.     !The Charnoch parameter CZC is varied from .011 to .018 
    1200. 

  1200.     !between 10-m wsp = 10 and 18. 
    1201. 

  1201. 

  1202.        IMPLICIT NONE
    1203. 

  1203.        REAL, INTENT(IN)  :: ustar, visc, wsp10
    1204. 

  1204.        REAL, INTENT(OUT) :: Z_0
    1205. 

  1205.        REAL, PARAMETER   :: G=9.81, CZO2=0.011
    1206. 

  1206.        REAL              :: CZC    !variable charnock "constant"   
    1207. 

  1207. 

  1208.        CZC = CZO2 + 0.007*MIN(MAX((wsp10-10.)/8., 0.), 1.0)
    1209. 

  1209.        Z_0 = CZC*ustar*ustar/G + (0.11*visc/MAX(ustar,0.1))
    1210. 

  1210.        Z_0 = MAX( Z_0, 1.27e-7)  !These max/mins were suggested by
    1211. 

  1211.        Z_0 = MIN( Z_0, 2.85e-3)  !Davis et al. (2008)
    1212. 

  1212. 

  1213.        return
    1214. 

  1214. 

  1215.    END SUBROUTINE charnock_1955
    1216. 

  1216. !--------------------------------------------------------------------
    1217. 

  1217.    SUBROUTINE garrat_1992(Zt,Zq,Z_0,Ren,landsea)
    1218. 

  1218. 

  1219.     !This formulation for the thermal and moisture roughness lengths 
    1220. 

  1220.     !(Zt and Zq) relates them to Z0 via the roughness Reynolds number (Ren).
    1221. 

  1221.     !This formula comes from Fairall et al. (2003). It is modified from
    1222. 

  1222.     !the original Garrat-Brutsaert model to better fit the COARE/HEXMAX
    1223. 

  1223.     !data. The formula for land uses a constant ratio (Z_0/7.4) taken
    1224. 

  1224.     !from Garrat (1992).
    1225. 

  1225. 

  1226.        IMPLICIT NONE
    1227. 

  1227.        REAL, INTENT(IN)  :: Ren, Z_0,landsea
    1228. 

  1228.        REAL, INTENT(OUT) :: Zt,Zq
    1229. 

  1229.        REAL :: Rq
    1230. 

  1230.        REAL, PARAMETER  :: e=2.71828183
    1231. 

  1231. 

  1232.        IF (landsea-1.5 .GT. 0) THEN    !WATER
    1233. 

  1233. 

  1234.           Zt = Z_0*EXP(2.0 - (2.48*(Ren**0.25)))
    1235. 

  1235.           Zq = Z_0*EXP(2.0 - (2.28*(Ren**0.25)))
    1236. 

  1236. 

  1237.           Zq = MIN( Zq, 5.0e-5)
    1238. 

  1238.           Zq = MAX( Zq, 2.0e-9)
    1239. 

  1239.           Zt = MIN( Zt, 5.0e-5)
    1240. 

  1240.           Zt = MAX( Zt, 2.0e-9) !same lower limit as ECMWF
    1241. 

  1241.        ELSE                            !LAND
    1242. 

  1242.           Zq = Z_0/(e**2.)      !taken from Garrat (1980,1992)
    1243. 

  1243.           Zt = Zq
    1244. 

  1244.        ENDIF
    1245. 

  1245.                    
    1246. 

  1246.        return
    1247. 

  1247. 

  1248.     END SUBROUTINE garrat_1992
    1249. 

  1249. !--------------------------------------------------------------------
    1250. 

  1250.     SUBROUTINE fairall_2001(Zt,Zq,Ren,ustar,visc)
    1251. 

  1251. 

  1252.     !This formulation for thermal and moisture roughness length (Zt and Zq) 
    1253. 

  1253.     !as a function of the roughness Reynolds number (Ren) comes from the 
    1254. 

  1254.     !COARE3.0 formulation, empirically derived from COARE and HEXMAX data
    1255. 

  1255.     ![Fairall et al. (2003)]. Edson et al. (2004; JGR) suspected that this
    1256. 

  1256.     !relationship overestimated roughness lengths for low Reynolds number 
    1257. 

  1257.     !flows, so a smooth flow relationship, taken from Garrat (1992, p. 102),
    1258. 

  1258.     !is used for flows with Ren < 2. 
    1259. 

  1259.     !
    1260. 

  1260.     !This is for use over water only.
    1261. 

  1261. 

  1262.        IMPLICIT NONE
    1263. 

  1263.        REAL, INTENT(IN)  :: Ren,ustar,visc
    1264. 

  1264.        REAL, INTENT(OUT) :: Zt,Zq
    1265. 

  1265. 

  1266.        IF (Ren .le. 2.) then
    1267. 

  1267. 

  1268.           Zt = (5.5e-5)*(Ren**(-0.63))
    1269. 

  1269.           !FOR SMOOTH SEAS, USE GARRAT FOR Zq
    1270. 

  1270.           Zq = 0.2*visc/MAX(ustar,0.1)
    1271. 

  1271. 

  1272.        ELSE
    1273. 

  1273.           
    1274. 

  1274.           !FOR ROUGH SEAS, USE FAIRALL
    1275. 

  1275.           Zt = (5.5e-5)*(Ren**(-0.63))
    1276. 

  1276.           Zq = Zt
    1277. 

  1277.  
    1278. 

  1278.        ENDIF
    1279. 

  1279. 

  1280.        Zt = MIN(Zt,5.5e-5)
    1281. 

  1281.        Zt = MAX(Zt,2.0e-9)
    1282. 

  1282. 

  1283.        Zq = MIN(Zt,5.5e-5)
    1284. 

  1284.        Zq = MAX(Zt,2.0e-9) 
    1285. 

  1285.                    
    1286. 

  1286.        return
    1287. 

  1287. 

  1288.     END SUBROUTINE fairall_2001
    1289. 

  1289. !--------------------------------------------------------------------
    1290. 

  1290.     SUBROUTINE Yang_2008(Z_0,Zt,Zq,ustar,tstar,Ren,visc,landsea)
    1291. 

  1291. 

  1292.     !This is a modified version of Yang et al (2002 QJRMS, 2008 JAMC) 
    1293. 

  1293.     !and Chen et al (2010, J of Hydromet). Although it was originally 
    1294. 

  1294.     !designed for arid regions with bare soil, it is modified 
    1295. 

  1295.     !here to perform over a broader spectrum of vegetation.
    1296. 

  1296.     !
    1297. 

  1297.     !The original formulation relates the thermal roughness length (Zt) 
    1298. 

  1298.     !to u* and T*:
    1299. 

  1299.     !  
    1300. 

  1300.     ! Zt = ht * EXP(-beta*(ustar**0.5)*(ABS(tstar)**0.25))
    1301. 

  1301.     !
    1302. 

  1302.     !where ht = Renc*visc/ustar and the critical Reynolds number 
    1303. 

  1303.     !(Renc) = 70. Beta was originally = 10 (2002 paper) but was revised 
    1304. 

  1304.     !to 7.2 (in 2008 paper). Their form typically varies the
    1305. 

  1305.     !ratio Z0/Zt by a few orders of magnitude (1-1E5).
    1306. 

  1306.     !
    1307. 

  1307.     !This modified form uses beta = 2.0, so zo/zt generally only varies by
    1308. 

  1308.     !10-300 over grass and cropland to about  20-2000 over forests.
    1309. 

  1309.     !
    1310. 

  1310.     !This should only be used over land!
    1311. 

  1311. 

  1312.        IMPLICIT NONE
    1313. 

  1313.        REAL, INTENT(IN)  :: Z_0, Ren, ustar, tstar, visc, landsea
    1314. 

  1314.        REAL              :: ht, tstar2
    1315. 

  1315.        REAL, INTENT(OUT) :: Zt,Zq
    1316. 

  1316.        REAL, PARAMETER  :: Renc=70., beta=2.0, e=2.71828183
    1317. 

  1317. 

  1318.        ht     = Renc*visc/MAX(ustar,0.01)
    1319. 

  1319.        tstar2 = MIN(tstar+0.2,0.0)
    1320. 

  1320. 

  1321.        Zt     = ht * EXP(-beta*(ustar**0.5)*(ABS(tstar2)**0.25))
    1322. 

  1322.        Zq     = Zt
    1323. 

  1323.                    
    1324. 

  1324.        Zt = MIN(Zt, Z_0/(e**2.))  !Garrat (1980,1992)
    1325. 

  1325.        Zq = MIN(Zq, Z_0/(e**2.)) 
    1326. 

  1326. 

  1327.        return
    1328. 

  1328. 

  1329.     END SUBROUTINE Yang_2008
    1330. 

  1330. !--------------------------------------------------------------------
    1331. 

  1331.     SUBROUTINE PSI_Hogstrom_1996(psi_m, psi_h, zL, Zt, Z_0, Za)
    1332. 

  1332. 

  1333.     ! This subroutine returns the stability functions based off
    1334. 

  1334.     ! of Hogstrom (1996).
    1335. 

  1335. 

  1336.        IMPLICIT NONE
    1337. 

  1337.        REAL, INTENT(IN)  :: zL, Zt, Z_0, Za
    1338. 

  1338.        REAL, INTENT(OUT) :: psi_m, psi_h
    1339. 

  1339.        REAL  :: x, x0, y, y0, zmL, zhL
    1340. 

  1340. 

  1341.        zmL = Z_0*zL/Za  
    1342. 

  1342.        zhL = Zt*zL/Za
    1343. 

  1343. 

  1344.        IF (zL .gt. 0.) THEN  !STABLE (not well tested - seem large)
    1345. 

  1345. 

  1346.           psi_m = -5.3*(zL - zmL)
    1347. 

  1347.           psi_h = -8.0*(zL - zhL)
    1348. 

  1348.  
    1349. 

  1349.        ELSE                 !UNSTABLE
    1350. 

  1350. 

  1351.           x = (1.-19.0*zL)**0.25
    1352. 

  1352.           x0= (1.-19.0*zmL)**0.25
    1353. 

  1353.           y = (1.-11.6*zL)**0.5
    1354. 

  1354.           y0= (1.-11.6*zhL)**0.5
    1355. 

  1355. 

  1356.           psi_m = 2.*LOG((1.+x)/(1.+x0)) + LOG((1.+x**2.)/(1.+x0**2.)) - &
    1357. 

  1357.                   2.*ATAN(x) + 2*ATAN(x0)
    1358. 

  1358.           psi_h = 2.*LOG((1.+y)/(1.+y0))
    1359. 

  1359. 

  1360.        ENDIF
    1361. 

  1361.                    
    1362. 

  1362.        return
    1363. 

  1363. 

  1364.     END SUBROUTINE PSI_Hogstrom_1996
    1365. 

  1365. !--------------------------------------------------------------------
    1366. 

  1366.     SUBROUTINE PSI_DyerHicks(psi_m, psi_h, zL, Zt, Z_0, Za)
    1367. 

  1367. 

  1368.     ! This subroutine returns the stability functions based off
    1369. 

  1369.     ! of Hogstrom (1996), but with different constants compatible
    1370. 

  1370.     ! with Dyer and Hicks (1970/74?). This formulation is used for
    1371. 

  1371.     ! testing/development by Nakanishi (personal communication).
    1372. 

  1372. 

  1373.        IMPLICIT NONE
    1374. 

  1374.        REAL, INTENT(IN)  :: zL, Zt, Z_0, Za
    1375. 

  1375.        REAL, INTENT(OUT) :: psi_m, psi_h
    1376. 

  1376.        REAL  :: x, x0, y, y0, zmL, zhL
    1377. 

  1377. 

  1378.        zmL = Z_0*zL/Za  
    1379. 

  1379.        zhL = Zt*zL/Za
    1380. 

  1380. 

  1381.        IF (zL .gt. 0.) THEN  !STABLE
    1382. 

  1382. 

  1383.           psi_m = -5.0*(zL - zmL)
    1384. 

  1384.           psi_h = -5.0*(zL - zhL)
    1385. 

  1385.  
    1386. 

  1386.        ELSE                 !UNSTABLE
    1387. 

  1387. 

  1388.           x = (1.-16.*zL)**0.25   
    1389. 

  1389.           x0= (1.-16.*zmL)**0.25
    1390. 

  1390.           !x = (1.-19.*zL)**0.25   !Hogstrom's form
    1391. 

  1391.           !x0= (1.-19.*zmL)**0.25
    1392. 

  1392. 

  1393.           y = (1.-16.*zL)**0.5
    1394. 

  1394.           y0= (1.-16.*zhL)**0.5
    1395. 

  1395. 

  1396.           psi_m = 2.*LOG((1.+x)/(1.+x0)) + LOG((1.+x**2.)/(1.+x0**2.)) - &
    1397. 

  1397.                   2.*ATAN(x) + 2*ATAN(x0)
    1398. 

  1398.           psi_h = 2.*LOG((1.+y)/(1.+y0))
    1399. 

  1399. 

  1400.        ENDIF
    1401. 

  1401.                    
    1402. 

  1402.        return
    1403. 

  1403. 

  1404.     END SUBROUTINE PSI_DyerHicks
    1405. 

  1405. !--------------------------------------------------------------------
    1406. 

  1406.     SUBROUTINE PSI_Beljaars_Holtslag_1991(psi_m, psi_h, zL)
    1407. 

  1407. 

  1408.     ! This subroutine returns the stability functions based off
    1409. 

  1409.     ! of Beljaar and Holtslag 1991, which is an extension of Holtslag
    1410. 

  1410.     ! and Debruin 1989.
    1411. 

  1411. 

  1412.        IMPLICIT NONE
    1413. 

  1413.        REAL, INTENT(IN)  :: zL
    1414. 

  1414.        REAL, INTENT(OUT) :: psi_m, psi_h
    1415. 

  1415.        REAL, PARAMETER  :: a=1., b=0.666, c=5., d=0.35
    1416. 

  1416. 

  1417.        IF (zL .lt. 0.) THEN  !UNSTABLE
    1418. 

  1418. 

  1419.           WRITE(*,*)"WARNING: Universal stability function from"
    1420. 

  1420.           WRITE(*,*)"        Beljaars and Holtslag (1991) should only"
    1421. 

  1421.           WRITE(*,*)"        be used in the stable regime!"
    1422. 

  1422.           psi_m = 0.
    1423. 

  1423.           psi_h = 0.
    1424. 

  1424.  
    1425. 

  1425.        ELSE                 !STABLE
    1426. 

  1426. 

  1427.           psi_m = -(a*zL + b*(zL -(c/d))*exp(-d*zL) + (b*c/d))
    1428. 

  1428.           psi_h = -((1.+.666*a*zL)**1.5 + &
    1429. 

  1429.                   b*(zL - (c/d))*exp(-d*zL) + (b*c/d) -1.)
    1430. 

  1430. 

  1431.        ENDIF
    1432. 

  1432.                    
    1433. 

  1433.        return
    1434. 

  1434. 

  1435.     END SUBROUTINE PSI_Beljaars_Holtslag_1991
    1436. 

  1436. !--------------------------------------------------------------------
    1437. 

  1437.     SUBROUTINE PSI_Zilitinkevich_Esau_2007(psi_m, psi_h, zL)
    1438. 

  1438. 

  1439.     ! This subroutine returns the stability functions come from
    1440. 

  1440.     ! Zilitinkevich and Esau (2007, BM), which are formulatioed from the
    1441. 

  1441.     ! "generalized similarity theory" and tuned to the LES DATABASE64
    1442. 

  1442.     ! to determine their dependence on z/L.
    1443. 

  1443. 

  1444.        IMPLICIT NONE
    1445. 

  1445.        REAL, INTENT(IN)  :: zL
    1446. 

  1446.        REAL, INTENT(OUT) :: psi_m, psi_h
    1447. 

  1447.        REAL, PARAMETER  :: Cm=3.0, Ct=2.5
    1448. 

  1448. 

  1449.        IF (zL .lt. 0.) THEN  !UNSTABLE
    1450. 

  1450. 

  1451.           WRITE(*,*)"WARNING: Universal stability function from"
    1452. 

  1452.           WRITE(*,*)"        Zilitinkevich and Esau (2007) should only"
    1453. 

  1453.           WRITE(*,*)"        be used in the stable regime!"
    1454. 

  1454.           psi_m = 0.
    1455. 

  1455.           psi_h = 0.
    1456. 

  1456.  
    1457. 

  1457.        ELSE                 !STABLE
    1458. 

  1458. 

  1459.           psi_m = -Cm*(zL**(5./6.))
    1460. 

  1460.           psi_h = -Ct*(zL**(4./5.))
    1461. 

  1461. 

  1462.        ENDIF
    1463. 

  1463.                    
    1464. 

  1464.        return
    1465. 

  1465. 

  1466.     END SUBROUTINE PSI_Zilitinkevich_Esau_2007
    1467. 

  1467. !--------------------------------------------------------------------
    1468. 

  1468.     SUBROUTINE PSI_Businger_1971(psi_m, psi_h, zL)
    1469. 

  1469. 

  1470.     ! This subroutine returns the flux-profile relationships
    1471. 

  1471.     ! of Businger el al. 1971.
    1472. 

  1472. 

  1473.        IMPLICIT NONE
    1474. 

  1474.        REAL, INTENT(IN)  :: zL
    1475. 

  1475.        REAL, INTENT(OUT) :: psi_m, psi_h
    1476. 

  1476.        REAL  :: x, y
    1477. 

  1477.        REAL, PARAMETER  ::  Pi180 = 3.14159265/180.
    1478. 

  1478. 

  1479.        IF (zL .lt. 0.) THEN  !UNSTABLE
    1480. 

  1480. 

  1481.           x = (1. - 15.0*zL)**0.25
    1482. 

  1482.           y = (1. - 9.0*zL)**0.5
    1483. 

  1483. 

  1484.           psi_m = LOG(((1.+x)/2.)**2.) + LOG((1.+x**2.)/2.) - &
    1485. 

  1485.                            2.*ATAN(x) + Pi180*90.
    1486. 

  1486.           psi_h = 2.*LOG((1.+y)/2.)
    1487. 

  1487. 

  1488.        ELSE                 !STABLE
    1489. 

  1489. 

  1490.           psi_m = -4.7*zL
    1491. 

  1491.           psi_h = -(4.7/0.74)*zL
    1492. 

  1492. 

  1493.        ENDIF
    1494. 

  1494.                    
    1495. 

  1495.        return
    1496. 

  1496. 

  1497.     END SUBROUTINE PSI_Businger_1971
    1498. 

  1498. !--------------------------------------------------------------------
    1499. 

  1499.     SUBROUTINE PSI_Suselj_Sood_2010(psi_m, psi_h, zL)
    1500. 

  1500. 

  1501.     !This subroutine returns flux-profile relatioships based off
    1502. 

  1502.     !of Lobocki (1993), which is derived from the MY-level 2 model.
    1503. 

  1503.     !Suselj and Sood (2010) applied the surface layer length scales
    1504. 

  1504.     !from Nakanishi (2001) to get this new relationship. These functions
    1505. 

  1505.     !are more agressive (larger magnitude) than most formulations. They
    1506. 

  1506.     !showed improvement over water, but untested over land.
    1507. 

  1507. 

  1508.        IMPLICIT NONE
    1509. 

  1509.        REAL, INTENT(IN)  :: zL
    1510. 

  1510.        REAL, INTENT(OUT) :: psi_m, psi_h
    1511. 

  1511.        REAL, PARAMETER  :: Rfc=0.19, Ric=0.183, PHIT=0.8
    1512. 

  1512. 

  1513.        IF (zL .gt. 0.) THEN  !STABLE
    1514. 

  1514. 

  1515.           psi_m = -(zL/Rfc + 1.1223*EXP(1.-1.6666/zL))
    1516. 

  1516.           !psi_h = -zL*Ric/((Rfc**2.)*PHIT) + 8.209*(zL**1.1091)
    1517. 

  1517.           !THEIR EQ FOR PSI_H CRASHES THE MODEL AND DOES NOT MATCH
    1518. 

  1518.           !THEIR FIG 1. THIS EQ (BELOW) MATCHES THEIR FIG 1 BETTER:
    1519. 

  1519.           psi_h = -(zL*Ric/((Rfc**2.)*5.) + 7.09*(zL**1.1091))
    1520. 

  1520.  
    1521. 

  1521.        ELSE                 !UNSTABLE
    1522. 

  1522. 

  1523.           psi_m = 0.9904*LOG(1. - 14.264*zL)
    1524. 

  1524.           psi_h = 1.0103*LOG(1. - 16.3066*zL)
    1525. 

  1525. 

  1526.        ENDIF
    1527. 

  1527.                    
    1528. 

  1528.        return
    1529. 

  1529. 

  1530.     END SUBROUTINE PSI_Suselj_Sood_2010
    1531. 

  1531. !--------------------------------------------------------------------
    1532. 

  1532.     SUBROUTINE Li_etal_2010(zL, Rib, zaz0, z0zt)
    1533. 

  1533. 

  1534.     !This subroutine returns a more robust z/L that best matches
    1535. 

  1535.     !the z/L from Hogstrom (1996) for unstable conditions and Beljaars
    1536. 

  1536.     !and Holtslag (1991) for stable conditions.
    1537. 

  1537. 

  1538.        IMPLICIT NONE
    1539. 

  1539.        REAL, INTENT(OUT)  :: zL
    1540. 

  1540.        REAL, INTENT(IN) :: Rib, zaz0, z0zt
    1541. 

  1541.        REAL :: alfa, beta, zaz02, z0zt2
    1542. 

  1542.        REAL, PARAMETER  :: au11=0.045, bu11=0.003, bu12=0.0059, &
    1543. 

  1543.                            bu21=-0.0828, bu22=0.8845, bu31=0.1739, &
    1544. 

  1544.                            bu32=-0.9213, bu33=-0.1057
    1545. 

  1545.        REAL, PARAMETER  :: aw11=0.5738, aw12=-0.4399, aw21=-4.901,&
    1546. 

  1546.                            aw22=52.50, bw11=-0.0539, bw12=1.540, &
    1547. 

  1547.                            bw21=-0.669, bw22=-3.282
    1548. 

  1548.        REAL, PARAMETER  :: as11=0.7529, as21=14.94, bs11=0.1569,&
    1549. 

  1549.                            bs21=-0.3091, bs22=-1.303
    1550. 

  1550.           
    1551. 

  1551.        !set limits according to Li et al (2010), p 157.
    1552. 

  1552.        zaz02=zaz0
    1553. 

  1553.        IF (zaz0 .lt. 100.0) zaz02=100.
    1554. 

  1554.        IF (zaz0 .gt. 100000.0) zaz02=100000.
    1555. 

  1555. 

  1556.        !set more limits according to Li et al (2010)
    1557. 

  1557.        z0zt2=z0zt
    1558. 

  1558.        IF (z0zt .lt. 0.5) z0zt2=0.5
    1559. 

  1559.        IF (z0zt .gt. 100.0) z0zt2=100.
    1560. 

  1560. 

  1561.        alfa = LOG(zaz02)
    1562. 

  1562.        beta = LOG(z0zt2)
    1563. 

  1563. 

  1564.        IF (Rib .le. 0.0) THEN
    1565. 

  1565.           zL = au11*alfa*Rib**2 + (                   &
    1566. 

  1566.                   (bu11*beta + bu12)*alfa**2 +        &
    1567. 

  1567.                   (bu21*beta + bu22)*alfa    +        &
    1568. 

  1568.                   (bu31*beta**2 + bu32*beta + bu33))*Rib
    1569. 

  1569.           !if(zL .LT. -15 .OR. zl .GT. 0.)print*,"VIOLATION Rib<0:",zL
    1570. 

  1570.           zL = MAX(zL,-15.) !LIMITS SET ACCORDING TO Li et al (2010)
    1571. 

  1571.           zL = MIN(zL,0.)   !Figure 1.
    1572. 

  1572.        ELSEIF (Rib .gt. 0.0 .AND. Rib .le. 0.2) THEN
    1573. 

  1573.           zL = ((aw11*beta + aw12)*alfa +             &
    1574. 

  1574.                 (aw21*beta + aw22))*Rib**2 +          &
    1575. 

  1575.                ((bw11*beta + bw12)*alfa +             &
    1576. 

  1576.                 (bw21*beta + bw22))*Rib
    1577. 

  1577.           !if(zL .LT. 0 .OR. zl .GT. 4)print*,"VIOLATION 0<Rib<0.2:",zL
    1578. 

  1578.           zL = MIN(zL,4.) !LIMITS APPROX SET ACCORDING TO Li et al (2010)
    1579. 

  1579.           zL = MAX(zL,0.) !THEIR FIGURE 1B.
    1580. 

  1580.        ELSE
    1581. 

  1581.           zL = (as11*alfa + as21)*Rib + bs11*alfa +   &
    1582. 

  1582.                 bs21*beta + bs22
    1583. 

  1583.           !if(zL .LE. 1 .OR. zl .GT. 23)print*,"VIOLATION Rib>0.2:",zL
    1584. 

  1584.           zL = MIN(zL,20.) !LIMITS ACCORDING TO Li et al (2010), THIER
    1585. 

  1585.                            !FIGUE 1C.
    1586. 

  1586.           zL = MAX(zL,1.)
    1587. 

  1587.        ENDIF
    1588. 

  1588. 

  1589.        return
    1590. 

  1590. 

  1591.     END SUBROUTINE Li_etal_2010
    1592. 

  1592. !--------------------------------------------------------------------
    1593. 

  1593. 

  1594. 

  1595. END MODULE module_sf_mynn
    1596. 

  1596.