PageRenderTime 81ms CodeModel.GetById 39ms RepoModel.GetById 1ms app.codeStats 0ms

/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

Large files files are truncated, but you can click here to view the full file

    

  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.

Large files files are truncated, but you can click here to view the full file