PageRenderTime 35ms CodeModel.GetById 0ms RepoModel.GetById 0ms app.codeStats 0ms

/wrfv2_fire/phys/module_sf_gfdl.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 1641 lines | 999 code | 185 blank | 457 comment | 68 complexity | 7b4d7cc1de1a9f359ac0c43e4bda38d3 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. !-- XLV latent heat of vaporization for water (J/kg)
  2. !
  3. MODULE module_sf_gfdl
  4. !real, dimension(-100:2000,-100:2000), save :: z00
  5. CONTAINS
  6. !-------------------------------------------------------------------
  7. SUBROUTINE SF_GFDL(U3D,V3D,T3D,QV3D,P3D, &
  8. CP,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2, CPM, &
  9. DT, SMOIS,num_soil_layers,ISLTYP,ZNT,UST,PSIM,PSIH, &
  10. XLAND,HFX,QFX,TAUX,TAUY,LH,GSW,GLW,TSK,FLHC,FLQC, & ! gopal's doing for Ocean coupling
  11. QGH,QSFC,U10,V10, &
  12. GZ1OZ0,WSPD,BR,ISFFLX, &
  13. EP1,EP2,KARMAN,NTSFLG,SFENTH, &
  14. ids,ide, jds,jde, kds,kde, &
  15. ims,ime, jms,jme, kms,kme, &
  16. its,ite, jts,jte, kts,kte )
  17. !-------------------------------------------------------------------
  18. USE MODULE_GFS_MACHINE, ONLY : kind_phys
  19. USE MODULE_GFS_FUNCPHYS , ONLY : gfuncphys,fpvs
  20. USE MODULE_GFS_PHYSCONS, grav => con_g
  21. !-------------------------------------------------------------------
  22. IMPLICIT NONE
  23. !-------------------------------------------------------------------
  24. !-- U3D 3D u-velocity interpolated to theta points (m/s)
  25. !-- V3D 3D v-velocity interpolated to theta points (m/s)
  26. !-- T3D temperature (K)
  27. !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
  28. !-- P3D 3D pressure (Pa)
  29. !-- DT time step (second)
  30. !-- CP heat capacity at constant pressure for dry air (J/kg/K)
  31. !-- ROVCP R/CP
  32. !-- R gas constant for dry air (J/kg/K)
  33. !-- XLV latent heat of vaporization for water (J/kg)
  34. !-- PSFC surface pressure (Pa)
  35. !-- ZNT roughness length (m)
  36. !-- MAVAIL surface moisture availability (between 0 and 1)
  37. !-- UST u* in similarity theory (m/s)
  38. !-- PSIM similarity stability function for momentum
  39. !-- PSIH similarity stability function for heat
  40. !-- XLAND land mask (1 for land, 2 for water)
  41. !-- HFX upward heat flux at the surface (W/m^2)
  42. !-- QFX upward moisture flux at the surface (kg/m^2/s)
  43. !-- TAUX RHO*U**2 (Kg/m/s^2) ! gopal's doing for Ocean coupling
  44. !-- TAUY RHO*U**2 (Kg/m/s^2) ! gopal's doing for Ocean coupling
  45. !-- LH net upward latent heat flux at surface (W/m^2)
  46. !-- GSW downward short wave flux at ground surface (W/m^2)
  47. !-- GLW downward long wave flux at ground surface (W/m^2)
  48. !-- TSK surface temperature (K)
  49. !-- FLHC exchange coefficient for heat (m/s)
  50. !-- FLQC exchange coefficient for moisture (m/s)
  51. !-- QGH lowest-level saturated mixing ratio
  52. !-- U10 diagnostic 10m u wind
  53. !-- V10 diagnostic 10m v wind
  54. !-- GZ1OZ0 log(z/z0) where z0 is roughness length
  55. !-- WSPD wind speed at lowest model level (m/s)
  56. !-- BR bulk Richardson number in surface layer
  57. !-- ISFFLX isfflx=1 for surface heat and moisture fluxes
  58. !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
  59. !-- KARMAN Von Karman constant
  60. !-- SFENTH enthalpy flux factor 0 zot via charnock ..>0 zot enhanced>15m/s
  61. !-- ids start index for i in domain
  62. !-- ide end index for i in domain
  63. !-- jds start index for j in domain
  64. !-- jde end index for j in domain
  65. !-- kds start index for k in domain
  66. !-- kde end index for k in domain
  67. !-- ims start index for i in memory
  68. !-- ime end index for i in memory
  69. !-- jms start index for j in memory
  70. !-- jme end index for j in memory
  71. !-- kms start index for k in memory
  72. !-- kme end index for k in memory
  73. !-- its start index for i in tile
  74. !-- ite end index for i in tile
  75. !-- jts start index for j in tile
  76. !-- jte end index for j in tile
  77. !-- kts start index for k in tile
  78. !-- kte end index for k in tile
  79. !-------------------------------------------------------------------
  80. INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
  81. ims,ime, jms,jme, kms,kme, &
  82. its,ite, jts,jte, kts,kte, &
  83. ISFFLX,NUM_SOIL_LAYERS,NTSFLG
  84. REAL, INTENT(IN) :: &
  85. CP, &
  86. EP1, &
  87. EP2, &
  88. KARMAN, &
  89. R, &
  90. ROVCP, &
  91. DT, &
  92. SFENTH, &
  93. XLV
  94. REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: &
  95. P3D, &
  96. QV3D, &
  97. T3D, &
  98. U3D, &
  99. V3D
  100. INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: ISLTYP
  101. REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), INTENT(INOUT):: SMOIS
  102. REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: &
  103. PSFC, &
  104. GLW, &
  105. GSW, &
  106. XLAND
  107. REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: &
  108. TSK, &
  109. BR, &
  110. CHS, &
  111. CHS2, &
  112. CPM, &
  113. CQS2, &
  114. FLHC, &
  115. FLQC, &
  116. GZ1OZ0, &
  117. HFX, &
  118. LH, &
  119. PSIM, &
  120. PSIH, &
  121. QFX, &
  122. QGH, &
  123. QSFC, &
  124. UST, &
  125. ZNT, &
  126. WSPD, &
  127. TAUX, & ! gopal's doing for Ocean coupling
  128. TAUY
  129. REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: &
  130. U10, &
  131. V10
  132. !--------------------------- LOCAL VARS ------------------------------
  133. REAL :: ESAT, &
  134. cpcgs, &
  135. smc, &
  136. smcdry, &
  137. smcmax
  138. REAL (kind=kind_phys) :: &
  139. RHOX
  140. REAL, DIMENSION(1:30) :: MAXSMC, &
  141. DRYSMC
  142. REAL (kind=kind_phys), DIMENSION(its:ite) :: &
  143. CH, &
  144. CM, &
  145. DDVEL, &
  146. DRAIN, &
  147. EP, &
  148. EVAP, &
  149. FH, &
  150. FH2, &
  151. FM, &
  152. HFLX, &
  153. PH, &
  154. PM, &
  155. PRSL1, &
  156. PRSLKI, &
  157. PS, &
  158. Q1, &
  159. Q2M, &
  160. QSS, &
  161. QSURF, &
  162. RB, &
  163. RCL, &
  164. RHO1, &
  165. SLIMSK, &
  166. STRESS, &
  167. T1, &
  168. T2M, &
  169. THGB, &
  170. THX, &
  171. TSKIN, &
  172. SHELEG, &
  173. U1, &
  174. U10M, &
  175. USTAR, &
  176. V1, &
  177. V10M, &
  178. WIND, &
  179. Z0RL, &
  180. Z1
  181. REAL, DIMENSION(kms:kme, ims:ime) :: &
  182. rpc, &
  183. tpc, &
  184. upc, &
  185. vpc
  186. REAL, DIMENSION(ims:ime) :: &
  187. pspc, &
  188. pkmax, &
  189. tstrc, &
  190. zoc, &
  191. wetc, &
  192. slwdc, &
  193. rib, &
  194. zkmax, &
  195. tkmax, &
  196. fxmx, &
  197. fxmy, &
  198. cdm, &
  199. fxh, &
  200. fxe, &
  201. xxfh, &
  202. xxfh2, &
  203. wind10, &
  204. tjloc
  205. INTEGER :: &
  206. I, &
  207. II, &
  208. IGPVS, &
  209. IM, &
  210. J, &
  211. K, &
  212. KM
  213. DATA MAXSMC/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
  214. 0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
  215. 0.439, 1.000, 0.200, 0.421, 0.000, 0.000, &
  216. 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
  217. 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
  218. DATA DRYSMC/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, &
  219. 0.067, 0.120, 0.103, 0.100, 0.126, 0.138, &
  220. 0.066, 0.000, 0.006, 0.028, 0.000, 0.000, &
  221. 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, &
  222. 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
  223. DATA IGPVS/0/
  224. save igpvs
  225. if(igpvs.eq.0) then
  226. ! call readzo(glat,glon,6,ims,ime,jms,jme,its,ite,jts,jte,z00)
  227. endif
  228. igpvs=1
  229. IM=ITE-ITS+1
  230. KM=KTE-KTS+1
  231. WRITE(0,*)'WITHIN THE GFDL SCHEME, NTSFLG=1 FOR GFDL SLAB 2010 UPGRADS',NTSFLG
  232. DO J=jts,jte
  233. DO i=its,ite
  234. DDVEL(I)=0.
  235. RCL(i)=1.
  236. PRSL1(i)=P3D(i,kts,j)*.001
  237. wetc(i)=1.0
  238. if(xland(i,j).lt.1.99) then
  239. smc=smois(i,1,j)
  240. smcdry=drysmc(isltyp(i,j))
  241. smcmax=maxsmc(isltyp(i,j))
  242. wetc(i)=(smc-smcdry)/(smcmax-smcdry)
  243. wetc(i)=amin1(1.,amax1(wetc(i),0.))
  244. endif
  245. ! convert from Pa to cgs...
  246. pspc(i)=PSFC(i,j)*10.
  247. pkmax(i)=P3D(i,kts,j)*10.
  248. PS(i)=PSFC(i,j)*.001
  249. Q1(I) = QV3D(i,kts,j)
  250. rpc(kts,i)=QV3D(i,kts,j)
  251. ! QSURF(I)=QSFC(I,J)
  252. QSURF(I)=0.
  253. SHELEG(I)=0.
  254. SLIMSK(i)=ABS(XLAND(i,j)-2.)
  255. TSKIN(i)=TSK(i,j)
  256. tstrc(i)=TSK(i,j)
  257. T1(I) = T3D(i,kts,j)
  258. tpc(kts,i)=T3D(i,kts,j)
  259. U1(I) = U3D(i,kts,j)
  260. upc(kts,i)=U3D(i,kts,j) * 100.
  261. USTAR(I) = UST(i,j)
  262. V1(I) = V3D(i,kts,j)
  263. vpc(kts,i)=v3D(i,kts,j) * 100.
  264. Z0RL(I) = ZNT(i,j)*100.
  265. zoc(i)=ZNT(i,j)*100.
  266. if(XLAND(i,j).gt.1.99) zoc(i)=- zoc(i)
  267. ! Z0RL(I) = z00(i,j)*100.
  268. ! slwdc... GFDL downward net flux in units of cal/(cm**2/min)
  269. ! also divide by 10**4 to convert from /m**2 to /cm**2
  270. slwdc(i)=gsw(i,j)+glw(i,j)
  271. slwdc(i)=0.239*60.*slwdc(i)*1.e-4
  272. tjloc(i)=float(j)
  273. ENDDO
  274. DO i=its,ite
  275. PRSLKI(i)=(PS(I)/PRSL1(I))**ROVCP
  276. THGB(I)=TSKIN(i)*(100./PS(I))**ROVCP
  277. THX(I)=T1(i)*(100./PRSL1(I))**ROVCP
  278. RHO1(I)=PRSL1(I)*1000./(R*T1(I)*(1.+EP1*Q1(I)))
  279. Q1(I)=Q1(I)/(1.+Q1(I))
  280. ENDDO
  281. ! if(j==2)then
  282. ! write(0,*)'--------------------------------------------'
  283. ! write(0,*) 'u, v, t, r, pkmax, pspc,wetc, tjloc,zoc,tstr'
  284. ! write(0,*)'--------------------------------------------'
  285. ! endif
  286. ! do i = its,ite
  287. ! WRITE(0,1010)i,j,upc(kts,i),vpc(kts,i),tpc(kts,i),rpc(kts,i), &
  288. ! pkmax(i),pspc(i),wetc(i),tjloc(i),zoc(i),tstrc(i)
  289. ! enddo
  290. CALL MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,tstrc, &
  291. pspc,pkmax,wetc,slwdc,tjloc, &
  292. upc,vpc,tpc,rpc,dt,J,wind10,xxfh2,ntsflg,SFENTH, &
  293. ids,ide, jds,jde, kds,kde, &
  294. ims,ime, jms,jme, kms,kme, &
  295. its,ite, jts,jte, kts,kte )
  296. ! if(j==2)then
  297. ! write(0,*)'--------------------------------------------'
  298. ! write(0,*) 'fxh, fxe, fxmx, fxmy, cdm, xxfh zoc,tstrc'
  299. ! write(0,*)'--------------------------------------------'
  300. ! endif
  301. ! do i = its,ite
  302. ! WRITE(0,1010)i,j,fxh(i),fxe(i),fxmx(i),fxmy(i),cdm(i),rib(i),xxfh(i),zoc(i),tstrc(i)
  303. ! enddo
  304. 1010 format(2I4,9F11.6)
  305. !GFDL CALL PROGTM(IM,KM,PS,U1,V1,T1,Q1, &
  306. !GFDL SHELEG,TSKIN,QSURF, &
  307. !WRF SMC,STC,DM,SOILTYP,SIGMAF,VEGTYPE,CANOPY,DLWFLX, &
  308. !WRF SLRAD,SNOWMT,DELT, &
  309. !GFDL Z0RL, &
  310. !WRF TG3,GFLUX,F10M, &
  311. !GFDL U10M,V10M,T2M,Q2M, &
  312. !WRF ZSOIL, &
  313. !GFDL CM,CH,RB, &
  314. !WRF RHSCNPY,RHSMC,AIM,BIM,CIM, &
  315. !GFDL RCL,PRSL1,PRSLKI,SLIMSK, &
  316. !GFDL DRAIN,EVAP,HFLX,STRESS,EP, &
  317. !GFDL FM,FH,USTAR,WIND,DDVEL, &
  318. !GFDL PM,PH,FH2,QSS,Z1 )
  319. DO i=its,ite
  320. ! update skin temp only when using GFDL slab...
  321. IF(NTSFLG==1) then
  322. tsk(i,j) = tstrc(i) ! gopal's doing
  323. !bugfix 4
  324. ! bob's doing patch tsk with neigboring values... are grid boundaries
  325. if(j.eq.jde) then
  326. tsk(i,j) = tsk(i,j-1)
  327. endif
  328. if(j.eq.jds) then
  329. tsk(i,j) = tsk(i,j+1)
  330. endif
  331. if(i.eq.ide) tsk(i,j) = tsk(i-1,j)
  332. if(i.eq.ids) tsk(i,j) = tsk(i+1,j)
  333. endif
  334. znt(i,j)= 0.01*abs(zoc(i))
  335. wspd(i,j) = SQRT(upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i))
  336. wspd(i,j) = amax1(wspd(i,j) ,100.)/100.
  337. u10m(i) = u1(i)*(wind10(i)/wspd(i,j))/100.
  338. v10m(i) = v1(i)*(wind10(i)/wspd(i,j))/100.
  339. ! br =0.0001*zfull(i,kmax)*dthv/
  340. ! & (gmks*theta(i,kmax)*wspd *wspd )
  341. ! zkmax = rgas*tpc(kmax,i)*qqlog(kmax)*og
  342. zkmax(i) = -R*tpc(kts,i)*alog(pkmax(i)/pspc(i))/grav
  343. !------------------------------------------------------------------------
  344. gz1oz0(i,j)=alog(zkmax(i)/znt(i,j))
  345. ustar (i)= 0.01*sqrt(cdm(i)* &
  346. (upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i)))
  347. ! convert from g/(cm*cm*sec) to kg/(m*m*sec)
  348. qfx (i,j)=-10.*fxe(i) ! BOB: qfx (i,1)=-10.*fxe(i)
  349. ! cpcgs = 1.00464e7
  350. ! convert from ergs/gram/K to J/kg/K cpmks=1004
  351. ! hfx (i,1)=-0.001*cpcgs*fxh(i)
  352. hfx (i,j)= -10.*CP*fxh(i) ! Bob: hfx (i,1)= -10.*CP*fxh(i)
  353. taux (i,j)= fxmx(i)/10. ! gopal's doing for Ocean coupling
  354. tauy (i,j)= fxmy(i)/10. ! gopal's doing for Ocean coupling
  355. fm(i) = karman/sqrt(cdm(i))
  356. fh(i) = karman*xxfh(i)
  357. PSIM(i,j)=GZ1OZ0(i,j)-FM(i)
  358. PSIH(i,j)=GZ1OZ0(i,j)-FH(i)
  359. fh2(i) = karman*xxfh2(i)
  360. ch(i) = karman*karman/(fm(i) * fh(i))
  361. cm(i) = cdm(i)
  362. U10(i,j)=U10M(i)
  363. V10(i,j)=V10M(i)
  364. BR(i,j)=rib(i)
  365. CHS(I,J)=CH(I)*wspd (i,j)
  366. CHS2(I,J)=USTAR(I)*KARMAN/FH2(I)
  367. CPM(I,J)=CP*(1.+0.8*QV3D(i,kts,j))
  368. esat = fpvs(t1(i))
  369. QGH(I,J)=ep2*esat/(1000.*ps(i)-esat)
  370. esat = fpvs(tskin(i))
  371. qss(i) = ep2*esat/(1000.*ps(i)-esat)
  372. QSFC(I,J)=qss(i)
  373. ! PSIH(i,j)=PH(i)
  374. ! PSIM(i,j)=PM(i)
  375. UST(i,j)=ustar(i)
  376. ! wspd (i,j) = SQRT(upc(kts,i)*upc(kts,i) + vpc(kts,i)*vpc(kts,i))
  377. ! wspd (i,j) = amax1(wspd (i,j) ,100.)/100.
  378. ! WSPD(i,j)=WIND(i)
  379. ! ZNT(i,j)=Z0RL(i)*.01
  380. ENDDO
  381. ! write(0,*)'fm,fh,cm,ch(125)', fm(125),fh(125),cm(125),ch(125)
  382. DO i=its,ite
  383. FLHC(i,j)=CPM(I,J)*RHO1(I)*CHS(I,J)
  384. FLQC(i,j)=RHO1(I)*CHS(I,J)
  385. ! GZ1OZ0(i,j)=LOG(Z1(I)/(Z0RL(I)*.01))
  386. CQS2(i,j)=CHS2(I,J)
  387. ENDDO
  388. IF (ISFFLX.EQ.0) THEN
  389. DO i=its,ite
  390. HFX(i,j)=0.
  391. LH(i,j)=0.
  392. QFX(i,j)=0.
  393. ENDDO
  394. ELSE
  395. DO i=its,ite
  396. IF(XLAND(I,J)-1.5.GT.0.)THEN
  397. ! HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
  398. ! cpcgs = 1.00464e7
  399. ! convert from ergs/gram/K to J/kg/K cpmks=1004
  400. ! hfx (i,j)=-0.001*cpcgs*fxh(i)
  401. hfx (i,j)= -10.*CP*fxh(i) ! Bob: hfx (i,1)= -10.*CP*fxh(i)
  402. ELSEIF(XLAND(I,J)-1.5.LT.0.)THEN
  403. ! HFX(I,J)=FLHC(I,J)*(THGB(I)-THX(I))
  404. ! cpcgs = 1.00464e7
  405. ! convert from ergs/gram/K to J/kg/K cpmks=1004
  406. ! hfx (i,j)=-0.001*cpcgs*fxh(i)
  407. hfx (i,j)= -10.*CP*fxh(i) ! Bob: hfx (i,j)= -10.*CP*fxh(i)
  408. HFX(I,J)=AMAX1(HFX(I,J),-250.)
  409. ENDIF
  410. ! QFX(I,J)=FLQC(I,J)*(QSFC(I,J)-Q1(I))
  411. ! convert from g/(cm*cm*sec) to kg/(m*m*sec)
  412. qfx(i,j)=-10.*fxe(i)
  413. QFX(I,J)=AMAX1(QFX(I,J),0.)
  414. LH(I,J)=XLV*QFX(I,J)
  415. ENDDO
  416. ENDIF
  417. ! if(j.eq.2) write(0,*) 'u3d ,ustar,cdm at end of gfdlsfcmod'
  418. ! write(0,*) j,(u3d(ii,1,j),ii=70,90)
  419. ! write(0,*) j,(ustar(ii),ii=70,90)
  420. ! write(0,*) j,(cdm(ii),ii=70,90)
  421. if(j.eq.jds.or.j.eq.jde) then
  422. write(0,*) "TSFC in gfdl sf mod,dt, its,ite,jts,jts", dt,its,ite,jts,jte,ids,ide,jds,jde
  423. write(0,*) "TSFC", (TSK(i,j),i=its,ite)
  424. endif
  425. ENDDO ! FOR THE J LOOP I PRESUME
  426. ! if(100.ge.its.and.100.le.ite.and.100.ge.jts.and.100.le.jte) then
  427. ! write(0,*) 'output vars of sf_gfdl at i,j=100'
  428. ! write(0,*) 'TSK',TSK(100,100)
  429. ! write(0,*) 'PSFC',PSFC(100,100)
  430. ! write(0,*) 'GLW',GLW(100,100)
  431. ! write(0,*) 'GSW',GSW(100,100)
  432. ! write(0,*) 'XLAND',XLAND(100,100)
  433. ! write(0,*) 'BR',BR(100,100)
  434. ! write(0,*) 'CHS',CHS(100,100)
  435. ! write(0,*) 'CHS2',CHS2(100,100)
  436. ! write(0,*) 'CPM',CPM(100,100)
  437. ! write(0,*) 'FLHC',FLHC(100,100)
  438. ! write(0,*) 'FLQC',FLQC(100,100)
  439. ! write(0,*) 'GZ1OZ0',GZ1OZ0(100,100)
  440. ! write(0,*) 'HFX',HFX(100,100)
  441. ! write(0,*) 'LH',LH(100,100)
  442. ! write(0,*) 'PSIM',PSIM(100,100)
  443. ! write(0,*) 'PSIH',PSIH(100,100)
  444. ! write(0,*) 'QFX',QFX(100,100)
  445. ! write(0,*) 'QGH',QGH(100,100)
  446. ! write(0,*) 'QSFC',QSFC(100,100)
  447. ! write(0,*) 'UST',UST(100,100)
  448. ! write(0,*) 'ZNT',ZNT(100,100)
  449. ! write(0,*) 'wet',wet(100)
  450. ! write(0,*) 'smois',smois(100,1,100)
  451. ! write(0,*) 'WSPD',WSPD(100,100)
  452. ! write(0,*) 'U10',U10(100,100)
  453. ! write(0,*) 'V10',V10(100,100)
  454. ! endif
  455. END SUBROUTINE SF_GFDL
  456. !-------------------------------------------------------------------
  457. SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,tstrc, &
  458. pspc,pkmax,wetc,slwdc,tjloc, &
  459. upc,vpc,tpc,rpc,dt,jfix,wind10,xxfh2,ntsflg,sfenth, &
  460. ids,ide, jds,jde, kds,kde, &
  461. ims,ime, jms,jme, kms,kme, &
  462. its,ite, jts,jte, kts,kte )
  463. !------------------------------------------------------------------------
  464. !
  465. ! MFLUX2 computes surface fluxes of momentum, heat,and moisture
  466. ! using monin-obukhov. the roughness length "z0" is prescribed
  467. ! over land and over ocean "z0" is computed using charnocks formula.
  468. ! the universal functions (from similarity theory approach) are
  469. ! those of hicks. This is Bob's doing.
  470. !
  471. !------------------------------------------------------------------------
  472. IMPLICIT NONE
  473. ! use allocate_mod
  474. ! use module_TLDATA , ONLY : tab,table,cp,g,rgas,og
  475. ! include 'RESOLUTION.h'
  476. ! include 'PARAMETERS.h'
  477. ! include 'STDUNITS.h' stdout
  478. ! include 'FLAGS.h'
  479. ! include 'BKINFO.h' nstep
  480. ! include 'CONSLEV.h'
  481. ! include 'CONMLEV.h'
  482. ! include 'ESTAB.h'
  483. ! include 'FILEC.h'
  484. ! include 'FILEPC.h'
  485. ! include 'GDINFO.h' ngd
  486. ! include 'LIMIT.h'
  487. ! include 'QLOGS.h'
  488. ! include 'TIME.h' dt(nnst)
  489. ! include 'WINDD.h'
  490. ! include 'ZLDATA.h' old MOBFLX?
  491. !-----------------------------------------------------------------------
  492. ! user interface variables
  493. !-----------------------------------------------------------------------
  494. integer,intent(in) :: ids,ide, jds,jde, kds,kde
  495. integer,intent(in) :: ims,ime, jms,jme, kms,kme
  496. integer,intent(in) :: its,ite, jts,jte, kts,kte
  497. integer,intent(in) :: jfix,ntsflg
  498. real, intent (out), dimension (ims :ime ) :: fxh
  499. real, intent (out), dimension (ims :ime ) :: fxe
  500. real, intent (out), dimension (ims :ime ) :: fxmx
  501. real, intent (out), dimension (ims :ime ) :: fxmy
  502. real, intent (out), dimension (ims :ime ) :: cdm
  503. ! real, intent (out), dimension (ims :ime ) :: cdm2
  504. real, intent (out), dimension (ims :ime ) :: rib
  505. real, intent (out), dimension (ims :ime ) :: xxfh
  506. real, intent (out), dimension (ims :ime ) :: xxfh2
  507. real, intent (out), dimension (ims :ime ) :: wind10
  508. real, intent ( inout), dimension (ims :ime ) :: zoc
  509. real, intent ( inout), dimension (ims :ime ) :: tstrc
  510. real, intent ( in) :: dt
  511. real, intent ( in) :: sfenth
  512. real, intent ( in), dimension (ims :ime ) :: pspc
  513. real, intent ( in), dimension (ims :ime ) :: pkmax
  514. real, intent ( in), dimension (ims :ime ) :: wetc
  515. real, intent ( in), dimension (ims :ime ) :: slwdc
  516. real, intent ( in), dimension (ims :ime ) :: tjloc
  517. real, intent ( in), dimension (kms:kme, ims :ime ) :: upc
  518. real, intent ( in), dimension (kms:kme, ims :ime ) :: vpc
  519. real, intent ( in), dimension (kms:kme, ims :ime ) :: tpc
  520. real, intent ( in), dimension (kms:kme, ims :ime ) :: rpc
  521. !-----------------------------------------------------------------------
  522. ! internal variables
  523. !-----------------------------------------------------------------------
  524. integer, parameter :: icntx = 30
  525. integer, dimension(1 :ime) :: ifz
  526. integer, dimension(1 :ime) :: indx
  527. integer, dimension(1 :ime) :: istb
  528. integer, dimension(1 :ime) :: it
  529. integer, dimension(1 :ime) :: iutb
  530. real, dimension(1 :ime) :: aap
  531. real, dimension(1 :ime) :: bq1
  532. real, dimension(1 :ime) :: bq1p
  533. real, dimension(1 :ime) :: delsrad
  534. real, dimension(1 :ime) :: ecof
  535. real, dimension(1 :ime) :: ecofp
  536. real, dimension(1 :ime) :: estso
  537. real, dimension(1 :ime) :: estsop
  538. real, dimension(1 :ime) :: fmz1
  539. real, dimension(1 :ime) :: fmz10
  540. real, dimension(1 :ime) :: fmz2
  541. real, dimension(1 :ime) :: fmzo1
  542. real, dimension(1 :ime) :: foft
  543. real, dimension(1 :ime) :: foftm
  544. real, dimension(1 :ime) :: frac
  545. real, dimension(1 :ime) :: land
  546. real, dimension(1 :ime) :: pssp
  547. real, dimension(1 :ime) :: qf
  548. real, dimension(1 :ime) :: rdiff
  549. real, dimension(1 :ime) :: rho
  550. real, dimension(1 :ime) :: rkmaxp
  551. real, dimension(1 :ime) :: rstso
  552. real, dimension(1 :ime) :: rstsop
  553. real, dimension(1 :ime) :: sf10
  554. real, dimension(1 :ime) :: sf2
  555. real, dimension(1 :ime) :: sfm
  556. real, dimension(1 :ime) :: sfzo
  557. real, dimension(1 :ime) :: sgzm
  558. real, dimension(1 :ime) :: slwa
  559. real, dimension(1 :ime) :: szeta
  560. real, dimension(1 :ime) :: szetam
  561. real, dimension(1 :ime) :: t1
  562. real, dimension(1 :ime) :: t2
  563. real, dimension(1 :ime) :: tab1
  564. real, dimension(1 :ime) :: tab2
  565. real, dimension(1 :ime) :: tempa1
  566. real, dimension(1 :ime) :: tempa2
  567. real, dimension(1 :ime) :: theta
  568. real, dimension(1 :ime) :: thetap
  569. real, dimension(1 :ime) :: tsg
  570. real, dimension(1 :ime) :: tsm
  571. real, dimension(1 :ime) :: tsp
  572. real, dimension(1 :ime) :: tss
  573. real, dimension(1 :ime) :: ucom
  574. real, dimension(1 :ime) :: uf10
  575. real, dimension(1 :ime) :: uf2
  576. real, dimension(1 :ime) :: ufh
  577. real, dimension(1 :ime) :: ufm
  578. real, dimension(1 :ime) :: ufzo
  579. real, dimension(1 :ime) :: ugzm
  580. real, dimension(1 :ime) :: uzeta
  581. real, dimension(1 :ime) :: uzetam
  582. real, dimension(1 :ime) :: vcom
  583. real, dimension(1 :ime) :: vrtkx
  584. real, dimension(1 :ime) :: vrts
  585. real, dimension(1 :ime) :: wind
  586. real, dimension(1 :ime) :: windp
  587. ! real, dimension(1 :ime) :: xxfh
  588. real, dimension(1 :ime) :: xxfm
  589. real, dimension(1 :ime) :: xxsh
  590. real, dimension(1 :ime) :: z10
  591. real, dimension(1 :ime) :: z2
  592. real, dimension(1 :ime) :: zeta
  593. real, dimension(1 :ime) :: zkmax
  594. real, dimension(1 :ime) :: pss
  595. real, dimension(1 :ime) :: tstar
  596. real, dimension(1 :ime) :: ukmax
  597. real, dimension(1 :ime) :: vkmax
  598. real, dimension(1 :ime) :: tkmax
  599. real, dimension(1 :ime) :: rkmax
  600. real, dimension(1 :ime) :: zot
  601. real, dimension(1 :ime) :: fhzo1
  602. real, dimension(1 :ime) :: sfh
  603. real :: ux13, yo, y,xo,x,ux21,ugzzo,ux11,ux12,uzetao,xnum,alll
  604. real :: ux1,ugz,x10,uzo,uq,ux2,ux3,xtan,xden,y10,uzet1o,ugz10
  605. real :: szet2, zal2,ugz2
  606. real :: rovcp,boycon,cmo2,psps1,zog,enrca,rca,cmo1,amask,en,ca,a,c
  607. real :: sgz,zal10,szet10,fmz,szo,sq,fmzo,rzeta1,zal1g,szetao,rzeta2,zal2g
  608. real :: hcap,xks,pith,teps,diffot,delten,alevp,psps2,alfus,nstep
  609. real :: shfx,sigt4,reflect
  610. real :: cor1,cor2,szetho,zal2gh,cons_p000001,cons_7,vis,ustar,restar,rat
  611. real :: wndm,ckg
  612. real :: yz,y1,y2,y3,y4,windmks,znott,znotm
  613. integer:: i,j,ii,iq,nnest,icnt,ngd,ip
  614. !-----------------------------------------------------------------------
  615. ! internal variables
  616. !-----------------------------------------------------------------------
  617. real, dimension (223) :: tab
  618. real, dimension (223) :: table
  619. real, dimension (101) :: tab11
  620. real, dimension (41) :: table4
  621. real, dimension (42) :: tab3
  622. real, dimension (54) :: table2
  623. real, dimension (54) :: table3
  624. real, dimension (74) :: table1
  625. real, dimension (80) :: tab22
  626. equivalence (tab(1),tab11(1))
  627. equivalence (tab(102),tab22(1))
  628. equivalence (tab(182),tab3(1))
  629. equivalence (table(1),table1(1))
  630. equivalence (table(75),table2(1))
  631. equivalence (table(129),table3(1))
  632. equivalence (table(183),table4(1))
  633. data amask/ -98.0/
  634. !-----------------------------------------------------------------------
  635. ! tables used to obtain the vapor pressures or saturated vapor
  636. ! pressure
  637. !-----------------------------------------------------------------------
  638. data tab11/21*0.01403,0.01719,0.02101,0.02561,0.03117,0.03784, &
  639. &.04584,.05542,.06685,.08049,.09672,.1160,.1388,.1658,.1977,.2353, &
  640. &.2796,.3316,.3925,.4638,.5472,.6444,.7577,.8894,1.042,1.220,1.425, &
  641. &1.662,1.936,2.252,2.615,3.032,3.511,4.060,4.688,5.406,6.225,7.159, &
  642. &8.223,9.432,10.80,12.36,14.13,16.12,18.38,20.92,23.80,27.03,30.67, &
  643. &34.76,39.35,44.49,50.26,56.71,63.93,71.98,80.97,90.98,102.1,114.5, &
  644. &128.3,143.6,160.6,179.4,200.2,223.3,248.8,276.9,307.9,342.1,379.8, &
  645. &421.3,466.9,517.0,572.0,632.3,698.5,770.9,850.2,937.0,1032./
  646. data tab22/1146.6,1272.0,1408.1,1556.7,1716.9,1890.3,2077.6,2279.6 &
  647. &,2496.7,2729.8,2980.0,3247.8,3534.1,3839.8,4164.8,4510.5,4876.9, &
  648. &5265.1,5675.2,6107.8,6566.2,7054.7,7575.3,8129.4,8719.2,9346.5, &
  649. &10013.,10722.,11474.,12272.,13119.,14017.,14969.,15977.,17044., &
  650. &18173.,19367.,20630.,21964.,23373.,24861.,26430.,28086.,29831., &
  651. &31671.,33608.,35649.,37796.,40055.,42430.,44927.,47551.,50307., &
  652. &53200.,56236.,59422.,62762.,66264.,69934.,73777.,77802.,82015., &
  653. &86423.,91034.,95855.,100890.,106160.,111660.,117400.,123400., &
  654. &129650.,136170.,142980.,150070.,157460.,165160.,173180.,181530., &
  655. &190220.,199260./
  656. data tab3/208670.,218450.,228610.,239180.,250160.,261560.,273400., &
  657. &285700.,298450.,311690.,325420.,339650.,354410.,369710.,385560., &
  658. &401980.,418980.,436590.,454810.,473670.,493170.,513350.,534220., &
  659. &555800.,578090.,601130.,624940.,649530.,674920.,701130.,728190., &
  660. &756110.,784920.,814630.,845280.,876880.,909450.,943020.,977610., &
  661. &1013250.,1049940.,1087740./
  662. data table1/20*0.0,.3160e-02,.3820e-02,.4600e-02,.5560e-02,.6670e-02, &
  663. & .8000e-02,.9580e-02,.1143e-01,.1364e-01,.1623e-01,.1928e-01, &
  664. &.2280e-01,.2700e-01,.3190e-01,.3760e-01,.4430e-01,.5200e-01, &
  665. &.6090e-01,.7130e-01,.8340e-01,.9720e-01,.1133e+00,.1317e-00, &
  666. &.1526e-00,.1780e-00,.2050e-00,.2370e-00,.2740e-00,.3160e-00, &
  667. &.3630e-00,.4170e-00,.4790e-00,.5490e-00,.6280e-00,.7180e-00, &
  668. &.8190e-00,.9340e-00,.1064e+01,.1209e+01,.1368e+01,.1560e+01, &
  669. &.1770e+01,.1990e+01,.2260e+01,.2540e+01,.2880e+01,.3230e+01, &
  670. &.3640e+01,.4090e+01,.4590e+01,.5140e+01,.5770e+01,.6450e+01, &
  671. &.7220e+01/
  672. data table2/.8050e+01,.8990e+01,.1001e+02,.1112e+02,.1240e+02, &
  673. &.1380e+02,.1530e+02,.1700e+02,.1880e+02,.2080e+02,.2310e+02, &
  674. &.2550e+02,.2810e+02,.3100e+02,.3420e+02,.3770e+02,.4150e+02, &
  675. &.4560e+02,.5010e+02,.5500e+02,.6030e+02,.6620e+02,.7240e+02, &
  676. &.7930e+02,.8680e+02,.9500e+02,.1146e+03,.1254e+03,.1361e+03, &
  677. &.1486e+03,.1602e+03,.1734e+03,.1873e+03,.2020e+03,.2171e+03, &
  678. &.2331e+03,.2502e+03,.2678e+03,.2863e+03,.3057e+03,.3250e+03, &
  679. &.3457e+03,.3664e+03,.3882e+03,.4101e+03,.4326e+03,.4584e+03, &
  680. &.4885e+03,.5206e+03,.5541e+03,.5898e+03,.6273e+03,.6665e+03, &
  681. &.7090e+03/
  682. data table3/.7520e+03,.7980e+03,.8470e+03,.8980e+03,.9520e+03, &
  683. &.1008e+04,.1067e+04,.1129e+04,.1194e+04,.1263e+04,.1334e+04, &
  684. &.1409e+04,.1488e+04,.1569e+04,.1656e+04,.1745e+04,.1840e+04, &
  685. &.1937e+04,.2041e+04,.2147e+04,.2259e+04,.2375e+04,.2497e+04, &
  686. &.2624e+04,.2756e+04,.2893e+04,.3036e+04,.3186e+04,.3340e+04, &
  687. &.3502e+04,.3670e+04,.3843e+04,.4025e+04,.4213e+04,.4408e+04, &
  688. &.4611e+04,.4821e+04,.5035e+04,.5270e+04,.5500e+04,.5740e+04, &
  689. &.6000e+04,.6250e+04,.6520e+04,.6810e+04,.7090e+04,.7390e+04, &
  690. &.7700e+04,.8020e+04,.8350e+04,.8690e+04,.9040e+04,.9410e+04, &
  691. &.9780e+04/
  692. data table4/.1016e+05,.1057e+05,.1098e+05,.1140e+05,.1184e+05, &
  693. &.1230e+05,.1275e+05,.1324e+05,.1373e+05,.1423e+05,.1476e+05, &
  694. &.1530e+05,.1585e+05,.1642e+05,.1700e+05,.1761e+05,.1822e+05, &
  695. &.1886e+05,.1950e+05,.2018e+05,.2087e+05,.2158e+05,.2229e+05, &
  696. &.2304e+05,.2381e+05,.2459e+05,.2539e+05,.2621e+05,.2706e+05, &
  697. &.2792e+05,.2881e+05,.2971e+05,.3065e+05,.3160e+05,.3257e+05, &
  698. &.3357e+05,.3459e+05,.3564e+05,.3669e+05,.3780e+05,.0000e+00/
  699. !
  700. ! spcify constants needed by MFLUX2
  701. !
  702. real,parameter :: cp = 1.00464e7
  703. real,parameter :: g = 980.6
  704. real,parameter :: rgas = 2.87e6
  705. real,parameter :: og = 1./g
  706. !
  707. ! character*10 routine
  708. ! routine = 'mflux2'
  709. !
  710. !------------------------------------------------------------------------
  711. ! set water availability constant "ecof" and land mask "land".
  712. ! limit minimum wind speed to 100 cm/s
  713. !------------------------------------------------------------------------
  714. j = IFIX(tjloc(its))
  715. ! constants for 10 m winds (correction for knots
  716. !
  717. cor1 = .120
  718. cor2 = 720.
  719. yz= -0.0001344
  720. y1= 3.015e-05
  721. y2= 1.517e-06
  722. y3= -3.567e-08
  723. y4= 2.046e-10
  724. do i = its,ite
  725. z10(i) = 1000.
  726. z2 (i) = 200.
  727. pss(i) = pspc(i)
  728. tstar(i) = tstrc(i)
  729. ukmax(i) = upc(1,i)
  730. vkmax(i) = vpc(1,i)
  731. tkmax(i) = tpc(1,i)
  732. rkmax(i) = rpc(1,i)
  733. enddo
  734. ! write(0,*)'z10,pss,tstar,u...rkmax(ite)', &
  735. ! z10(ite), pss(ite),tstar(ite),ukmax(ite), &
  736. ! vkmax(ite),tkmax(ite),rkmax(ite)
  737. do i = its,ite
  738. windp(i) = SQRT(ukmax(i)*ukmax(i) + vkmax(i)*vkmax(i))
  739. wind (i) = amax1(windp(i),100.)
  740. if (zoc(i) .LT. amask) zoc(i) = -0.0185*0.001*wind(i)*wind(i)*og
  741. if (zoc(i) .GT. 0.0) then
  742. ecof(i) = wetc(i)
  743. land(i) = 1.0
  744. zot (i) = zoc(i)
  745. else
  746. ecof(i) = wetc(i)
  747. land(i) = 0.0
  748. zot (i) = zoc(i)
  749. ! now use 2 regime fit for znot thermal
  750. windmks=wind(i)*.01
  751. znott=0.2375*exp(-0.5250*windmks) + 0.0025*exp(-0.0211*windmks)
  752. znott=0.01*znott
  753. ! go back to moon et al for below 7m/s
  754. if(windmks.le. 7.) &
  755. znott = (0.0185/9.8*(7.59e-8*wind(i)**2+ &
  756. 2.46e-4*wind(i))**2)
  757. ! back to cgs
  758. zot (i) = 100.*znott
  759. ! end of kwon correction....
  760. ! in hwrf, thermal znot(zot) is passed as argument zoc
  761. ! in hwrf, momentum znot is recalculated internally
  762. zoc(i)=-(0.0185/9.8*(7.59e-8*wind(i)**2+ &
  763. 2.46e-4*wind(i))**2)*100.
  764. if(wind(i).ge.1250.0) &
  765. zoc(i)=-(.000739793 * wind(i) -0.58)/10
  766. if(wind(i).ge.3000.) then
  767. windmks=wind(i)*.01
  768. ! kwon znotm
  769. znotm = yz +windmks*y1 +windmks**2*y2 +windmks**3*y3 +windmks**4*y4 !powell 2003
  770. ! back to cgs
  771. zoc(i) = 100.*znotm
  772. endif
  773. endif
  774. !------------------------------------------------------------------------
  775. ! where necessary modify zo values over ocean.
  776. !------------------------------------------------------------------------
  777. enddo
  778. !------------------------------------------------------------------------
  779. ! define constants:
  780. ! a and c = constants used in evaluating universal function for
  781. ! stable case
  782. ! ca = karmen constant
  783. ! cm01 = constant part of vertical integral of universal
  784. ! function; stable case ( 0.5 < zeta < or = 10.0)
  785. ! cm02 = constant part of vertical integral of universal
  786. ! function; stable case ( zeta > 10.0)
  787. !------------------------------------------------------------------------
  788. en = 2.
  789. c = .76
  790. a = 5.
  791. ca = .4
  792. cmo1 = .5*a - 1.648
  793. cmo2 = 17.193 + .5*a - 10.*c
  794. boycon = .61
  795. rovcp=rgas/cp
  796. ! write(0,*)'rgas,cp,rovcp ', rgas,cp,rovcp
  797. ! write(0,*)'--------------------------------------------------'
  798. ! write(0,*)'pkmax, pspc, theta, zkmax, zoc'
  799. ! write(0,*)'--------------------------------------------------'
  800. do i = its,ite
  801. ! theta(i) = tkmax(i)*rqc9
  802. theta(i) = tkmax(i)/((pkmax(i)/pspc(i))**rovcp)
  803. vrtkx(i) = 1.0 + boycon*rkmax(i)
  804. ! zkmax(i) = rgas*tkmax(i)*qqlog(kmax)*og
  805. zkmax(i) = -rgas*tkmax(i)*alog(pkmax(i)/pspc(i))*og
  806. ! IF(I==78)write(0,*)I,JFIX,pkmax(i),pspc(i),theta(i),zkmax(i),zoc(i)
  807. enddo
  808. ! write(0,*)'pkmax,pspc ', pkmax,pspc
  809. ! write(0,*)'theta, zkmax, zoc ', theta, zkmax, zoc
  810. !------------------------------------------------------------------------
  811. ! get saturation mixing ratios at surface
  812. !------------------------------------------------------------------------
  813. do i = its,ite
  814. tsg (i) = tstar(i)
  815. tab1 (i) = tstar(i) - 153.16
  816. it (i) = IFIX(tab1(i))
  817. tab2 (i) = tab1(i) - FLOAT(it(i))
  818. t1 (i) = tab(it(i) + 1)
  819. t2 (i) = table(it(i) + 1)
  820. estso(i) = t1(i) + tab2(i)*t2(i)
  821. psps1 = (pss(i) - estso(i))
  822. if(psps1 .EQ. 0.0)then
  823. psps1 = .1
  824. endif
  825. rstso(i) = 0.622*estso(i)/psps1
  826. vrts (i) = 1. + boycon*ecof(i)*rstso(i)
  827. enddo
  828. !------------------------------------------------------------------------
  829. ! check if consideration of virtual temperature changes stability.
  830. ! if so, set "dthetav" to near neutral value (1.0e-4). also check
  831. ! for very small lapse rates; if ABS(tempa1) <1.0e-4 then
  832. ! tempa1=1.0e-4
  833. !------------------------------------------------------------------------
  834. do i = its,ite
  835. tempa1(i) = theta(i)*vrtkx(i) - tstar(i)*vrts(i)
  836. tempa2(i) = tempa1(i)*(theta(i) - tstar(i))
  837. if (tempa2(i) .LT. 0.) tempa1(i) = 1.0e-4
  838. tab1(i) = ABS(tempa1(i))
  839. if (tab1(i) .LT. 1.0e-4) tempa1(i) = 1.0e-4
  840. !------------------------------------------------------------------------
  841. ! compute bulk richardson number "rib" at each point. if "rib"
  842. ! exceeds 95% of critical richardson number "tab1" then "rib = tab1"
  843. !------------------------------------------------------------------------
  844. rib (i) = g*zkmax(i)*tempa1(i)/ &
  845. (tkmax(i)*vrtkx(i)*wind(i)*wind(i))
  846. tab2(i) = ABS(zoc(i))
  847. tab1(i) = 0.95/(c*(1. - tab2(i)/zkmax(i)))
  848. if (rib(i) .GT. tab1(i)) rib(i) = tab1(i)
  849. enddo
  850. do i = its,ite
  851. zeta(i) = ca*rib(i)/0.03
  852. enddo
  853. ! write(0,*)'rib,zeta,vrtkx,vrts(ite) ', rib(ite),zeta(ite), &
  854. ! vrtkx(ite),vrts(ite)
  855. !------------------------------------------------------------------------
  856. ! begin looping through points on line, solving wegsteins iteration
  857. ! for zeta at each point, and using hicks functions
  858. !------------------------------------------------------------------------
  859. !------------------------------------------------------------------------
  860. ! set initial guess of zeta=non - dimensional height "szeta" for
  861. ! stable points
  862. !------------------------------------------------------------------------
  863. rca = 1./ca
  864. enrca = en*rca
  865. ! turn off interfacial layer by zeroing out enrca
  866. enrca = 0.0
  867. zog = .0185*og
  868. !------------------------------------------------------------------------
  869. ! stable points
  870. !------------------------------------------------------------------------
  871. ip = 0
  872. do i = its,ite
  873. if (zeta(i) .GE. 0.0) then
  874. ip = ip + 1
  875. istb(ip) = i
  876. endif
  877. enddo
  878. if (ip .EQ. 0) go to 170
  879. do i = 1,ip
  880. szetam(i) = 1.0e+30
  881. sgzm(i) = 0.0e+00
  882. szeta(i) = zeta(istb(i))
  883. ifz(i) = 1
  884. enddo
  885. !------------------------------------------------------------------------
  886. ! begin wegstein iteration for "zeta" at stable points using
  887. ! hicks(1976)
  888. !------------------------------------------------------------------------
  889. do icnt = 1,icntx
  890. do i = 1,ip
  891. if (ifz(i) .EQ. 0) go to 80
  892. zal1g = ALOG(szeta(i))
  893. if (szeta(i) .LE. 0.5) then
  894. fmz1(i) = (zal1g + a*szeta(i))*rca
  895. else if (szeta(i) .GT. 0.5 .AND. szeta(i) .LE. 10.) then
  896. rzeta1 = 1./szeta(i)
  897. fmz1(i) = (8.*zal1g + 4.25*rzeta1 - &
  898. 0.5*rzeta1*rzeta1 + cmo1)*rca
  899. else if (szeta(i) .GT. 10.) then
  900. fmz1(i) = (c*szeta(i) + cmo2)*rca
  901. endif
  902. szetao = ABS(zoc(istb(i)))/zkmax(istb(i))*szeta(i)
  903. zal2g = ALOG(szetao)
  904. if (szetao .LE. 0.5) then
  905. fmzo1(i) = (zal2g + a*szetao)*rca
  906. sfzo (i) = 1. + a*szetao
  907. else if (szetao .GT. 0.5 .AND. szetao .LE. 10.) then
  908. rzeta2 = 1./szetao
  909. fmzo1(i) = (8.*zal2g + 4.25*rzeta2 - &
  910. 0.5*rzeta2*rzeta2 + cmo1)*rca
  911. sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
  912. else if (szetao .GT. 10.) then
  913. fmzo1(i) = (c*szetao + cmo2)*rca
  914. sfzo (i) = c*szetao
  915. endif
  916. ! compute heat & moisture parts of zot.. for calculation of sfh
  917. szetho = ABS(zot(istb(i)))/zkmax(istb(i))*szeta(i)
  918. zal2gh = ALOG(szetho)
  919. if (szetho .LE. 0.5) then
  920. fhzo1(i) = (zal2gh + a*szetho)*rca
  921. sfzo (i) = 1. + a*szetho
  922. else if (szetho .GT. 0.5 .AND. szetho .LE. 10.) then
  923. rzeta2 = 1./szetho
  924. fhzo1(i) = (8.*zal2gh + 4.25*rzeta2 - &
  925. 0.5*rzeta2*rzeta2 + cmo1)*rca
  926. sfzo (i) = 8.0 - 4.25*rzeta2 + rzeta2*rzeta2
  927. else if (szetho .GT. 10.) then
  928. fhzo1(i) = (c*szetho + cmo2)*rca
  929. sfzo (i) = c*szetho
  930. endif
  931. !------------------------------------------------------------------------
  932. ! compute universal function at 10 meters for diagnostic purposes
  933. !------------------------------------------------------------------------
  934. !!!! if (ngd .EQ. nNEST) then
  935. szet10 = ABS(z10(istb(i)))/zkmax(istb(i))*szeta(i)
  936. zal10 = ALOG(szet10)
  937. if (szet10 .LE. 0.5) then
  938. fmz10(i) = (zal10 + a*szet10)*rca
  939. else if (szet10 .GT. 0.5 .AND. szet10 .LE. 10.) then
  940. rzeta2 = 1./szet10
  941. fmz10(i) = (8.*zal10 + 4.25*rzeta2 - &
  942. 0.5*rzeta2*rzeta2 + cmo1)*rca
  943. else if (szet10 .GT. 10.) then
  944. fmz10(i) = (c*szet10 + cmo2)*rca
  945. endif
  946. sf10(i) = fmz10(i) - fmzo1(i)
  947. ! compute 2m values for diagnostics in HWRF
  948. szet2 = ABS(z2 (istb(i)))/zkmax(istb(i))*szeta(i)
  949. zal2 = ALOG(szet2 )
  950. if (szet2 .LE. 0.5) then
  951. fmz2 (i) = (zal2 + a*szet2 )*rca
  952. else if (szet2 .GT. 0.5 .AND. szet2 .LE. 2.) then
  953. rzeta2 = 1./szet2
  954. fmz2 (i) = (8.*zal2 + 4.25*rzeta2 - &
  955. 0.5*rzeta2*rzeta2 + cmo1)*rca
  956. else if (szet2 .GT. 2.) then
  957. fmz2 (i) = (c*szet2 + cmo2)*rca
  958. endif
  959. sf2 (i) = fmz2 (i) - fmzo1(i)
  960. !!!! endif
  961. sfm(i) = fmz1(i) - fmzo1(i)
  962. sfh(i) = fmz1(i) - fhzo1(i)
  963. sgz = ca*rib(istb(i))*sfm(i)*sfm(i)/ &
  964. (sfh(i) + enrca*sfzo(i))
  965. fmz = (sgz - szeta(i))/szeta(i)
  966. fmzo = ABS(fmz)
  967. if (fmzo .GE. 5.0e-5) then
  968. sq = (sgz - sgzm(i))/(szeta(i) - szetam(i))
  969. if(sq .EQ. 1) then
  970. write(0,*)'NCO ERROR DIVIDE BY ZERO IN MFLUX2 (STABLE CASE)'
  971. write(0,*)'sq is 1 ',fmzo,sgz,sgzm(i),szeta(i),szetam(i)
  972. endif
  973. szetam(i) = szeta(i)
  974. szeta (i) = (sgz - szeta(i)*sq)/(1.0 - sq)
  975. sgzm (i) = sgz
  976. else
  977. ifz(i) = 0
  978. endif
  979. 80 continue
  980. enddo
  981. enddo
  982. do i = 1,ip
  983. if (ifz(i) .GE. 1) go to 110
  984. enddo
  985. go to 130
  986. 110 continue
  987. write(6,120)
  988. 120 format(2X, ' NON-CONVERGENCE FOR STABLE ZETA IN ROW ')
  989. ! call MPI_CLOSE(1,routine)
  990. !------------------------------------------------------------------------
  991. ! update "zo" for ocean points. "zo"cannot be updated within the
  992. ! wegsteins iteration as the scheme (for the near neutral case)
  993. ! can become unstable
  994. !------------------------------------------------------------------------
  995. 130 continue
  996. do i = 1,ip
  997. szo = zoc(istb(i))
  998. if (szo .LT. 0.0) then
  999. wndm=wind(istb(i))*0.01
  1000. if(wndm.lt.15.0) then
  1001. ckg=0.0185*og
  1002. else
  1003. !! ckg=(0.000308*wndm+0.00925)*og
  1004. !! ckg=(0.000616*wndm)*og
  1005. ckg=(sfenth*(4*0.000308*wndm) + (1.-sfenth)*0.0185 )*og
  1006. endif
  1007. szo = - ckg*wind(istb(i))*wind(istb(i))/ &
  1008. (sfm(i)*sfm(i))
  1009. cons_p000001 = .000001
  1010. cons_7 = 7.
  1011. vis = 1.4E-1
  1012. ustar = sqrt( -szo / zog)
  1013. restar = -ustar * szo / vis
  1014. restar = max(restar,cons_p000001)
  1015. ! Rat taken from Zeng, Zhao and Dickinson 1997
  1016. rat = 2.67 * restar ** .25 - 2.57
  1017. rat = min(rat ,cons_7) !constant
  1018. rat=0.
  1019. zot(istb(i)) = szo * exp(-rat)
  1020. else
  1021. zot(istb(i)) = zoc(istb(i))
  1022. endif
  1023. ! in hwrf thermal zn

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