PageRenderTime 74ms CodeModel.GetById 34ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_sf_sfclayrev.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1091 lines | 624 code | 98 blank | 369 comment | 18 complexity | 5fd19da281f6a973f0b95d62160137a0 MD5 | raw file
Possible License(s): AGPL-1.0 
    

  1. !WRF:MODEL_LAYER:PHYSICS
    2. 

  2. !
    3. 

  3. MODULE module_sf_sfclayrev
    4. 

  4. 

  5.  REAL    , PARAMETER ::  VCONVC=1.
    6. 

  6.  REAL    , PARAMETER ::  CZO=0.0185
    7. 

  7.  REAL    , PARAMETER ::  OZO=1.59E-5
    8. 

  8. 

  9.  REAL,   DIMENSION(0:1000 ),SAVE          :: PSIMTB,PSIHTB
    10. 

  10. 

  11. CONTAINS
    12. 

  12. 

  13. !-------------------------------------------------------------------
    14. 

  14.    SUBROUTINE SFCLAYREV(U3D,V3D,T3D,QV3D,P3D,dz8w,                    &
    15. 

  15.                      CP,G,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM,      &
    16. 

  16.                      ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
    17. 

  17.                      XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, &
    18. 

  18.                      U10,V10,TH2,T2,Q2,                            &
    19. 

  19.                      GZ1OZ0,WSPD,BR,ISFFLX,DX,                     &
    20. 

  20.                      SVP1,SVP2,SVP3,SVPT0,EP1,EP2,                 &
    21. 

  21.                      KARMAN,EOMEG,STBOLT,                          &
    22. 

  22.                      P1000mb,                                      &
    23. 

  23.                      ids,ide, jds,jde, kds,kde,                    &
    24. 

  24.                      ims,ime, jms,jme, kms,kme,                    &
    25. 

  25.                      its,ite, jts,jte, kts,kte,                    &
    26. 

  26.                      ustm,ck,cka,cd,cda,isftcflx,iz0tlnd           )
    27. 

  27. !-------------------------------------------------------------------
    28. 

  28.       IMPLICIT NONE
    29. 

  29. !-------------------------------------------------------------------
    30. 

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

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

  32. !-- T3D         temperature (K)
    33. 

  33. !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
    34. 

  34. !-- P3D         3D pressure (Pa)
    35. 

  35. !-- dz8w        dz between full levels (m)
    36. 

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

  37. !-- G           acceleration due to gravity (m/s^2)
    38. 

  38. !-- ROVCP       R/CP
    39. 

  39. !-- R           gas constant for dry air (J/kg/K)
    40. 

  40. !-- XLV         latent heat of vaporization for water (J/kg)
    41. 

  41. !-- PSFC        surface pressure (Pa)
    42. 

  42. !-- ZNT         roughness length (m)
    43. 

  43. !-- UST         u* in similarity theory (m/s)
    44. 

  44. !-- USTM        u* in similarity theory (m/s) without vconv correction
    45. 

  45. !               used to couple with TKE scheme
    46. 

  46. !-- PBLH        PBL height from previous time (m)
    47. 

  47. !-- MAVAIL      surface moisture availability (between 0 and 1)
    48. 

  48. !-- ZOL         z/L height over Monin-Obukhov length
    49. 

  49. !-- MOL         T* (similarity theory) (K)
    50. 

  50. !-- REGIME      flag indicating PBL regime (stable, unstable, etc.)
    51. 

  51. !-- PSIM        similarity stability function for momentum
    52. 

  52. !-- PSIH        similarity stability function for heat
    53. 

  53. !-- XLAND       land mask (1 for land, 2 for water)
    54. 

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

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

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

  57. !-- TSK         surface temperature (K)
    58. 

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

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

  60. !-- CHS         heat/moisture exchange coefficient for LSM (m/s)
    61. 

  61. !-- QGH         lowest-level saturated mixing ratio
    62. 

  62. !-- QSFC        ground saturated mixing ratio
    63. 

  63. !-- U10         diagnostic 10m u wind
    64. 

  64. !-- V10         diagnostic 10m v wind
    65. 

  65. !-- TH2         diagnostic 2m theta (K)
    66. 

  66. !-- T2          diagnostic 2m temperature (K)
    67. 

  67. !-- Q2          diagnostic 2m mixing ratio (kg/kg)
    68. 

  68. !-- GZ1OZ0      log(z/z0) where z0 is roughness length
    69. 

  69. !-- WSPD        wind speed at lowest model level (m/s)
    70. 

  70. !-- BR          bulk Richardson number in surface layer
    71. 

  71. !-- ISFFLX      isfflx=1 for surface heat and moisture fluxes
    72. 

  72. !-- DX          horizontal grid size (m)
    73. 

  73. !-- SVP1        constant for saturation vapor pressure (kPa)
    74. 

  74. !-- SVP2        constant for saturation vapor pressure (dimensionless)
    75. 

  75. !-- SVP3        constant for saturation vapor pressure (K)
    76. 

  76. !-- SVPT0       constant for saturation vapor pressure (K)
    77. 

  77. !-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
    78. 

  78. !-- EP2         constant for specific humidity calculation 
    79. 

  79. !               (R_d/R_v) (dimensionless)
    80. 

  80. !-- KARMAN      Von Karman constant
    81. 

  81. !-- EOMEG       angular velocity of earth's rotation (rad/s)
    82. 

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

  83. !-- ck          enthalpy exchange coeff at 10 meters
    84. 

  84. !-- cd          momentum exchange coeff at 10 meters
    85. 

  85. !-- cka         enthalpy exchange coeff at the lowest model level
    86. 

  86. !-- cda         momentum exchange coeff at the lowest model level
    87. 

  87. !-- isftcflx    =0, (Charnock and Carlson-Boland); =1, AHW Ck, Cd, =2 Garratt
    88. 

  88. !-- iz0tlnd     =0 Carlson-Boland, =1 Czil_new
    89. 

  89. !-- ids         start index for i in domain
    90. 

  90. !-- ide         end index for i in domain
    91. 

  91. !-- jds         start index for j in domain
    92. 

  92. !-- jde         end index for j in domain
    93. 

  93. !-- kds         start index for k in domain
    94. 

  94. !-- kde         end index for k in domain
    95. 

  95. !-- ims         start index for i in memory
    96. 

  96. !-- ime         end index for i in memory
    97. 

  97. !-- jms         start index for j in memory
    98. 

  98. !-- jme         end index for j in memory
    99. 

  99. !-- kms         start index for k in memory
    100. 

  100. !-- kme         end index for k in memory
    101. 

  101. !-- its         start index for i in tile
    102. 

  102. !-- ite         end index for i in tile
    103. 

  103. !-- jts         start index for j in tile
    104. 

  104. !-- jte         end index for j in tile
    105. 

  105. !-- kts         start index for k in tile
    106. 

  106. !-- kte         end index for k in tile
    107. 

  107. !-------------------------------------------------------------------
    108. 

  108.       INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
    109. 

  109.                                         ims,ime, jms,jme, kms,kme, &
    110. 

  110.                                         its,ite, jts,jte, kts,kte
    111. 

  111. !                                                               
    112. 

  112.       INTEGER,  INTENT(IN )   ::        ISFFLX
    113. 

  113.       REAL,     INTENT(IN )   ::        SVP1,SVP2,SVP3,SVPT0
    114. 

  114.       REAL,     INTENT(IN )   ::        EP1,EP2,KARMAN,EOMEG,STBOLT
    115. 

  115.       REAL,     INTENT(IN )   ::        P1000mb
    116. 

  116. !
    117. 

  117.       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
    118. 

  118.                 INTENT(IN   )   ::                           dz8w
    119. 

  119.                                         
    120. 

  120.       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
    121. 

  121.                 INTENT(IN   )   ::                           QV3D, &
    122. 

  122.                                                               P3D, &
    123. 

  123.                                                               T3D
    124. 

  124. 

  125.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    126. 

  126.                 INTENT(IN   )               ::             MAVAIL, &
    127. 

  127.                                                              PBLH, &
    128. 

  128.                                                             XLAND, &
    129. 

  129.                                                               TSK
    130. 

  130.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    131. 

  131.                 INTENT(OUT  )               ::                U10, &
    132. 

  132.                                                               V10, &
    133. 

  133.                                                               TH2, &
    134. 

  134.                                                                T2, &
    135. 

  135.                                                                Q2, &
    136. 

  136.                                                              QSFC
    137. 

  137. 

  138. !
    139. 

  139.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    140. 

  140.                 INTENT(INOUT)               ::             REGIME, &
    141. 

  141.                                                               HFX, &
    142. 

  142.                                                               QFX, &
    143. 

  143.                                                                LH, &
    144. 

  144.                                                           MOL,RMOL
    145. 

  145. !m the following 5 are change to memory size
    146. 

  146. !
    147. 

  147.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    148. 

  148.                 INTENT(INOUT)   ::                 GZ1OZ0,WSPD,BR, &
    149. 

  149.                                                         PSIM,PSIH
    150. 

  150. 

  151.       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
    152. 

  152.                 INTENT(IN   )   ::                            U3D, &
    153. 

  153.                                                               V3D
    154. 

  154.                                         
    155. 

  155.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    156. 

  156.                 INTENT(IN   )               ::               PSFC
    157. 

  157. 

  158.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    159. 

  159.                 INTENT(INOUT)   ::                            ZNT, &
    160. 

  160.                                                               ZOL, &
    161. 

  161.                                                               UST, &
    162. 

  162.                                                               CPM, &
    163. 

  163.                                                              CHS2, &
    164. 

  164.                                                              CQS2, &
    165. 

  165.                                                               CHS
    166. 

  166. 

  167.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    168. 

  168.                 INTENT(INOUT)   ::                      FLHC,FLQC
    169. 

  169. 

  170.       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
    171. 

  171.                 INTENT(INOUT)   ::                                 &
    172. 

  172.                                                               QGH
    173. 

  173. 

  174. 

  175.                                     
    176. 

  176.       REAL,     INTENT(IN   )               ::   CP,G,ROVCP,R,XLV,DX
    177. 

  177.  
    178. 

  178.       REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme )              , &
    179. 

  179.                 INTENT(OUT)     ::              ck,cka,cd,cda,ustm
    180. 

  180. 

  181.       INTEGER,  OPTIONAL,  INTENT(IN )   ::     ISFTCFLX, IZ0TLND
    182. 

  182. ! LOCAL VARS
    183. 

  183. 

  184.       REAL,     DIMENSION( its:ite ) ::                       U1D, &
    185. 

  185.                                                               V1D, &
    186. 

  186.                                                              QV1D, &
    187. 

  187.                                                               P1D, &
    188. 

  188.                                                               T1D
    189. 

  189. 

  190.       REAL,     DIMENSION( its:ite ) ::                    dz8w1d
    191. 

  191. 

  192.       INTEGER ::  I,J
    193. 

  193. 

  194.       DO J=jts,jte
    195. 

  195.         DO i=its,ite
    196. 

  196.           dz8w1d(I) = dz8w(i,1,j)
    197. 

  197.         ENDDO
    198. 

  198.    
    199. 

  199.         DO i=its,ite
    200. 

  200.            U1D(i) =U3D(i,1,j)
    201. 

  201.            V1D(i) =V3D(i,1,j)
    202. 

  202.            QV1D(i)=QV3D(i,1,j)
    203. 

  203.            P1D(i) =P3D(i,1,j)
    204. 

  204.            T1D(i) =T3D(i,1,j)
    205. 

  205.         ENDDO
    206. 

  206. 

  207.         !  Sending array starting locations of optional variables may cause
    208. 

  208.         !  troubles, so we explicitly change the call.
    209. 

  209. 

  210.         CALL SFCLAYREV1D(J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,               &
    211. 

  211.                 CP,G,ROVCP,R,XLV,PSFC(ims,j),CHS(ims,j),CHS2(ims,j),&
    212. 

  212.                 CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j),   &
    213. 

  213.                 ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j),    &
    214. 

  214.                 MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j),  &
    215. 

  215.                 XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j),     &
    216. 

  216.                 U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j),        &
    217. 

  217.                 Q2(ims,j),FLHC(ims,j),FLQC(ims,j),QGH(ims,j),      &
    218. 

  218.                 QSFC(ims,j),LH(ims,j),                             &
    219. 

  219.                 GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX,     &
    220. 

  220.                 SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,  &
    221. 

  221.                 P1000mb,                                           &
    222. 

  222.                 ids,ide, jds,jde, kds,kde,                         &
    223. 

  223.                 ims,ime, jms,jme, kms,kme,                         &
    224. 

  224.                 its,ite, jts,jte, kts,kte                          &
    225. 

  225. #if ( EM_CORE == 1 )
    226. 

  226.                 ,isftcflx,iz0tlnd,                                 &
    227. 

  227.                 USTM(ims,j),CK(ims,j),CKA(ims,j),                  &
    228. 

  228.                 CD(ims,j),CDA(ims,j)                               &
    229. 

  229. #endif
    230. 

  230.                                                                    )
    231. 

  231.       ENDDO
    232. 

  232. 

  233. 

  234.    END SUBROUTINE SFCLAYREV
    235. 

  235. 

  236. 

  237. !-------------------------------------------------------------------
    238. 

  238.    SUBROUTINE SFCLAYREV1D(J,UX,VX,T1D,QV1D,P1D,dz8w1d,                &
    239. 

  239.                      CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, &
    240. 

  240.                      ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH,      &
    241. 

  241.                      XLAND,HFX,QFX,TSK,                            &
    242. 

  242.                      U10,V10,TH2,T2,Q2,FLHC,FLQC,QGH,              &
    243. 

  243.                      QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX,             &
    244. 

  244.                      SVP1,SVP2,SVP3,SVPT0,EP1,EP2,                 &
    245. 

  245.                      KARMAN,EOMEG,STBOLT,                          &
    246. 

  246.                      P1000mb,                                      &
    247. 

  247.                      ids,ide, jds,jde, kds,kde,                    &
    248. 

  248.                      ims,ime, jms,jme, kms,kme,                    &
    249. 

  249.                      its,ite, jts,jte, kts,kte,                    &
    250. 

  250.                      isftcflx, iz0tlnd,                            &
    251. 

  251.                      ustm,ck,cka,cd,cda                            )
    252. 

  252. !-------------------------------------------------------------------
    253. 

  253.       IMPLICIT NONE
    254. 

  254. !-------------------------------------------------------------------
    255. 

  255.       REAL,     PARAMETER     ::        XKA=2.4E-5
    256. 

  256.       REAL,     PARAMETER     ::        PRT=1.
    257. 

  257. 

  258.       INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
    259. 

  259.                                         ims,ime, jms,jme, kms,kme, &
    260. 

  260.                                         its,ite, jts,jte, kts,kte, &
    261. 

  261.                                         J
    262. 

  262. !                                                               
    263. 

  263.       INTEGER,  INTENT(IN )   ::        ISFFLX
    264. 

  264.       REAL,     INTENT(IN )   ::        SVP1,SVP2,SVP3,SVPT0
    265. 

  265.       REAL,     INTENT(IN )   ::        EP1,EP2,KARMAN,EOMEG,STBOLT
    266. 

  266.       REAL,     INTENT(IN )   ::        P1000mb
    267. 

  267. 

  268. !
    269. 

  269.       REAL,     DIMENSION( ims:ime )                             , &
    270. 

  270.                 INTENT(IN   )               ::             MAVAIL, &
    271. 

  271.                                                              PBLH, &
    272. 

  272.                                                             XLAND, &
    273. 

  273.                                                               TSK
    274. 

  274. !
    275. 

  275.       REAL,     DIMENSION( ims:ime )                             , &
    276. 

  276.                 INTENT(IN   )               ::             PSFCPA
    277. 

  277. 

  278.       REAL,     DIMENSION( ims:ime )                             , &
    279. 

  279.                 INTENT(INOUT)               ::             REGIME, &
    280. 

  280.                                                               HFX, &
    281. 

  281.                                                               QFX, &
    282. 

  282.                                                          MOL,RMOL
    283. 

  283. !m the following 5 are changed to memory size---
    284. 

  284. !
    285. 

  285.       REAL,     DIMENSION( ims:ime )                             , &
    286. 

  286.                 INTENT(INOUT)   ::                 GZ1OZ0,WSPD,BR, &
    287. 

  287.                                                         PSIM,PSIH
    288. 

  288. 

  289.       REAL,     DIMENSION( ims:ime )                             , &
    290. 

  290.                 INTENT(INOUT)   ::                            ZNT, &
    291. 

  291.                                                               ZOL, &
    292. 

  292.                                                               UST, &
    293. 

  293.                                                               CPM, &
    294. 

  294.                                                              CHS2, &
    295. 

  295.                                                              CQS2, &
    296. 

  296.                                                               CHS
    297. 

  297. 

  298.       REAL,     DIMENSION( ims:ime )                             , &
    299. 

  299.                 INTENT(INOUT)   ::                      FLHC,FLQC
    300. 

  300. 

  301.       REAL,     DIMENSION( ims:ime )                             , &
    302. 

  302.                 INTENT(INOUT)   ::                                 &
    303. 

  303.                                                               QGH
    304. 

  304. 

  305.       REAL,     DIMENSION( ims:ime )                             , &
    306. 

  306.                 INTENT(OUT)     ::                        U10,V10, &
    307. 

  307.                                                 TH2,T2,Q2,QSFC,LH
    308. 

  308. 

  309.                                     
    310. 

  310.       REAL,     INTENT(IN   )               ::   CP,G,ROVCP,R,XLV,DX
    311. 

  311. 

  312. ! MODULE-LOCAL VARIABLES, DEFINED IN SUBROUTINE SFCLAY
    313. 

  313.       REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::  dz8w1d
    314. 

  314. 

  315.       REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::      UX, &
    316. 

  316.                                                                VX, &
    317. 

  317.                                                              QV1D, &
    318. 

  318.                                                               P1D, &
    319. 

  319.                                                               T1D
    320. 

  320.  
    321. 

  321.       REAL, OPTIONAL, DIMENSION( ims:ime )                       , &
    322. 

  322.                 INTENT(OUT)     ::              ck,cka,cd,cda,ustm
    323. 

  323. 

  324.       INTEGER,  OPTIONAL,  INTENT(IN )   ::     ISFTCFLX, IZ0TLND
    325. 

  325. 

  326. ! LOCAL VARS
    327. 

  327. 

  328.       REAL,     DIMENSION( its:ite )        ::                 ZA, &
    329. 

  329.                                                         THVX,ZQKL, &
    330. 

  330.                                                            ZQKLP1, &
    331. 

  331.                                                            THX,QX, &
    332. 

  332.                                                             PSIH2, &
    333. 

  333.                                                             PSIM2, &
    334. 

  334.                                                            PSIH10, &
    335. 

  335.                                                            PSIM10, &
    336. 

  336.                                                            DENOMQ, &
    337. 

  337.                                                           DENOMQ2, &
    338. 

  338.                                                           DENOMT2, &
    339. 

  339.                                                             WSPDI, &
    340. 

  340.                                                            GZ2OZ0, &
    341. 

  341.                                                            GZ10OZ0
    342. 

  342. !
    343. 

  343.       REAL,     DIMENSION( its:ite )        ::                     &
    344. 

  344.                                                       RHOX,GOVRTH, &
    345. 

  345.                                                             TGDSA
    346. 

  346. !
    347. 

  347.       REAL,     DIMENSION( its:ite)         ::          SCR3,SCR4
    348. 

  348.       REAL,     DIMENSION( its:ite )        ::         THGB, PSFC
    349. 

  349. !
    350. 

  350.       INTEGER                               ::                 KL
    351. 

  351. 

  352.       INTEGER ::  N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10
    353. 

  353. 

  354.       REAL    ::  PL,THCON,TVCON,E1
    355. 

  355.       REAL    ::  ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10
    356. 

  356.       REAL    ::  DTG,PSIX,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2,PSIQ10
    357. 

  357.       REAL    ::  FLUXC,VSGD,Z0Q,VISC,RESTAR,CZIL,RESTAR2
    358. 

  358. !
    359. 

  359. ! .... paj ...
    360. 

  360. !
    361. 

  361.       REAL    :: zolri,zolri2,zolzz,zol0
    362. 

  362.       REAL    :: psih_stable,psim_stable,psih_unstable,psim_unstable
    363. 

  363.       REAL    :: zl2,zl10,z0t
    364. 

  364.       REAL,     DIMENSION( its:ite )        ::         pq,pq2,pq10
    365. 

  365. 

  366. 

  367. !-------------------------------------------------------------------
    368. 

  368.       KL=kte
    369. 

  369. 

  370.       DO i=its,ite
    371. 

  371. ! PSFC cb
    372. 

  372.          PSFC(I)=PSFCPA(I)/1000.
    373. 

  373.       ENDDO
    374. 

  374. !                                                      
    375. 

  375. !----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:  
    376. 

  376. !                                                            
    377. 

  377.       DO 5 I=its,ite                                   
    378. 

  378.         TGDSA(I)=TSK(I)                                    
    379. 

  379. ! PSFC cb
    380. 

  380. !        THGB(I)=TSK(I)*(100./PSFC(I))**ROVCP                
    381. 

  381.         THGB(I)=TSK(I)*(P1000mb/PSFCPA(I))**ROVCP   
    382. 

  382.     5 CONTINUE                                               
    383. 

  383. !                                                            
    384. 

  384. !-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,
    385. 

  385. !     T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1.  
    386. 

  386. !                                                                 
    387. 

  387. !     *** NOTE ***                                           
    388. 

  388. !         THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,         
    389. 

  389. !         SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE 
    390. 

  390. !         TENDENCIES.                             
    391. 

  391. !                                                           
    392. 

  392.    10 CONTINUE                                                     
    393. 

  393. 

  394. !     DO 24 I=its,ite
    395. 

  395. !        UX(I)=U1D(I)
    396. 

  396. !        VX(I)=V1D(I)
    397. 

  397. !  24 CONTINUE                                             
    398. 

  398.                                                              
    399. 

  399.    26 CONTINUE                                               
    400. 

  400.                                                    
    401. 

  401. !.....SCR3(I,K) STORE TEMPERATURE,                           
    402. 

  402. !     SCR4(I,K) STORE VIRTUAL TEMPERATURE.                                       
    403. 

  403.                                                                                  
    404. 

  404.       DO 30 I=its,ite
    405. 

  405. ! PL cb
    406. 

  406.          PL=P1D(I)/1000.
    407. 

  407.          SCR3(I)=T1D(I)                                                   
    408. 

  408. !         THCON=(100./PL)**ROVCP                                                 
    409. 

  409.          THCON=(P1000mb*0.001/PL)**ROVCP
    410. 

  410.          THX(I)=SCR3(I)*THCON                                               
    411. 

  411.          SCR4(I)=SCR3(I)                                                    
    412. 

  412.          THVX(I)=THX(I)                                                     
    413. 

  413.          QX(I)=0.                                                             
    414. 

  414.    30 CONTINUE                                                                 
    415. 

  415. !                                                                                
    416. 

  416.       DO I=its,ite
    417. 

  417.          QGH(I)=0.                                                                
    418. 

  418.          FLHC(I)=0.                                                               
    419. 

  419.          FLQC(I)=0.                                                               
    420. 

  420.          CPM(I)=CP                                                                
    421. 

  421.       ENDDO
    422. 

  422. !                                                                                
    423. 

  423. !     IF(IDRY.EQ.1)GOTO 80                                                   
    424. 

  424.       DO 50 I=its,ite
    425. 

  425.          QX(I)=QV1D(I)                                                    
    426. 

  426.          TVCON=(1.+EP1*QX(I))                                      
    427. 

  427.          THVX(I)=THX(I)*TVCON                                               
    428. 

  428.          SCR4(I)=SCR3(I)*TVCON                                              
    429. 

  429.    50 CONTINUE                                                                 
    430. 

  430. !                                                                                
    431. 

  431.       DO 60 I=its,ite
    432. 

  432.         E1=SVP1*EXP(SVP2*(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))                       
    433. 

  433. !  for land points QSFC can come from previous time step
    434. 

  434.         if(xland(i).gt.1.5.or.qsfc(i).le.0.0)QSFC(I)=EP2*E1/(PSFC(I)-E1)                                                 
    435. 

  435. ! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP CONSISTENT WITH MYJSFC CHANGE
    436. 

  436. ! Q2SAT = QGH IN LSM
    437. 

  437.         E1=SVP1*EXP(SVP2*(T1D(I)-SVPT0)/(T1D(I)-SVP3))                       
    438. 

  438.         PL=P1D(I)/1000.
    439. 

  439.         QGH(I)=EP2*E1/(PL-E1)                                                 
    440. 

  440.         CPM(I)=CP*(1.+0.8*QX(I))                                   
    441. 

  441.    60 CONTINUE                                                                   
    442. 

  442.    80 CONTINUE
    443. 

  443.                                                                                  
    444. 

  444. !-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND             
    445. 

  445. !     LEVEL, AND THE LAYER THICKNESSES.                                          
    446. 

  446.                                                                                  
    447. 

  447.       DO 90 I=its,ite
    448. 

  448.         ZQKLP1(I)=0.
    449. 

  449.         RHOX(I)=PSFC(I)*1000./(R*SCR4(I))                                       
    450. 

  450.    90 CONTINUE                                                                   
    451. 

  451. !                                                                                
    452. 

  452.       DO 110 I=its,ite                                                   
    453. 

  453.            ZQKL(I)=dz8w1d(I)+ZQKLP1(I)
    454. 

  454.   110 CONTINUE                                                                 
    455. 

  455. !                                                                                
    456. 

  456.       DO 120 I=its,ite
    457. 

  457.          ZA(I)=0.5*(ZQKL(I)+ZQKLP1(I))                                        
    458. 

  458.   120 CONTINUE                                                                 
    459. 

  459. !                                                                                
    460. 

  460.       DO 160 I=its,ite
    461. 

  461.         GOVRTH(I)=G/THX(I)                                                    
    462. 

  462.   160 CONTINUE                                                                   
    463. 

  463.                                                                                  
    464. 

  464. !-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO               
    465. 

  465. !     AKB(1976), EQ(12).                                                         
    466. 

  466.                    
    467. 

  467.       DO 260 I=its,ite
    468. 

  468.         GZ1OZ0(I)=ALOG((ZA(I)+ZNT(I))/ZNT(I))   ! log((z+z0)/z0)                                     
    469. 

  469.         GZ2OZ0(I)=ALOG((2.+ZNT(I))/ZNT(I))      ! log((2+z0)/z0)                           
    470. 

  470.         GZ10OZ0(I)=ALOG((10.+ZNT(I))/ZNT(I))    ! log((10+z0)z0)                    
    471. 

  471.         IF((XLAND(I)-1.5).GE.0)THEN                                            
    472. 

  472.           ZL=ZNT(I)                                                            
    473. 

  473.         ELSE                                                                     
    474. 

  474.           ZL=0.01                                                                
    475. 

  475.         ENDIF                                                                    
    476. 

  476.         WSPD(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))                        
    477. 

  477. 

  478.         TSKV=THGB(I)*(1.+EP1*QSFC(I))                     
    479. 

  479.         DTHVDZ=(THVX(I)-TSKV)                                                 
    480. 

  480. !  Convective velocity scale Vc and subgrid-scale velocity Vsg
    481. 

  481. !  following Beljaars (1995, QJRMS) and Mahrt and Sun (1995, MWR)
    482. 

  482. !                                ... HONG Aug. 2001
    483. 

  483. !
    484. 

  484. !       VCONV = 0.25*sqrt(g/tskv*pblh(i)*dthvm)
    485. 

  485. !      Use Beljaars over land, old MM5 (Wyngaard) formula over water
    486. 

  486.         if (xland(i).lt.1.5) then
    487. 

  487.         fluxc = max(hfx(i)/rhox(i)/cp                    &
    488. 

  488.               + ep1*tskv*qfx(i)/rhox(i),0.)
    489. 

  489.         VCONV = vconvc*(g/tgdsa(i)*pblh(i)*fluxc)**.33
    490. 

  490.         else
    491. 

  491.         IF(-DTHVDZ.GE.0)THEN
    492. 

  492.           DTHVM=-DTHVDZ
    493. 

  493.         ELSE
    494. 

  494.           DTHVM=0.
    495. 

  495.         ENDIF
    496. 

  496.         VCONV = 2.*SQRT(DTHVM)
    497. 

  497.         endif
    498. 

  498. ! Mahrt and Sun low-res correction
    499. 

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

  500.         WSPD(I)=SQRT(WSPD(I)*WSPD(I)+VCONV*VCONV+vsgd*vsgd)
    501. 

  501.         WSPD(I)=AMAX1(WSPD(I),0.1)
    502. 

  502.         BR(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD(I)*WSPD(I))                        
    503. 

  503. !  IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2
    504. 

  504.         IF(MOL(I).LT.0.)BR(I)=AMIN1(BR(I),0.0)
    505. 

  505. !jdf
    506. 

  506.         RMOL(I)=-GOVRTH(I)*DTHVDZ*ZA(I)*KARMAN
    507. 

  507. !jdf
    508. 

  508. 

  509.   260 CONTINUE                                                                   
    510. 

  510. 

  511. !                                                                                
    512. 

  512. !-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:            
    513. 

  513. !                                                                                
    514. 

  514. !                                                                                
    515. 

  515. !     THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)           
    516. 

  516. !     AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).                              
    517. 

  517. !                                                                                
    518. 

  518. !     CRITERIA FOR THE CLASSES ARE AS FOLLOWS:                                   
    519. 

  519. !                                                                                
    520. 

  520. !        1. BR .GE. 0.0;                                                         
    521. 

  521. !               REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),               
    522. 

  522. !                                                                                
    523. 

  523. !        3. BR .EQ. 0.0                                                          
    524. 

  524. !               REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),              
    525. 

  525. !                                                                                
    526. 

  526. !        4. BR .LT. 0.0                                                          
    527. 

  527. !               REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).                
    528. 

  528. !                                                                                
    529. 

  529. !CCCCC                                                                           
    530. 

  530. 

  531.       DO 320 I=its,ite
    532. 

  532. !                                                                           
    533. 

  533.       if (br(I).gt.0) then
    534. 

  534.         if (br(I).gt.250.0) then
    535. 

  535.         zol(I)=zolri(250.0,ZA(I),ZNT(I))
    536. 

  536.         else
    537. 

  537.         zol(I)=zolri(br(I),ZA(I),ZNT(I))
    538. 

  538.         endif
    539. 

  539.       endif
    540. 

  540. !
    541. 

  541.       if (br(I).lt.0) then
    542. 

  542.        IF(UST(I).LT.0.001)THEN
    543. 

  543.           ZOL(I)=BR(I)*GZ1OZ0(I)
    544. 

  544.         ELSE
    545. 

  545.         if (br(I).lt.-250.0) then
    546. 

  546.         zol(I)=zolri(-250.0,ZA(I),ZNT(I))
    547. 

  547.         else
    548. 

  548.         zol(I)=zolri(br(I),ZA(I),ZNT(I))
    549. 

  549.         endif
    550. 

  550.        ENDIF
    551. 

  551.       endif
    552. 

  552. !
    553. 

  553. ! ... paj: compute integrated similarity functions.
    554. 

  554. !
    555. 

  555.         zolzz=zol(I)*(za(I)+znt(I))/za(I) ! (z+z0/L
    556. 

  556.         zol10=zol(I)*(10.+znt(I))/za(I)   ! (10+z0)/L
    557. 

  557.         zol2=zol(I)*(2.+znt(I))/za(I)     ! (2+z0)/L
    558. 

  558.         zol0=zol(I)*znt(I)/za(I)          ! z0/L
    559. 

  559.         ZL2=(2.)/ZA(I)*ZOL(I)             ! 2/L      
    560. 

  560.         ZL10=(10.)/ZA(I)*ZOL(I)           ! 10/L
    561. 

  561. 

  562.         IF((XLAND(I)-1.5).LT.0.)THEN
    563. 

  563.         ZL=(0.01)/ZA(I)*ZOL(I)   ! (0.01)/L     
    564. 

  564.         ELSE
    565. 

  565.         ZL=ZOL0                     ! z0/L
    566. 

  566.         ENDIF
    567. 

  567. 

  568.         IF(BR(I).LT.0.)GOTO 310  ! go to unstable regime (class 4)
    569. 

  569.         IF(BR(I).EQ.0.)GOTO 280  ! go to neutral regime (class 3)
    570. 

  570. !                                                                                
    571. 

  571. !-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:                                    
    572. 

  572. !
    573. 

  573.         REGIME(I)=1.
    574. 

  574. !
    575. 

  575. ! ... paj: psim and psih. Follows Cheng and Brutsaert 2005 (CB05).
    576. 

  576. !
    577. 

  577.         psim(I)=psim_stable(zolzz)-psim_stable(zol0)
    578. 

  578.         psih(I)=psih_stable(zolzz)-psih_stable(zol0)
    579. 

  579. !
    580. 

  580.         psim10(I)=psim_stable(zol10)-psim_stable(zol0)
    581. 

  581.         psih10(I)=psih_stable(zol10)-psih_stable(zol0)
    582. 

  582. !
    583. 

  583.         psim2(I)=psim_stable(zol2)-psim_stable(zol0)
    584. 

  584.         psih2(I)=psih_stable(zol2)-psih_stable(zol0)
    585. 

  585. !
    586. 

  586. ! ... paj: preparations to compute PSIQ. Follows CB05+Carlson Boland JAM 1978.
    587. 

  587. !
    588. 

  588.         pq(I)=psih_stable(zol(I))-psih_stable(zl)
    589. 

  589.         pq2(I)=psih_stable(zl2)-psih_stable(zl)
    590. 

  590.         pq10(I)=psih_stable(zl10)-psih_stable(zl)
    591. 

  591. !
    592. 

  592. !       1.0 over Monin-Obukhov length
    593. 

  593.         RMOL(I)=ZOL(I)/ZA(I) 
    594. 

  594. !                                                                                
    595. 

  595. 

  596.         GOTO 320                                                                 
    597. 

  597. !                                                                                
    598. 

  598. !-----CLASS 3; FORCED CONVECTION:                                                
    599. 

  599. !                                                                                
    600. 

  600.   280   REGIME(I)=3.                                                           
    601. 

  601.         PSIM(I)=0.0                                                              
    602. 

  602.         PSIH(I)=PSIM(I)                                                          
    603. 

  603.         PSIM10(I)=0.                                                   
    604. 

  604.         PSIH10(I)=PSIM10(I)                                           
    605. 

  605.         PSIM2(I)=0.                                                  
    606. 

  606.         PSIH2(I)=PSIM2(I)                                           
    607. 

  607. !
    608. 

  608. ! paj: preparations to compute PSIQ.
    609. 

  609. !
    610. 

  610.         pq(I)=PSIH(I)
    611. 

  611.         pq2(I)=PSIH2(I)
    612. 

  612.         pq10(I)=0.
    613. 

  613. !
    614. 

  614.         ZOL(I)=0.                                             
    615. 

  615.         RMOL(I) = ZOL(I)/ZA(I)  
    616. 

  616. 

  617.         GOTO 320                                                                 
    618. 

  618. !                                                                                
    619. 

  619. !-----CLASS 4; FREE CONVECTION:                                                  
    620. 

  620. !                                                                                
    621. 

  621.   310   CONTINUE                                                                 
    622. 

  622.         REGIME(I)=4.                                                           
    623. 

  623. !
    624. 

  624. ! ... paj: PSIM and PSIH ...
    625. 

  625. !
    626. 

  626.         psim(I)=psim_unstable(zolzz)-psim_unstable(zol0)
    627. 

  627.         psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
    628. 

  628. !
    629. 

  629.         psim10(I)=psim_unstable(zol10)-psim_unstable(zol0)
    630. 

  630.         psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
    631. 

  631. !
    632. 

  632.         psim2(I)=psim_unstable(zol2)-psim_unstable(zol0)
    633. 

  633.         psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
    634. 

  634. !
    635. 

  635. ! ... paj: preparations to compute PSIQ 
    636. 

  636. !
    637. 

  637.         pq(I)=psih_unstable(zol(I))-psih_unstable(zl)
    638. 

  638.         pq2(I)=psih_unstable(zl2)-psih_unstable(zl)
    639. 

  639.         pq10(I)=psih_unstable(zl10)-psih_unstable(zl)
    640. 

  640. !
    641. 

  641. !---LIMIOT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS            
    642. 

  642. !---  THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL                 
    643. 

  643.         PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
    644. 

  644.         PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
    645. 

  645.         PSIH2(I)=AMIN1(PSIH2(I),0.9*GZ2OZ0(I))
    646. 

  646.         PSIM10(I)=AMIN1(PSIM10(I),0.9*GZ10OZ0(I))
    647. 

  647. !
    648. 

  648. ! AHW: mods to compute ck, cd
    649. 

  649.         PSIH10(I)=AMIN1(PSIH10(I),0.9*GZ10OZ0(I))
    650. 

  650. 

  651.         RMOL(I) = ZOL(I)/ZA(I)  
    652. 

  652. 

  653.   320 CONTINUE                                                                   
    654. 

  654. !                                                                                
    655. 

  655. !-----COMPUTE THE FRICTIONAL VELOCITY:                                           
    656. 

  656. !     ZA(1982) EQS(2.60),(2.61).                                                 
    657. 

  657. !                                                                                
    658. 

  658.       DO 330 I=its,ite
    659. 

  659.         DTG=THX(I)-THGB(I)                                                   
    660. 

  660.         PSIX=GZ1OZ0(I)-PSIM(I)                                                   
    661. 

  661.         PSIX10=GZ10OZ0(I)-PSIM10(I)
    662. 

  662. 

  663. !     LOWER LIMIT ADDED TO PREVENT LARGE FLHC IN SOIL MODEL
    664. 

  664. !     ACTIVATES IN UNSTABLE CONDITIONS WITH THIN LAYERS OR HIGH Z0
    665. 

  665. !       PSIT=AMAX1(GZ1OZ0(I)-PSIH(I),2.)
    666. 

  666.        PSIT=GZ1OZ0(I)-PSIH(I)
    667. 

  667.        PSIT2=GZ2OZ0(I)-PSIH2(I)
    668. 

  668. !
    669. 

  669.         IF((XLAND(I)-1.5).GE.0)THEN                                            
    670. 

  670.           ZL=ZNT(I)                                                            
    671. 

  671.         ELSE                                                                     
    672. 

  672.           ZL=0.01                                                                
    673. 

  673.         ENDIF                                                                    
    674. 

  674. !
    675. 

  675.         PSIQ=ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-pq(I)
    676. 

  676.         PSIQ2=ALOG(KARMAN*UST(I)*2./XKA+2./ZL)-pq2(I)
    677. 

  677. 

  678. ! AHW: mods to compute ck, cd
    679. 

  679.         PSIQ10=ALOG(KARMAN*UST(I)*10./XKA+10./ZL)-pq10(I)
    680. 

  680.         IF ( PRESENT(ISFTCFLX) ) THEN
    681. 

  681.            IF ( ISFTCFLX.EQ.1 .AND. (XLAND(I)-1.5).GE.0. ) THEN
    682. 

  682. ! v3.1
    683. 

  683. !             Z0Q = 1.e-4 + 1.e-3*(MAX(0.,UST(I)-1.))**2
    684. 

  684. ! hfip1
    685. 

  685. !             Z0Q = 0.62*2.0E-5/UST(I) + 1.E-3*(MAX(0.,UST(I)-1.5))**2
    686. 

  686. ! v3.2
    687. 

  687.               Z0Q = 1.e-4
    688. 

  688. !
    689. 

  689. ! ... paj: recompute psih for z0q
    690. 

  690. !
    691. 

  691.            zolzz=zol(I)*(za(I)+z0q)/za(I)    ! (z+z0q)/L
    692. 

  692.            zol10=zol(I)*(10.+z0q)/za(I)   ! (10+z0q)/L
    693. 

  693.            zol2=zol(I)*(2.+z0q)/za(I)     ! (2+z0q)/L
    694. 

  694.            zol0=zol(I)*z0q/za(I)          ! z0q/L
    695. 

  695. !
    696. 

  696.               if (zol(I).gt.0.) then
    697. 

  697.               psih(I)=psih_stable(zolzz)-psih_stable(zol0)
    698. 

  698.               psih10(I)=psih_stable(zol10)-psih_stable(zol0)
    699. 

  699.               psih2(I)=psih_stable(zol2)-psih_stable(zol0)
    700. 

  700.               else
    701. 

  701.                 if (zol(I).eq.0) then
    702. 

  702.                 psih(I)=0.
    703. 

  703.                 psih10(I)=0.
    704. 

  704.                 psih2(I)=0.
    705. 

  705.                 else
    706. 

  706.                 psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
    707. 

  707.                 psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
    708. 

  708.                 psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
    709. 

  709.                 endif
    710. 

  710.               endif
    711. 

  711. !
    712. 

  712.               PSIQ=ALOG((ZA(I)+z0q)/Z0Q)-PSIH(I)
    713. 

  713.               PSIT=PSIQ
    714. 

  714.               PSIQ2=ALOG((2.+z0q)/Z0Q)-PSIH2(I)
    715. 

  715.               PSIQ10=ALOG((10.+z0q)/Z0Q)-PSIH10(I)
    716. 

  716.               PSIT2=PSIQ2
    717. 

  717.            ENDIF
    718. 

  718.           IF ( ISFTCFLX.EQ.2 .AND. (XLAND(I)-1.5).GE.0. ) THEN
    719. 

  719. ! AHW: Garratt formula: Calculate roughness Reynolds number
    720. 

  720. !        Kinematic viscosity of air (linear approc to
    721. 

  721. !                 temp dependence at sea levle)
    722. 

  722.               VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5
    723. 

  723. !!            VISC=1.5E-5
    724. 

  724.               RESTAR=UST(I)*ZNT(I)/VISC
    725. 

  725.               RESTAR2=2.48*SQRT(SQRT(RESTAR))-2.
    726. 

  726. !
    727. 

  727. ! ... paj: compute psih for z0t for temperature ...
    728. 

  728. !
    729. 

  729.               z0t=znt(I)/exp(RESTAR2)
    730. 

  730. !
    731. 

  731.            zolzz=zol(I)*(za(I)+z0t)/za(I)    ! (z+z0t)/L
    732. 

  732.            zol10=zol(I)*(10.+z0t)/za(I)   ! (10+z0t)/L
    733. 

  733.            zol2=zol(I)*(2.+z0t)/za(I)     ! (2+z0t)/L
    734. 

  734.            zol0=zol(I)*z0t/za(I)          ! z0t/L
    735. 

  735. !
    736. 

  736.               if (zol(I).gt.0.) then
    737. 

  737.               psih(I)=psih_stable(zolzz)-psih_stable(zol0)
    738. 

  738.               psih10(I)=psih_stable(zol10)-psih_stable(zol0)
    739. 

  739.               psih2(I)=psih_stable(zol2)-psih_stable(zol0)
    740. 

  740.               else
    741. 

  741.                 if (zol(I).eq.0) then
    742. 

  742.                 psih(I)=0.
    743. 

  743.                 psih10(I)=0.
    744. 

  744.                 psih2(I)=0.
    745. 

  745.                 else
    746. 

  746.                 psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
    747. 

  747.                 psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
    748. 

  748.                 psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
    749. 

  749.                 endif
    750. 

  750.               endif
    751. 

  751. !
    752. 

  752. !              PSIT=GZ1OZ0(I)-PSIH(I)+RESTAR2
    753. 

  753. !              PSIT2=GZ2OZ0(I)-PSIH2(I)+RESTAR2
    754. 

  754.               PSIT=ALOG((ZA(I)+z0t)/Z0t)-PSIH(I)
    755. 

  755.               PSIT2=ALOG((2.+z0t)/Z0t)-PSIH2(I)
    756. 

  756. !
    757. 

  757.               z0t=znt(I)/exp(2.28*SQRT(SQRT(RESTAR))-2.)
    758. 

  758. !
    759. 

  759.            zolzz=zol(I)*(za(I)+z0t)/za(I)    ! (z+z0t)/L
    760. 

  760.            zol10=zol(I)*(10.+z0t)/za(I)   ! (10+z0t)/L
    761. 

  761.            zol2=zol(I)*(2.+z0t)/za(I)     ! (2+z0t)/L
    762. 

  762.            zol0=zol(I)*z0t/za(I)          ! z0t/L
    763. 

  763. !
    764. 

  764.               if (zol(I).gt.0.) then
    765. 

  765.               psih(I)=psih_stable(zolzz)-psih_stable(zol0)
    766. 

  766.               psih10(I)=psih_stable(zol10)-psih_stable(zol0)
    767. 

  767.               psih2(I)=psih_stable(zol2)-psih_stable(zol0)
    768. 

  768.               else
    769. 

  769.                 if (zol(I).eq.0) then
    770. 

  770.                 psih(I)=0.
    771. 

  771.                 psih10(I)=0.
    772. 

  772.                 psih2(I)=0.
    773. 

  773.                 else
    774. 

  774.                 psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
    775. 

  775.                 psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
    776. 

  776.                 psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
    777. 

  777.                 endif
    778. 

  778.               endif
    779. 

  779. !
    780. 

  780.               PSIQ=ALOG((ZA(I)+z0t)/Z0t)-PSIH(I)
    781. 

  781.               PSIQ2=ALOG((2.+z0t)/Z0t)-PSIH2(I)
    782. 

  782.               PSIQ10=ALOG((10.+z0t)/Z0t)-PSIH10(I)
    783. 

  783. !              PSIQ=GZ1OZ0(I)-PSIH(I)+2.28*SQRT(SQRT(RESTAR))-2.
    784. 

  784. !              PSIQ2=GZ2OZ0(I)-PSIH2(I)+2.28*SQRT(SQRT(RESTAR))-2.
    785. 

  785. !              PSIQ10=GZ10OZ0(I)-PSIH(I)+2.28*SQRT(SQRT(RESTAR))-2.
    786. 

  786.            ENDIF
    787. 

  787.         ENDIF
    788. 

  788.         IF(PRESENT(ck) .and. PRESENT(cd) .and. PRESENT(cka) .and. PRESENT(cda)) THEN
    789. 

  789.            Ck(I)=(karman/psix10)*(karman/psiq10)
    790. 

  790.            Cd(I)=(karman/psix10)*(karman/psix10)
    791. 

  791.            Cka(I)=(karman/psix)*(karman/psiq)
    792. 

  792.            Cda(I)=(karman/psix)*(karman/psix)
    793. 

  793.         ENDIF
    794. 

  794.         IF ( PRESENT(IZ0TLND) ) THEN
    795. 

  795.            IF ( IZ0TLND.EQ.1 .AND. (XLAND(I)-1.5).LE.0. ) THEN
    796. 

  796.               ZL=ZNT(I)
    797. 

  797. !             CZIL RELATED CHANGES FOR LAND
    798. 

  798.               VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5
    799. 

  799.               RESTAR=UST(I)*ZL/VISC
    800. 

  800. !             Modify CZIL according to Chen & Zhang, 2009
    801. 

  801. 

  802.               CZIL = 10.0 ** ( -0.40 * ( ZL / 0.07 ) )
    803. 

  803. !
    804. 

  804. ! ... paj: compute phish for z0t over land
    805. 

  805. !
    806. 

  806.               z0t=znt(I)/exp(CZIL*KARMAN*SQRT(RESTAR))
    807. 

  807. !
    808. 

  808.            zolzz=zol(I)*(za(I)+z0t)/za(I)    ! (z+z0t)/L
    809. 

  809.            zol10=zol(I)*(10.+z0t)/za(I)   ! (10+z0t)/L
    810. 

  810.            zol2=zol(I)*(2.+z0t)/za(I)     ! (2+z0t)/L
    811. 

  811.            zol0=zol(I)*z0t/za(I)          ! z0t/L
    812. 

  812. !
    813. 

  813.               if (zol(I).gt.0.) then
    814. 

  814.               psih(I)=psih_stable(zolzz)-psih_stable(zol0)
    815. 

  815.               psih10(I)=psih_stable(zol10)-psih_stable(zol0)
    816. 

  816.               psih2(I)=psih_stable(zol2)-psih_stable(zol0)
    817. 

  817.               else
    818. 

  818.                 if (zol(I).eq.0) then
    819. 

  819.                 psih(I)=0.
    820. 

  820.                 psih10(I)=0.
    821. 

  821.                 psih2(I)=0.
    822. 

  822.                 else
    823. 

  823.                 psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
    824. 

  824.                 psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
    825. 

  825.                 psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
    826. 

  826.                 endif
    827. 

  827.               endif
    828. 

  828. !
    829. 

  829.               PSIQ=ALOG((ZA(I)+z0t)/Z0t)-PSIH(I)
    830. 

  830.               PSIQ2=ALOG((2.+z0t)/Z0t)-PSIH2(I)
    831. 

  831.               PSIT=PSIQ
    832. 

  832.               PSIT2=PSIQ2
    833. 

  833. !
    834. 

  834. !              PSIT=GZ1OZ0(I)-PSIH(I)+CZIL*KARMAN*SQRT(RESTAR)
    835. 

  835. !              PSIQ=GZ1OZ0(I)-PSIH(I)+CZIL*KARMAN*SQRT(RESTAR)
    836. 

  836. !              PSIT2=GZ2OZ0(I)-PSIH2(I)+CZIL*KARMAN*SQRT(RESTAR)
    837. 

  837. !              PSIQ2=GZ2OZ0(I)-PSIH2(I)+CZIL*KARMAN*SQRT(RESTAR)
    838. 

  838. 

  839.            ENDIF
    840. 

  840.         ENDIF
    841. 

  841. ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE 
    842. 

  842.         UST(I)=0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX                                             
    843. 

  843. ! TKE coupling: compute ust without vconv for use in tke scheme
    844. 

  844.         WSPDI(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))
    845. 

  845.         IF ( PRESENT(USTM) ) THEN
    846. 

  846.         USTM(I)=0.5*USTM(I)+0.5*KARMAN*WSPDI(I)/PSIX
    847. 

  847.         ENDIF
    848. 

  848. 

  849.         U10(I)=UX(I)*PSIX10/PSIX                                    
    850. 

  850.         V10(I)=VX(I)*PSIX10/PSIX                                   
    851. 

  851.         TH2(I)=THGB(I)+DTG*PSIT2/PSIT                                
    852. 

  852.         Q2(I)=QSFC(I)+(QX(I)-QSFC(I))*PSIQ2/PSIQ                   
    853. 

  853.         T2(I) = TH2(I)*(PSFCPA(I)/P1000mb)**ROVCP                     
    854. 

  854. !                                                                                
    855. 

  855.         IF((XLAND(I)-1.5).LT.0.)THEN                                            
    856. 

  856.           UST(I)=AMAX1(UST(I),0.001)
    857. 

  857.         ENDIF                                                                    
    858. 

  858.         MOL(I)=KARMAN*DTG/PSIT/PRT                              
    859. 

  859.         DENOMQ(I)=PSIQ
    860. 

  860.         DENOMQ2(I)=PSIQ2
    861. 

  861.         DENOMT2(I)=PSIT2
    862. 

  862.   330 CONTINUE                                                                   
    863. 

  863. !                                                                                
    864. 

  864.   335 CONTINUE                                                                   
    865. 

  865.                                                                                   
    866. 

  866. !-----COMPUTE THE SURFACE SENSIBLE AND LATENT HEAT FLUXES:                       
    867. 

  867.                                                                                  
    868. 

  868.       DO i=its,ite
    869. 

  869.         QFX(i)=0.                                                              
    870. 

  870.         HFX(i)=0.                                                              
    871. 

  871.       ENDDO
    872. 

  872. 

  873.       IF (ISFFLX.EQ.0) GOTO 410                                                
    874. 

  874.                                                                                  
    875. 

  875. !-----OVER WATER, ALTER ROUGHNESS LENGTH (ZNT) ACCORDING TO WIND (UST).          
    876. 

  876.                                                                                  
    877. 

  877.       DO 360 I=its,ite
    878. 

  878.         IF((XLAND(I)-1.5).GE.0)THEN                                            
    879. 

  879.           ZNT(I)=CZO*UST(I)*UST(I)/G+OZO                                   
    880. 

  880. ! AHW: change roughness length, and hence the drag coefficients Ck and Cd
    881. 

  881.           IF ( PRESENT(ISFTCFLX) ) THEN
    882. 

  882.              IF ( ISFTCFLX.NE.0 ) THEN
    883. 

  883. !               ZNT(I)=10.*exp(-9.*UST(I)**(-.3333))
    884. 

  884.                 ZNT(I)=10.*exp(-9.5*UST(I)**(-.3333))
    885. 

  885.                 ZNT(I)=ZNT(I) + 0.11*1.5E-5/AMAX1(UST(I),0.01)
    886. 

  886.                 ZNT(I)=MIN(ZNT(I),2.85e-3)
    887. 

  887.                 ZNT(I)=MAX(ZNT(I),1.27e-7)
    888. 

  888.              ENDIF
    889. 

  889.           ENDIF
    890. 

  890.           ZL = ZNT(I)
    891. 

  891.         ELSE
    892. 

  892.           ZL = 0.01
    893. 

  893.         ENDIF                                                                    
    894. 

  894.         FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/DENOMQ(I)
    895. 

  895. !       FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/(   &
    896. 

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

  897.         DTTHX=ABS(THX(I)-THGB(I))                                            
    898. 

  898.         IF(DTTHX.GT.1.E-5)THEN                                                   
    899. 

  899.           FLHC(I)=CPM(I)*RHOX(I)*UST(I)*MOL(I)/(THX(I)-THGB(I))          
    900. 

  900. !         write(*,1001)FLHC(I),CPM(I),RHOX(I),UST(I),MOL(I),THX(I),THGB(I),I
    901. 

  901.  1001   format(f8.5,2x,f12.7,2x,f12.10,2x,f12.10,2x,f13.10,2x,f12.8,f12.8,2x,i3)
    902. 

  902.         ELSE                                                                     
    903. 

  903.           FLHC(I)=0.                                                             
    904. 

  904.         ENDIF                                                                    
    905. 

  905.   360 CONTINUE                                                                   
    906. 

  906. 

  907. !                                                                                
    908. 

  908. !-----COMPUTE SURFACE MOIST FLUX:                                                
    909. 

  909. !                                                                                
    910. 

  910. !     IF(IDRY.EQ.1)GOTO 390                                                
    911. 

  911. !                                                                                
    912. 

  912.       DO 370 I=its,ite
    913. 

  913.         QFX(I)=FLQC(I)*(QSFC(I)-QX(I))                                     
    914. 

  914.         QFX(I)=AMAX1(QFX(I),0.)                                            
    915. 

  915.         LH(I)=XLV*QFX(I)
    916. 

  916.   370 CONTINUE                                                                 
    917. 

  917.                                                                                 
    918. 

  918. !-----COMPUTE SURFACE HEAT FLUX:                                                 
    919. 

  919. !                                                                                
    920. 

  920.   390 CONTINUE                                                                 
    921. 

  921.       DO 400 I=its,ite
    922. 

  922.         IF(XLAND(I)-1.5.GT.0.)THEN                                           
    923. 

  923.           HFX(I)=FLHC(I)*(THGB(I)-THX(I)) 
    924. 

  924.           IF ( PRESENT(ISFTCFLX) ) THEN
    925. 

  925.              IF ( ISFTCFLX.NE.0 ) THEN
    926. 

  926. ! AHW: add dissipative heating term
    927. 

  927.                 HFX(I)=HFX(I)+RHOX(I)*USTM(I)*USTM(I)*WSPDI(I)
    928. 

  928.              ENDIF
    929. 

  929.           ENDIF 
    930. 

  930.         ELSEIF(XLAND(I)-1.5.LT.0.)THEN                                       
    931. 

  931.           HFX(I)=FLHC(I)*(THGB(I)-THX(I))                                
    932. 

  932.           HFX(I)=AMAX1(HFX(I),-250.)                                       
    933. 

  933.         ENDIF                                                                  
    934. 

  934.   400 CONTINUE                                                                 
    935. 

  935.          
    936. 

  936.       DO I=its,ite
    937. 

  937.          IF((XLAND(I)-1.5).GE.0)THEN
    938. 

  938.            ZL=ZNT(I)
    939. 

  939.          ELSE
    940. 

  940.            ZL=0.01
    941. 

  941.          ENDIF
    942. 

  942. !v3.1.1
    943. 

  943. !         CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) &
    944. 

  944. !                /XKA+ZA(I)/ZL)-PSIH(I))
    945. 

  945.          CHS(I)=UST(I)*KARMAN/DENOMQ(I)
    946. 

  946. !        GZ2OZ0(I)=ALOG(2./ZNT(I))
    947. 

  947. !        PSIM2(I)=-10.*GZ2OZ0(I)
    948. 

  948. !        PSIM2(I)=AMAX1(PSIM2(I),-10.)
    949. 

  949. !        PSIH2(I)=PSIM2(I)
    950. 

  950. ! v3.1.1
    951. 

  951. !         CQS2(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*2.0  &
    952. 

  952. !               /XKA+2.0/ZL)-PSIH2(I))
    953. 

  953. !         CHS2(I)=UST(I)*KARMAN/(GZ2OZ0(I)-PSIH2(I))
    954. 

  954.          CQS2(I)=UST(I)*KARMAN/DENOMQ2(I)
    955. 

  955.          CHS2(I)=UST(I)*KARMAN/DENOMT2(I)
    956. 

  956.       ENDDO
    957. 

  957.                                                                         
    958. 

  958.   410 CONTINUE                                                                   
    959. 

  959. !jdf
    960. 

  960. !     DO I=its,ite
    961. 

  961. !       IF(UST(I).GE.0.1) THEN
    962. 

  962. !         RMOL(I)=RMOL(I)*(-FLHC(I))/(UST(I)*UST(I)*UST(I))
    963. 

  963. !       ELSE
    964. 

  964. !         RMOL(I)=RMOL(I)*(-FLHC(I))/(0.1*0.1*0.1)
    965. 

  965. !       ENDIF
    966. 

  966. !     ENDDO
    967. 

  967. !jdf
    968. 

  968. 

  969. !                                                                                
    970. 

  970.    END SUBROUTINE SFCLAYREV1D
    971. 

  971. 

  972. !====================================================================
    973. 

  973.    SUBROUTINE sfclayrevinit( allowed_to_read )         
    974. 

  974. 

  975. !   LOGICAL , INTENT(IN)      ::      allowed_to_read
    976. 

  976. !   INTEGER                   ::      N
    977. 

  977. !   REAL                      ::      ZOLN,X,Y
    978. 

  978. 

  979. !   DO N=0,1000
    980. 

  980. !      ZOLN=-FLOAT(N)*0.1
    981. 

  981. !      X=(1-16.*ZOLN)**0.25
    982. 

  982. !      PSIMTB(N)=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))- &
    983. 

  983. !                2.*ATAN(X)+2.*ATAN(1.)
    984. 

  984. !      Y=(1-16*ZOLN)**0.5
    985. 

  985. !      PSIHTB(N)=2*ALOG(0.5*(1+Y))
    986. 

  986. !   ENDDO
    987. 

  987. 

  988.    END SUBROUTINE sfclayrevinit
    989. 

  989. 

  990. !-------------------------------------------------------------------          
    991. 

  991. 

  992. END MODULE module_sf_sfclayrev
    993. 

  993. 

  994. !
    995. 

  995. ! ----------------------------------------------------------
    996. 

  996. !
    997. 

  997.       function zolri(ri,z,z0)
    998. 

  998. !
    999. 

  999.       real right,left,midpoint
    1000. 

  1000. !
    1001. 

  1001.       if (ri.lt.0.)then
    1002. 

  1002.       left=-1200.
    1003. 

  1003.       right=0.
    1004. 

  1004.       else
    1005. 

  1005.       left=0.
    1006. 

  1006.       right=20000.
    1007. 

  1007.       endif
    1008. 

  1008. !
    1009. 

  1009.       Do While (abs(right - left) > 0.01)
    1010. 

  1010.       midpoint=(right+left)/2.
    1011. 

  1011.       a=zolri2(left,ri,z,z0)
    1012. 

  1012.       b=zolri2(midpoint,ri,z,z0)
    1013. 

  1013.       c=zolri2(right,ri,z,z0)
    1014. 

  1014. !
    1015. 

  1015.       if ((a * b) < 0) then
    1016. 

  1016.       right = midpoint
    1017. 

  1017.       else
    1018. 

  1018.             if ((c * b) < 0) then
    1019. 

  1019.             left = midpoint
    1020. 

  1020.             else
    1021. 

  1021.             goto 11
    1022. 

  1022.             endif
    1023. 

  1023.       endif
    1024. 

  1024. !
    1025. 

  1025.       enddo
    1026. 

  1026.  11   continue
    1027. 

  1027. !
    1028. 

  1028.       zolri=midpoint
    1029. 

  1029. 

  1030.       return
    1031. 

  1031.       end
    1032. 

  1032. 

  1033. !
    1034. 

  1034. ! -----------------------------------------------------------------------
    1035. 

  1035. !
    1036. 

  1036.       function zolri2(zol2,ri2,z,z0)
    1037. 

  1037. !
    1038. 

  1038.       zol20=zol2*z0/z ! z0/L
    1039. 

  1039.       zol3=zol2+zol20 ! (z+z0)/L
    1040. 

  1040. !
    1041. 

  1041.       if (ri2.lt.0) then
    1042. 

  1042.       psix2=log((z+z0)/z0)-(psim_unstable(zol3)-psim_unstable(zol20))
    1043. 

  1043.       psih2=log((z+z0)/z0)-(psih_unstable(zol3)-psih_unstable(zol20))
    1044. 

  1044.       else
    1045. 

  1045.       psix2=log((z+z0)/z0)-(psim_stable(zol3)-psim_stable(zol20))
    1046. 

  1046.       psih2=log((z+z0)/z0)-(psih_stable(zol3)-psih_stable(zol20))
    1047. 

  1047.       endif
    1048. 

  1048. !
    1049. 

  1049.       zolri2=zol2*psih2/psix2**2-ri2
    1050. 

  1050. !
    1051. 

  1051.       return
    1052. 

  1052.       end
    1053. 

  1053. !
    1054. 

  1054. ! ... integrated similarity functions ...
    1055. 

  1055. !
    1056. 

  1056.       function psim_stable(zolf)
    1057. 

  1057.         psim_stable=-6.1*log(zolf+(1+zolf**2.5)**(1./2.5))
    1058. 

  1058.       return
    1059. 

  1059.       end
    1060. 

  1060. 

  1061.       function psih_stable(zolf)
    1062. 

  1062.         psih_stable=-5.3*log(zolf+(1+zolf**1.1)**(1./1.1))
    1063. 

  1063.       return
    1064. 

  1064.       end
    1065. 

  1065.       
    1066. 

  1066.       function psim_unstable(zolf)
    1067. 

  1067.         x=(1.-16.*zolf)**.25
    1068. 

  1068.         psimk=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))-2.*ATAN(X)+2.*ATAN(1.)
    1069. 

  1069. !
    1070. 

  1070.         ym=(1.-10.*zolf)**0.33
    1071. 

  1071.         psimc=(3./2.)*log((ym**2.+ym+1.)/3.)-sqrt(3.)*ATAN((2.*ym+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.)
    1072. 

  1072. !
    1073. 

  1073.         psim_unstable=(psimk+zolf**2*(psimc))/(1+zolf**2.)
    1074. 

  1074. 

  1075.       return
    1076. 

  1076.       end
    1077. 

  1077. 

  1078.       function psih_unstable(zolf)
    1079. 

  1079.         y=(1.-16.*zolf)**.5
    1080. 

  1080.         psihk=2.*log((1+y)/2.)
    1081. 

  1081. !
    1082. 

  1082.         yh=(1.-34.*zolf)**0.33
    1083. 

  1083.         psihc=(3./2.)*log((yh**2.+yh+1.)/3.)-sqrt(3.)*ATAN((2.*yh+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.)
    1084. 

  1084. !
    1085. 

  1085.         psih_unstable=(psihk+zolf**2*(psihc))/(1+zolf**2.)
    1086. 

  1086. 

  1087.       return
    1088. 

  1088.       end
    1089. 

  1089. 

  1090. 

  1091. 

  1092.