PageRenderTime 47ms CodeModel.GetById 18ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/share/module_soil_pre.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 4250 lines | 3067 code | 692 blank | 491 comment | 30 complexity | 9a91b75be6fa25ff4e561715c40b39be MD5 | raw file
Possible License(s): AGPL-1.0
  1. #if ( ! NMM_CORE == 1 )
  2. MODULE module_soil_pre
  3. USE module_date_time
  4. USE module_state_description
  5. CHARACTER (LEN=3) :: num_cat_count
  6. INTEGER , PARAMETER , DIMENSION(0:300) :: ints = &
  7. (/ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, &
  8. 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, &
  9. 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, &
  10. 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, &
  11. 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, &
  12. 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, &
  13. 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, &
  14. 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, &
  15. 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, &
  16. 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, &
  17. 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, &
  18. 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, &
  19. 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, &
  20. 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, &
  21. 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, &
  22. 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, &
  23. 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, &
  24. 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, &
  25. 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, &
  26. 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, &
  27. 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, &
  28. 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, &
  29. 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, &
  30. 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, &
  31. 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, &
  32. 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, &
  33. 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, &
  34. 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, &
  35. 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, &
  36. 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300 /)
  37. ! Excluded middle processing
  38. LOGICAL , SAVE :: hold_ups
  39. INTEGER , SAVE :: em_width
  40. LOGICAL , EXTERNAL :: skip_middle_points_t
  41. CONTAINS
  42. SUBROUTINE adjust_for_seaice_pre ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
  43. xland , landusef , isltyp , soilcat , soilctop , &
  44. soilcbot , tmn , &
  45. seaice_threshold , &
  46. fractional_seaice, &
  47. num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
  48. iswater , isice , &
  49. scheme , &
  50. ids , ide , jds , jde , kds , kde , &
  51. ims , ime , jms , jme , kms , kme , &
  52. its , ite , jts , jte , kts , kte )
  53. IMPLICIT NONE
  54. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  55. ims , ime , jms , jme , kms , kme , &
  56. its , ite , jts , jte , kts , kte , &
  57. iswater , isice
  58. INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
  59. REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
  60. REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
  61. REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
  62. INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
  63. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
  64. vegcat, xland , soilcat , tmn
  65. REAL , INTENT(IN) :: seaice_threshold
  66. INTEGER :: i , j , num_seaice_changes , loop
  67. CHARACTER (LEN=132) :: message
  68. INTEGER, INTENT(IN) :: fractional_seaice
  69. REAL :: XICE_THRESHOLD
  70. IF ( FRACTIONAL_SEAICE == 0 ) THEN
  71. xice_threshold = 0.5
  72. ELSEIF ( FRACTIONAL_SEAICE == 1 ) THEN
  73. xice_threshold = 0.02
  74. ENDIF
  75. num_seaice_changes = 0
  76. fix_seaice : SELECT CASE ( scheme )
  77. CASE ( SLABSCHEME )
  78. DO j = jts , MIN(jde-1,jte)
  79. DO i = its , MIN(ide-1,ite)
  80. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  81. IF ( xice(i,j) .GT. 200.0 ) THEN
  82. xice(i,j) = 0.
  83. num_seaice_changes = num_seaice_changes + 1
  84. END IF
  85. END DO
  86. END DO
  87. IF ( num_seaice_changes .GT. 0 ) THEN
  88. WRITE ( message , FMT='(A,I6)' ) &
  89. 'Total pre number of sea ice locations removed (due to FLAG values) = ', &
  90. num_seaice_changes
  91. CALL wrf_debug ( 0 , message )
  92. END IF
  93. num_seaice_changes = 0
  94. DO j = jts , MIN(jde-1,jte)
  95. DO i = its , MIN(ide-1,ite)
  96. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  97. IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
  98. ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
  99. IF ( FRACTIONAL_SEAICE == 0 ) THEN
  100. xice(i,j) = 1.0
  101. ENDIF
  102. num_seaice_changes = num_seaice_changes + 1
  103. if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
  104. vegcat(i,j)=isice
  105. ivgtyp(i,j)=isice
  106. lu_index(i,j)=isice
  107. landmask(i,j)=1.
  108. xland(i,j)=1.
  109. DO loop=1,num_veg_cat
  110. landusef(i,loop,j)=0.
  111. END DO
  112. landusef(i,ivgtyp(i,j),j)=1.
  113. isltyp(i,j) = 16
  114. soilcat(i,j)=isltyp(i,j)
  115. DO loop=1,num_soil_top_cat
  116. soilctop(i,loop,j)=0
  117. END DO
  118. DO loop=1,num_soil_bot_cat
  119. soilcbot(i,loop,j)=0
  120. END DO
  121. soilctop(i,isltyp(i,j),j)=1.
  122. soilcbot(i,isltyp(i,j),j)=1.
  123. ELSE
  124. xice(i,j) = 0.0
  125. END IF
  126. END DO
  127. END DO
  128. IF ( num_seaice_changes .GT. 0 ) THEN
  129. WRITE ( message , FMT='(A,I6)' ) &
  130. 'Total pre number of sea ice location changes (water to land) = ', num_seaice_changes
  131. CALL wrf_debug ( 0 , message )
  132. END IF
  133. CASE ( LSMSCHEME , NOAHMPSCHEME , RUCLSMSCHEME , SSIBSCHEME) !mchen add for ssib
  134. num_seaice_changes = 0
  135. DO j = jts , MIN(jde-1,jte)
  136. DO i = its , MIN(ide-1,ite)
  137. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  138. IF ( landmask(i,j) .GT. 0.5 ) THEN
  139. if (xice(i,j).gt.0) num_seaice_changes = num_seaice_changes + 1
  140. xice(i,j) = 0.
  141. END IF
  142. END DO
  143. END DO
  144. IF ( num_seaice_changes .GT. 0 ) THEN
  145. WRITE ( message , FMT='(A,I6)' ) &
  146. 'Total pre number of land location changes (seaice set to zero) = ', num_seaice_changes
  147. CALL wrf_debug ( 0 , message )
  148. END IF
  149. END SELECT fix_seaice
  150. END SUBROUTINE adjust_for_seaice_pre
  151. SUBROUTINE adjust_for_seaice_post ( xice , landmask , tsk_old , tsk , ivgtyp , vegcat , lu_index , &
  152. xland , landusef , isltyp , soilcat , soilctop , &
  153. soilcbot , tmn , vegfra , &
  154. tslb , smois , sh2o , &
  155. seaice_threshold , &
  156. sst , flag_sst , &
  157. fractional_seaice, &
  158. num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
  159. num_soil_layers , &
  160. iswater , isice , &
  161. scheme , &
  162. ids , ide , jds , jde , kds , kde , &
  163. ims , ime , jms , jme , kms , kme , &
  164. its , ite , jts , jte , kts , kte )
  165. IMPLICIT NONE
  166. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  167. ims , ime , jms , jme , kms , kme , &
  168. its , ite , jts , jte , kts , kte , &
  169. iswater , isice
  170. INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
  171. INTEGER , INTENT(IN) :: num_soil_layers
  172. REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
  173. REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
  174. REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
  175. REAL , DIMENSION(ims:ime,1:num_soil_layers,jms:jme) , INTENT(INOUT):: tslb , smois , sh2o
  176. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN):: sst
  177. INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
  178. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
  179. vegcat, xland , soilcat , tmn , &
  180. tsk_old , vegfra
  181. INTEGER , INTENT(IN) :: flag_sst
  182. REAL , INTENT(IN) :: seaice_threshold
  183. REAL :: total_depth , mid_point_depth
  184. INTEGER :: i , j , num_seaice_changes , loop
  185. CHARACTER (LEN=132) :: message
  186. INTEGER, INTENT(IN) :: fractional_seaice
  187. real :: xice_threshold
  188. IF ( FRACTIONAL_SEAICE == 0 ) THEN
  189. xice_threshold = 0.5
  190. ELSEIF ( FRACTIONAL_SEAICE == 1 ) THEN
  191. xice_threshold = 0.02
  192. ENDIF
  193. num_seaice_changes = 0
  194. fix_seaice : SELECT CASE ( scheme )
  195. CASE ( SLABSCHEME )
  196. CASE ( LSMSCHEME , NOAHMPSCHEME , SSIBSCHEME ) !mchen add for ssib
  197. DO j = jts , MIN(jde-1,jte)
  198. DO i = its , MIN(ide-1,ite)
  199. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  200. IF ( xice(i,j) .GT. 200.0 ) THEN
  201. xice(i,j) = 0.
  202. num_seaice_changes = num_seaice_changes + 1
  203. END IF
  204. END DO
  205. END DO
  206. IF ( num_seaice_changes .GT. 0 ) THEN
  207. WRITE ( message , FMT='(A,I6)' ) &
  208. 'Total post number of sea ice locations removed (due to FLAG values) = ', &
  209. num_seaice_changes
  210. CALL wrf_debug ( 0 , message )
  211. END IF
  212. num_seaice_changes = 0
  213. DO j = jts , MIN(jde-1,jte)
  214. DO i = its , MIN(ide-1,ite)
  215. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  216. IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
  217. ( ( tsk_old(i,j) .GT. 170 ) .AND. ( tsk_old(i,j) .LT. 400 ) ) )THEN
  218. tsk(i,j) = tsk_old(i,j)
  219. END IF
  220. IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
  221. ( ( tsk_old(i,j) .LT. 170 ) .OR. ( tsk_old(i,j) .GT. 400 ) ) )THEN
  222. print *,'TSK woes in seaice post, i,j=',i,j,' tsk = ',tsk(i,j), tsk_old(i,j)
  223. CALL wrf_error_fatal ( 'TSK is unrealistic, problems for seaice post')
  224. ELSE IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
  225. ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
  226. IF ( FRACTIONAL_SEAICE == 0 ) THEN
  227. xice(i,j) = 1.0
  228. ENDIF
  229. num_seaice_changes = num_seaice_changes + 1
  230. if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
  231. vegcat(i,j)=isice
  232. ivgtyp(i,j)=isice
  233. lu_index(i,j)=isice
  234. landmask(i,j)=1.
  235. xland(i,j)=1.
  236. vegfra(i,j)=0.
  237. DO loop=1,num_veg_cat
  238. landusef(i,loop,j)=0.
  239. END DO
  240. landusef(i,ivgtyp(i,j),j)=1.
  241. tsk_old(i,j) = tsk(i,j)
  242. isltyp(i,j) = 16
  243. soilcat(i,j)=isltyp(i,j)
  244. DO loop=1,num_soil_top_cat
  245. soilctop(i,loop,j)=0
  246. END DO
  247. DO loop=1,num_soil_bot_cat
  248. soilcbot(i,loop,j)=0
  249. END DO
  250. soilctop(i,isltyp(i,j),j)=1.
  251. soilcbot(i,isltyp(i,j),j)=1.
  252. total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
  253. DO loop = 1,num_soil_layers
  254. mid_point_depth=(total_depth/num_soil_layers)/2. + &
  255. (loop-1)*(total_depth/num_soil_layers)
  256. tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
  257. mid_point_depth*tmn(i,j) ) / total_depth
  258. END DO
  259. DO loop=1,num_soil_layers
  260. smois(i,loop,j) = 1.0
  261. sh2o(i,loop,j) = 0.0
  262. END DO
  263. ELSE IF ( xice(i,j) .LT. xice_threshold ) THEN
  264. xice(i,j) = 0.
  265. END IF
  266. END DO
  267. END DO
  268. IF ( num_seaice_changes .GT. 0 ) THEN
  269. WRITE ( message , FMT='(A,I6)' ) &
  270. 'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
  271. CALL wrf_debug ( 0 , message )
  272. END IF
  273. CASE ( RUCLSMSCHEME )
  274. DO j = jts , MIN(jde-1,jte)
  275. DO i = its , MIN(ide-1,ite)
  276. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  277. IF ( xice(i,j) .GT. 200.0 ) THEN
  278. xice(i,j) = 0.
  279. num_seaice_changes = num_seaice_changes + 1
  280. END IF
  281. END DO
  282. END DO
  283. IF ( num_seaice_changes .GT. 0 ) THEN
  284. WRITE ( message , FMT='(A,I6)' ) &
  285. 'Total post number of sea ice locations removed (due to FLAG values) = ', &
  286. num_seaice_changes
  287. CALL wrf_debug ( 0 , message )
  288. END IF
  289. num_seaice_changes = 0
  290. DO j = jts , MIN(jde-1,jte)
  291. DO i = its , MIN(ide-1,ite)
  292. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  293. IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
  294. ( ( tsk_old(i,j) .GT. 170 ) .AND. ( tsk_old(i,j) .LT. 400 ) ) )THEN
  295. tsk(i,j) = tsk_old(i,j)
  296. END IF
  297. IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
  298. ( ( tsk_old(i,j) .LT. 170 ) .OR. ( tsk_old(i,j) .GT. 400 ) ) )THEN
  299. print *,'TSK woes in seaice post, i,j=',i,j,' tsk = ',tsk(i,j), tsk_old(i,j)
  300. CALL wrf_error_fatal ( 'TSK is unrealistic, problems for seaice post')
  301. ELSE IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
  302. ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
  303. IF ( FRACTIONAL_SEAICE == 0 ) THEN
  304. xice(i,j) = 1.0
  305. ELSE
  306. xice(i,j)=max(0.25,xice(i,j))
  307. ENDIF
  308. num_seaice_changes = num_seaice_changes + 1
  309. if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
  310. vegcat(i,j)=isice
  311. ivgtyp(i,j)=isice
  312. lu_index(i,j)=isice
  313. landmask(i,j)=1.
  314. xland(i,j)=1.
  315. vegfra(i,j)=0.
  316. DO loop=1,num_veg_cat
  317. landusef(i,loop,j)=0.
  318. END DO
  319. landusef(i,ivgtyp(i,j),j)=1.
  320. !tgs - compute blended sea ice/water skin temperature
  321. if(flag_sst.eq.1) then
  322. tsk(i,j) = xice(i,j)*(min(seaice_threshold,tsk(i,j))) &
  323. +(1-xice(i,j))*sst(i,j)
  324. else
  325. tsk(i,j) = xice(i,j)*(min(seaice_threshold,tsk(i,j))) &
  326. +(1-xice(i,j))*tsk(i,j)
  327. endif
  328. tsk_old(i,j) = tsk(i,j)
  329. isltyp(i,j) = 16
  330. soilcat(i,j)=isltyp(i,j)
  331. DO loop=1,num_soil_top_cat
  332. soilctop(i,loop,j)=0
  333. END DO
  334. DO loop=1,num_soil_bot_cat
  335. soilcbot(i,loop,j)=0
  336. END DO
  337. soilctop(i,isltyp(i,j),j)=1.
  338. soilcbot(i,isltyp(i,j),j)=1.
  339. total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
  340. tslb(i,1,j) = tsk(i,j)
  341. tslb(i,num_soil_layers,j) = tmn(i,j)
  342. DO loop = 2,num_soil_layers-1
  343. mid_point_depth=(total_depth/num_soil_layers)/4. + &
  344. (loop-2)*(total_depth/num_soil_layers)
  345. tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
  346. mid_point_depth*tmn(i,j) ) / total_depth
  347. END DO
  348. DO loop=1,num_soil_layers
  349. smois(i,loop,j) = 1.0
  350. sh2o(i,loop,j) = 0.0
  351. END DO
  352. ELSE IF ( xice(i,j) .LT. xice_threshold ) THEN
  353. xice(i,j) = 0.
  354. END IF
  355. END DO
  356. END DO
  357. IF ( num_seaice_changes .GT. 0 ) THEN
  358. WRITE ( message , FMT='(A,I6)' ) &
  359. 'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
  360. CALL wrf_debug ( 0 , message )
  361. END IF
  362. END SELECT fix_seaice
  363. END SUBROUTINE adjust_for_seaice_post
  364. SUBROUTINE process_percent_cat_new ( landmask , &
  365. landuse_frac , soil_top_cat , soil_bot_cat , &
  366. isltyp , ivgtyp , &
  367. num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
  368. ids , ide , jds , jde , kds , kde , &
  369. ims , ime , jms , jme , kms , kme , &
  370. its , ite , jts , jte , kts , kte , &
  371. iswater )
  372. IMPLICIT NONE
  373. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  374. ims , ime , jms , jme , kms , kme , &
  375. its , ite , jts , jte , kts , kte , &
  376. iswater
  377. INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
  378. REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landuse_frac
  379. REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(IN):: soil_top_cat
  380. REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(IN):: soil_bot_cat
  381. INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
  382. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask
  383. INTEGER :: i , j , l , ll, dominant_index
  384. REAL :: dominant_value
  385. #ifdef WRF_CHEM
  386. ! REAL :: lwthresh = .99
  387. REAL :: lwthresh = .50
  388. #else
  389. REAL :: lwthresh = .50
  390. #endif
  391. INTEGER , PARAMETER :: iswater_soil = 14
  392. INTEGER :: iforce
  393. CHARACTER (LEN=132) :: message
  394. CHARACTER(LEN=256) :: mminlu
  395. LOGICAL :: aggregate_lu
  396. integer :: change_water , change_land
  397. change_water = 0
  398. change_land = 0
  399. ! Sanity check on the 50/50 points
  400. DO j = jts , MIN(jde-1,jte)
  401. DO i = its , MIN(ide-1,ite)
  402. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  403. dominant_value = landuse_frac(i,iswater,j)
  404. IF ( dominant_value .EQ. lwthresh ) THEN
  405. DO l = 1 , num_veg_cat
  406. IF ( l .EQ. iswater ) CYCLE
  407. IF ( ( landuse_frac(i,l,j) .EQ. lwthresh ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
  408. PRINT *,i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
  409. landuse_frac(i,l,j) = lwthresh - .01
  410. landuse_frac(i,iswater,j) = lwthresh + 0.01
  411. ELSE IF ( ( landuse_frac(i,l,j) .EQ. lwthresh ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
  412. PRINT *,i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
  413. landuse_frac(i,l,j) = lwthresh + .01
  414. landuse_frac(i,iswater,j) = lwthresh - 0.01
  415. END IF
  416. END DO
  417. END IF
  418. END DO
  419. END DO
  420. ! Compute the aggregate of the vegetation/land use categories. Lump all of the grasses together,
  421. ! all of the shrubs, all of the trees, etc. Choose the correct set of available land use
  422. ! categories. Also, make sure that the user wants to actually avail themselves of aforementioned
  423. ! opportunity, as mayhaps they don't.
  424. CALL nl_get_mminlu ( 1 , mminlu )
  425. CALL nl_get_aggregate_lu ( 1 , aggregate_lu )
  426. IF ( aggregate_lu ) THEN
  427. DO j = jts , MIN(jde-1,jte)
  428. DO i = its , MIN(ide-1,ite)
  429. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  430. CALL aggregate_categories_part1 ( landuse_frac , iswater , num_veg_cat , mminlu(1:4) )
  431. END DO
  432. END DO
  433. END IF
  434. ! Compute the dominant VEGETATION INDEX.
  435. DO j = jts , MIN(jde-1,jte)
  436. DO i = its , MIN(ide-1,ite)
  437. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  438. dominant_value = landuse_frac(i,1,j)
  439. dominant_index = 1
  440. DO l = 2 , num_veg_cat
  441. IF ( l .EQ. iswater ) THEN
  442. ! wait a bit
  443. ELSE IF ( ( l .NE. iswater ) .AND. ( landuse_frac(i,l,j) .GT. dominant_value ) ) THEN
  444. dominant_value = landuse_frac(i,l,j)
  445. dominant_index = l
  446. END IF
  447. END DO
  448. IF ( landuse_frac(i,iswater,j) .GT. lwthresh ) THEN
  449. dominant_value = landuse_frac(i,iswater,j)
  450. dominant_index = iswater
  451. ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
  452. ( landmask(i,j) .LT. 0.5) .AND. &
  453. ( dominant_value .EQ. lwthresh) ) THEN
  454. dominant_value = landuse_frac(i,iswater,j)
  455. dominant_index = iswater
  456. ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
  457. ( landmask(i,j) .GT. 0.5) .AND. &
  458. ( dominant_value .EQ. lwthresh) ) THEN
  459. !no op
  460. ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh ) .AND. &
  461. ( dominant_value .LT. lwthresh ) ) THEN
  462. dominant_value = landuse_frac(i,iswater,j)
  463. dominant_index = iswater
  464. END IF
  465. IF ( dominant_index .EQ. iswater ) THEN
  466. if(landmask(i,j).gt.lwthresh) then
  467. !print *,'changing to water at point ',i,j
  468. !WRITE ( num_cat_count , FMT = '(I3)' ) num_veg_cat
  469. !WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) ints(1:num_veg_cat)
  470. !CALL wrf_debug(1,message)
  471. !WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) nint(landuse_frac(i,:,j)*100)
  472. !CALL wrf_debug(1,message)
  473. change_water=change_water+1
  474. endif
  475. landmask(i,j) = 0
  476. ELSE IF ( dominant_index .NE. iswater ) THEN
  477. if(landmask(i,j).lt.lwthresh) then
  478. !print *,'changing to land at point ',i,j
  479. !WRITE ( num_cat_count , FMT = '(I3)' ) num_veg_cat
  480. !WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) ints(1:num_veg_cat)
  481. !CALL wrf_debug(1,message)
  482. !WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) nint(landuse_frac(i,:,j)*100)
  483. !CALL wrf_debug(1,message)
  484. change_land=change_land+1
  485. endif
  486. landmask(i,j) = 1
  487. END IF
  488. ivgtyp(i,j) = dominant_index
  489. END DO
  490. END DO
  491. ! Compute the dominant SOIL TEXTURE INDEX, TOP.
  492. iforce = 0
  493. DO i = its , MIN(ide-1,ite)
  494. DO j = jts , MIN(jde-1,jte)
  495. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  496. dominant_value = soil_top_cat(i,1,j)
  497. dominant_index = 1
  498. IF ( landmask(i,j) .GT. lwthresh ) THEN
  499. DO l = 2 , num_soil_top_cat
  500. IF ( ( l .NE. iswater_soil ) .AND. ( soil_top_cat(i,l,j) .GT. dominant_value ) ) THEN
  501. dominant_value = soil_top_cat(i,l,j)
  502. dominant_index = l
  503. END IF
  504. END DO
  505. IF ( dominant_value .LT. 0.01 ) THEN
  506. iforce = iforce + 1
  507. WRITE ( message , FMT = '(A,I4,I4)' ) &
  508. 'based on landuse, changing soil to land at point ',i,j
  509. CALL wrf_debug(1,message)
  510. WRITE ( num_cat_count , FMT = '(I3)' ) num_soil_top_cat
  511. WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) (ints(l),l=1,num_soil_top_cat)
  512. CALL wrf_debug(1,message)
  513. WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) &
  514. ((nint(soil_top_cat(i,ints(l),j)*100)), l=1,num_soil_top_cat)
  515. CALL wrf_debug(1,message)
  516. dominant_index = 8
  517. END IF
  518. ELSE
  519. dominant_index = iswater_soil
  520. END IF
  521. isltyp(i,j) = dominant_index
  522. END DO
  523. END DO
  524. if(iforce.ne.0)then
  525. WRITE(message,FMT='(A,I4,A,I6)' ) &
  526. 'forcing artificial silty clay loam at ',iforce,' points, out of ',&
  527. (MIN(ide-1,ite)-its+1)*(MIN(jde-1,jte)-jts+1)
  528. CALL wrf_debug(0,message)
  529. endif
  530. print *,'LAND CHANGE = ',change_land
  531. print *,'WATER CHANGE = ',change_water
  532. END SUBROUTINE process_percent_cat_new
  533. SUBROUTINE process_soil_real ( tsk , tmn , tavgsfc, &
  534. landmask , sst , ht, toposoil, &
  535. st_input , sm_input , sw_input , &
  536. st_levels_input , sm_levels_input , sw_levels_input , &
  537. zs , dzs , tslb , smois , sh2o , &
  538. flag_sst , flag_tavgsfc, flag_soilhgt, &
  539. flag_soil_layers, flag_soil_levels, &
  540. ids , ide , jds , jde , kds , kde , &
  541. ims , ime , jms , jme , kms , kme , &
  542. its , ite , jts , jte , kts , kte , &
  543. sf_surface_physics , num_soil_layers , real_data_init_type , &
  544. num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  545. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
  546. IMPLICIT NONE
  547. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  548. ims , ime , jms , jme , kms , kme , &
  549. its , ite , jts , jte , kts , kte , &
  550. sf_surface_physics , num_soil_layers , real_data_init_type , &
  551. num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  552. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc
  553. INTEGER , INTENT(IN) :: flag_sst, flag_tavgsfc
  554. INTEGER , INTENT(IN) :: flag_soil_layers, flag_soil_levels, flag_soilhgt
  555. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
  556. INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
  557. INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
  558. INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
  559. REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
  560. REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
  561. REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
  562. REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
  563. REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
  564. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tavgsfc, ht, toposoil
  565. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk, tmn
  566. INTEGER :: i , j , k, l , dominant_index , num_soil_cat , num_veg_cat, closest_layer
  567. REAL :: dominant_value, closest_depth, diff_cm
  568. REAL , ALLOCATABLE , DIMENSION(:) :: depth ! Soil layer thicknesses (cm)
  569. REAL, PARAMETER :: get_temp_closest_to = 30. ! use soil temperature closest to this depth (cm)
  570. REAL, PARAMETER :: something_big = 1.e6 ! Initialize closest depth as something big (cm)
  571. INTEGER :: something_far = 1000 ! Soil array index far away
  572. CHARACTER (LEN=132) :: message
  573. ! Case statement for tmn initialization
  574. ! Need to have a reasonable default value for annual mean deeeeep soil temperature
  575. ! For sf_surface_physics = 1, we want to use close to a 30 cm value
  576. ! for the bottom level of the soil temps.
  577. ! NOTE: We are assuming that soil_layers are the same for each grid point
  578. fix_bottom_level_for_temp : SELECT CASE ( sf_surface_physics )
  579. CASE (SLABSCHEME)
  580. IF ( flag_tavgsfc .EQ. 1 ) THEN
  581. CALL wrf_debug ( 0 , 'Using average surface temperature for tmn')
  582. DO j = jts , MIN(jde-1,jte)
  583. DO i = its , MIN(ide-1,ite)
  584. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  585. tmn(i,j) = tavgsfc(i,j)
  586. END DO
  587. END DO
  588. ELSE
  589. ! Look for soil temp close to 30 cm
  590. closest_layer = something_far
  591. closest_depth = something_big
  592. DO k = 1, num_st_levels_input
  593. diff_cm = abs( st_levels_input(k) - get_temp_closest_to )
  594. IF ( diff_cm < closest_depth ) THEN
  595. closest_depth = diff_cm
  596. closest_layer = k
  597. END IF
  598. END DO
  599. IF ( closest_layer == something_far ) THEN
  600. CALL wrf_debug ( 0 , 'No soil temperature data for grid%tmn near 30 cm')
  601. CALL wrf_debug ( 0 , 'Using 1 degree static annual mean temps' )
  602. ELSE
  603. write(message, FMT='(A,F7.2,A,I3)')&
  604. 'Soil temperature closest to ',get_temp_closest_to, &
  605. ' at level ',st_levels_input(closest_layer)
  606. CALL wrf_debug ( 0 , message )
  607. DO j = jts , MIN(jde-1,jte)
  608. DO i = its , MIN(ide-1,ite)
  609. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  610. tmn(i,j) = st_input(i,closest_layer+1,j)
  611. END DO
  612. END DO
  613. END IF
  614. END IF
  615. CASE (LSMSCHEME)
  616. CASE (NOAHMPSCHEME)
  617. CASE (RUCLSMSCHEME)
  618. CASE (PXLSMSCHEME)
  619. ! When the input data from met_em is in layers, there is an additional level added to the beginning
  620. ! of the array to define the surface, which is why we add the extra value (closest_layer+1)
  621. IF ( flag_tavgsfc .EQ. 1 ) THEN
  622. CALL wrf_debug ( 0 , 'Using average surface temperature for tmn')
  623. DO j = jts , MIN(jde-1,jte)
  624. DO i = its , MIN(ide-1,ite)
  625. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  626. tmn(i,j) = tavgsfc(i,j)
  627. END DO
  628. END DO
  629. ELSE
  630. ! Look for soil temp close to 30 cm
  631. closest_layer = num_st_levels_input+1
  632. closest_depth = something_big
  633. DO k = 1, num_st_levels_input
  634. diff_cm = abs( st_levels_input(k) - get_temp_closest_to )
  635. IF ( diff_cm < closest_depth ) THEN
  636. closest_depth = diff_cm
  637. closest_layer = k
  638. END IF
  639. END DO
  640. IF ( closest_layer == num_st_levels_input + 1 ) THEN
  641. CALL wrf_debug ( 0 , 'No soil temperature data for grid%tmn near 30 cm')
  642. CALL wrf_debug ( 0 , 'Using 1 degree static annual mean temps' )
  643. ELSE
  644. write(message, FMT='(A,F7.2,A,I3)')&
  645. 'Soil temperature closest to ',get_temp_closest_to, &
  646. ' at level ',st_levels_input(closest_layer)
  647. CALL wrf_debug ( 0 , message )
  648. DO j = jts , MIN(jde-1,jte)
  649. DO i = its , MIN(ide-1,ite)
  650. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  651. tmn(i,j) = st_input(i,closest_layer+1,j)
  652. END DO
  653. END DO
  654. END IF
  655. END IF
  656. #if 0
  657. ! Loop over layers and do a weighted mean
  658. IF ( ALLOCATED ( depth ) ) DEALLOCATE ( depth )
  659. ALLOCATE ( depth(num_st_levels_input) )
  660. IF ( flag_soil_layers == 1 ) THEN
  661. DO k = num_st_levels_input, 2, -1
  662. depth(k) = st_levels_input(k) - st_levels_input(k-1)
  663. END DO
  664. depth(1) = st_levels_input(1)
  665. ELSE IF ( flag_soil_levels == 1 ) THEN
  666. DO k = 2, num_st_levels_input
  667. depth(k) = st_levels_input(k) - st_levels_input(k-1)
  668. END DO
  669. depth(1) = 0.
  670. END IF
  671. DO j = jts , MIN(jde-1,jte)
  672. DO i = its , MIN(ide-1,ite)
  673. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  674. tmn(i,j) = 0.
  675. DO k = 1, num_st_levels_input
  676. tmn(i,j) = tmn(i,j) + depth(k) * st_input(i,k,j)
  677. END DO
  678. END DO
  679. END DO
  680. DEALLOCATE ( depth )
  681. #endif
  682. END SELECT fix_bottom_level_for_temp
  683. ! Adjust the various soil temperature values depending on the difference in
  684. ! elevation between the current model's elevation and the incoming data's
  685. ! orography.
  686. adjust_soil : SELECT CASE ( sf_surface_physics )
  687. CASE ( SLABSCHEME , LSMSCHEME , NOAHMPSCHEME , RUCLSMSCHEME, PXLSMSCHEME, SSIBSCHEME )
  688. CALL adjust_soil_temp_new ( tmn , sf_surface_physics , tsk , ht , &
  689. toposoil , landmask , st_input, st_levels_input, &
  690. flag_soilhgt , flag_tavgsfc , &
  691. flag_soil_layers , flag_soil_levels, &
  692. num_st_levels_input, num_st_levels_alloc, &
  693. ids , ide , jds , jde , kds , kde , &
  694. ims , ime , jms , jme , kms , kme , &
  695. its , ite , jts , jte , kts , kte )
  696. END SELECT adjust_soil
  697. ! Initialize the soil depth, and the soil temperature and moisture.
  698. IF ( ( sf_surface_physics .EQ. SLABSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  699. CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
  700. CALL init_soil_1_real ( tsk , tmn , tslb , zs , dzs , num_soil_layers , real_data_init_type , &
  701. landmask , sst , flag_sst , &
  702. ids , ide , jds , jde , kds , kde , &
  703. ims , ime , jms , jme , kms , kme , &
  704. its , ite , jts , jte , kts , kte )
  705. ELSE IF ( ( sf_surface_physics .EQ. LSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  706. CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
  707. CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
  708. st_input , sm_input , sw_input , landmask , sst , &
  709. zs , dzs , &
  710. st_levels_input , sm_levels_input , sw_levels_input , &
  711. num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  712. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
  713. flag_sst , flag_soil_layers , flag_soil_levels , &
  714. ids , ide , jds , jde , kds , kde , &
  715. ims , ime , jms , jme , kms , kme , &
  716. its , ite , jts , jte , kts , kte )
  717. ELSE IF ( ( sf_surface_physics .EQ. NOAHMPSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  718. CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
  719. CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
  720. st_input , sm_input , sw_input , landmask , sst , &
  721. zs , dzs , &
  722. st_levels_input , sm_levels_input , sw_levels_input , &
  723. num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  724. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
  725. flag_sst , flag_soil_layers , flag_soil_levels , &
  726. ids , ide , jds , jde , kds , kde , &
  727. ims , ime , jms , jme , kms , kme , &
  728. its , ite , jts , jte , kts , kte )
  729. ELSE IF ( ( sf_surface_physics .EQ. RUCLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  730. CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
  731. CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
  732. st_input , sm_input , landmask , sst , &
  733. zs , dzs , &
  734. st_levels_input , sm_levels_input , &
  735. num_soil_layers , num_st_levels_input , num_sm_levels_input , &
  736. num_st_levels_alloc , num_sm_levels_alloc , &
  737. flag_sst , flag_soil_layers , flag_soil_levels , &
  738. ids , ide , jds , jde , kds , kde , &
  739. ims , ime , jms , jme , kms , kme , &
  740. its , ite , jts , jte , kts , kte )
  741. ELSE IF ( ( sf_surface_physics .EQ. PXLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  742. CALL init_soil_depth_7 ( zs , dzs , num_soil_layers )
  743. CALL init_soil_7_real ( tsk , tmn , smois , sh2o, tslb , &
  744. st_input , sm_input , sw_input, landmask , sst , &
  745. zs , dzs , &
  746. st_levels_input , sm_levels_input , sw_levels_input, &
  747. num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  748. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
  749. flag_sst , flag_soil_layers , flag_soil_levels , &
  750. ids , ide , jds , jde , kds , kde , &
  751. ims , ime , jms , jme , kms , kme , &
  752. its , ite , jts , jte , kts , kte )
  753. ELSE IF ( ( sf_surface_physics .EQ. 8 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  754. CALL init_soil_depth_8 ( zs , dzs , num_soil_layers )
  755. CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
  756. st_input , sm_input , landmask , sst , &
  757. zs , dzs , &
  758. st_levels_input , sm_levels_input , &
  759. num_soil_layers , num_st_levels_input , num_sm_levels_input , &
  760. num_st_levels_alloc , num_sm_levels_alloc , &
  761. flag_sst , flag_soil_layers , flag_soil_levels , &
  762. ids , ide , jds , jde , kds , kde , &
  763. ims , ime , jms , jme , kms , kme , &
  764. its , ite , jts , jte , kts , kte )
  765. END IF
  766. END SUBROUTINE process_soil_real
  767. SUBROUTINE process_soil_ideal ( xland,xice,vegfra,snow,canwat, &
  768. ivgtyp,isltyp,tslb,smois, &
  769. tsk,tmn,zs,dzs, &
  770. num_soil_layers, &
  771. sf_surface_physics , &
  772. ids,ide, jds,jde, kds,kde,&
  773. ims,ime, jms,jme, kms,kme,&
  774. its,ite, jts,jte, kts,kte )
  775. IMPLICIT NONE
  776. INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
  777. ims,ime, jms,jme, kms,kme, &
  778. its,ite, jts,jte, kts,kte
  779. INTEGER, INTENT(IN) :: num_soil_layers , sf_surface_physics
  780. REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(INOUT) :: smois, tslb
  781. REAL, DIMENSION(num_soil_layers), INTENT(OUT) :: dzs,zs
  782. REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(INOUT) :: tsk, tmn
  783. REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra
  784. INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
  785. ! Local variables.
  786. INTEGER :: itf,jtf
  787. itf=MIN(ite,ide-1)
  788. jtf=MIN(jte,jde-1)
  789. IF ( ( sf_surface_physics .EQ. SLABSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  790. CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
  791. CALL init_soil_1_ideal(tsk,tmn,tslb,xland, &
  792. ivgtyp,zs,dzs,num_soil_layers, &
  793. ids,ide, jds,jde, kds,kde, &
  794. ims,ime, jms,jme, kms,kme, &
  795. its,ite, jts,jte, kts,kte )
  796. ELSE IF ( ( sf_surface_physics .EQ. LSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  797. CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
  798. CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
  799. ivgtyp,isltyp,tslb,smois,tmn, &
  800. num_soil_layers, &
  801. ids,ide, jds,jde, kds,kde, &
  802. ims,ime, jms,jme, kms,kme, &
  803. its,ite, jts,jte, kts,kte )
  804. ELSE IF ( ( sf_surface_physics .EQ. NOAHMPSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  805. CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
  806. CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
  807. ivgtyp,isltyp,tslb,smois,tmn, &
  808. num_soil_layers, &
  809. ids,ide, jds,jde, kds,kde, &
  810. ims,ime, jms,jme, kms,kme, &
  811. its,ite, jts,jte, kts,kte )
  812. ELSE IF ( ( sf_surface_physics .EQ. RUCLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  813. CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
  814. END IF
  815. END SUBROUTINE process_soil_ideal
  816. SUBROUTINE adjust_soil_temp_new ( tmn , sf_surface_physics , tsk , ter , &
  817. toposoil , landmask , st_input , st_levels_input, &
  818. flag_toposoil , flag_tavgsfc , &
  819. flag_soil_layers , flag_soil_levels, &
  820. num_st_levels_input, num_st_levels_alloc, &
  821. ids , ide , jds , jde , kds , kde , &
  822. ims , ime , jms , jme , kms , kme , &
  823. its , ite , jts , jte , kts , kte )
  824. IMPLICIT NONE
  825. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  826. ims , ime , jms , jme , kms , kme , &
  827. its , ite , jts , jte , kts , kte
  828. INTEGER , INTENT(IN) :: num_st_levels_input, num_st_levels_alloc
  829. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: ter , toposoil , landmask
  830. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tmn , tsk
  831. REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
  832. INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(IN) :: st_levels_input
  833. INTEGER , INTENT(IN) :: sf_surface_physics , flag_toposoil , flag_tavgsfc
  834. INTEGER , INTENT(IN) :: flag_soil_layers , flag_soil_levels
  835. INTEGER :: i , j, k , st_near_sfc
  836. REAL :: soil_elev_min_val , soil_elev_max_val , soil_elev_min_dif , soil_elev_max_dif
  837. ! Adjust the annual mean temperature as if it is based on from a sea-level elevation
  838. ! if the value used is from the actual annula mean data set. If the input field to
  839. ! be used for tmn is one of the first-guess input temp fields, need to do an adjustment
  840. ! only on the diff in topo from the model terrain and the first-guess terrain.
  841. SELECT CASE ( sf_surface_physics )
  842. CASE ( LSMSCHEME , NOAHMPSCHEME )
  843. DO j = jts , MIN(jde-1,jte)
  844. DO i = its , MIN(ide-1,ite)
  845. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  846. IF (landmask(i,j) .GT. 0.5 ) THEN
  847. tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
  848. END IF
  849. END DO
  850. END DO
  851. CASE (RUCLSMSCHEME)
  852. DO j = jts , MIN(jde-1,jte)
  853. DO i = its , MIN(ide-1,ite)
  854. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  855. IF (landmask(i,j) .GT. 0.5 ) THEN
  856. tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
  857. END IF
  858. END DO
  859. END DO
  860. END SELECT
  861. ! Do we have a soil field with which to modify soil temperatures?
  862. IF ( flag_toposoil .EQ. 1 ) THEN
  863. DO j = jts , MIN(jde-1,jte)
  864. DO i = its , MIN(ide-1,ite)
  865. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  866. ! Is the toposoil field OK, or is it a subversive soil elevation field. We can tell
  867. ! usually by looking at values. Anything less than -1000 m (lower than the Dead Sea) is
  868. ! bad. Anything larger than 10 km (taller than Everest) is toast. Also, anything where
  869. ! the difference between the soil elevation and the terrain is greater than 3 km means
  870. ! that the soil data is either all zeros or that the data are inconsistent. Any of these
  871. ! three conditions is grievous enough to induce a WRF fatality. However, if they are at
  872. ! a water point, then we can safely ignore them.
  873. soil_elev_min_val = toposoil(i,j)
  874. soil_elev_max_val = toposoil(i,j)
  875. soil_elev_min_dif = ter(i,j) - toposoil(i,j)
  876. soil_elev_max_dif = ter(i,j) - toposoil(i,j)
  877. IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
  878. CYCLE
  879. ELSE IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
  880. !print *,'no soil temperature elevation adjustment, soil height too small = ',toposoil(i,j)
  881. cycle
  882. ! CALL wrf_error_fatal ( 'TOPOSOIL values have large negative values < -1000 m, unrealistic.' )
  883. ENDIF
  884. IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
  885. CYCLE
  886. ELSE IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
  887. print *,'no soil temperature elevation adjustment, soil height too high = ',toposoil(i,j)
  888. cycle
  889. CALL wrf_error_fatal ( 'TOPOSOIL values have large positive values > 10,000 m , unrealistic.' )
  890. ENDIF
  891. IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
  892. ( landmask(i,j) .LT. 0.5 ) ) THEN
  893. CYCLE
  894. ELSE IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
  895. ( landmask(i,j) .GT. 0.5 ) ) THEN
  896. print *,'no soil temperature elevation adjustment, diff of soil height and terrain = ',ter(i,j) - toposoil(i,j)
  897. cycle
  898. CALL wrf_error_fatal ( 'TOPOSOIL difference with terrain elevation differs by more than 3000 m, unrealistic' )
  899. ENDIF
  900. ! For each of the fields that we would like to modify, check to see if it came in from the SI.
  901. ! If so, then use a -6.5 K/km lapse rate (based on the elevation diffs). We only adjust when we
  902. ! are not at a water point.
  903. IF (landmask(i,j) .GT. 0.5 ) THEN
  904. IF ( sf_surface_physics .EQ. SLABSCHEME ) THEN
  905. st_near_sfc = 0 ! Check if there is a soil layer above 40 cm
  906. DO k = 1, num_st_levels_input
  907. IF ( st_levels_input(k) .LE. 40 ) THEN
  908. st_near_sfc = 1
  909. END IF
  910. END DO
  911. IF ( ( flag_tavgsfc == 1 ) .OR. ( st_near_sfc == 1 ) ) THEN
  912. tmn(i,j) = tmn(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  913. ELSE
  914. tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
  915. END IF
  916. END IF
  917. tsk(i,j) = tsk(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  918. IF ( flag_soil_layers == 1 ) THEN
  919. DO k = 2, num_st_levels_input+1
  920. st_input(i,k,j) = st_input(i,k,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  921. END DO
  922. ELSE
  923. DO k = 1, num_st_levels_input
  924. st_input(i,k,j) = st_input(i,k,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  925. END DO
  926. END IF
  927. END IF
  928. END DO
  929. END DO
  930. END IF
  931. END SUBROUTINE adjust_soil_temp_new
  932. SUBROUTINE init_soil_depth_1 ( zs , dzs , num_soil_layers )
  933. IMPLICIT NONE
  934. INTEGER, INTENT(IN) :: num_soil_layers
  935. REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
  936. INTEGER :: l
  937. ! Define layers (top layer = 0.01 m). Double the thicknesses at each step (dzs values).
  938. ! The distance from the ground level to the midpoint of the layer is given by zs.
  939. ! ------- Ground Level ---------- || || || ||
  940. ! || || || || zs(1) = 0.005 m
  941. ! -- -- -- -- -- -- -- -- -- || || || \/
  942. ! || || ||
  943. ! ----------------------------------- || || || \/ dzs(1) = 0.01 m
  944. ! || || ||
  945. ! || || || zs(2) = 0.02
  946. ! -- -- -- -- -- -- -- -- -- || || \/
  947. ! || ||
  948. ! || ||
  949. ! ----------------------------------- || || \/ dzs(2) = 0.02 m
  950. ! || ||
  951. ! || ||
  952. ! || ||
  953. ! || || zs(3) = 0.05
  954. ! -- -- -- -- -- -- -- -- -- || \/
  955. ! ||
  956. ! ||
  957. ! ||
  958. ! ||
  959. ! ----------------------------------- \/ dzs(3) = 0.04 m
  960. IF ( num_soil_layers .NE. 5 ) THEN
  961. PRINT '(A)','Usually, the 5-layer diffusion uses 5 layers. Change this in the namelist.'
  962. CALL wrf_error_fatal ( '5-layer_diffusion_uses_5_layers' )
  963. END IF
  964. dzs(1)=.01
  965. zs(1)=.5*dzs(1)
  966. DO l=2,num_soil_layers
  967. dzs(l)=2*dzs(l-1)
  968. zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
  969. ENDDO
  970. END SUBROUTINE init_soil_depth_1
  971. SUBROUTINE init_soil_depth_2 ( zs , dzs , num_soil_layers )
  972. IMPLICIT NONE
  973. INTEGER, INTENT(IN) :: num_soil_layers
  974. REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
  975. INTEGER :: l
  976. dzs = (/ 0.1 , 0.3 , 0.6 , 1.0 /)
  977. IF ( num_soil_layers .NE. 4 ) THEN
  978. PRINT '(A)','Usually, the LSM uses 4 layers. Change this in the namelist.'
  979. CALL wrf_error_fatal ( 'LSM_uses_4_layers' )
  980. END IF
  981. zs(1)=.5*dzs(1)
  982. DO l=2,num_soil_layers
  983. zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
  984. ENDDO
  985. END SUBROUTINE init_soil_depth_2
  986. SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_layers )
  987. IMPLICIT NONE
  988. INTEGER, INTENT(IN) :: num_soil_layers
  989. REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
  990. INTEGER :: l
  991. CHARACTER (LEN=132) :: message
  992. ! in RUC LSM ZS - soil levels, and DZS - soil layer thicknesses, not used
  993. ! ZS is specified in the namelist: num_soil_layers = 6 or 9.
  994. ! Other options with number of levels are possible, but
  995. ! WRF users should change consistently the namelist entry with the
  996. ! ZS array in this subroutine.
  997. IF ( num_soil_layers .EQ. 6) THEN
  998. zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 1.60 , 3.00 /)
  999. ! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
  1000. ELSEIF ( num_soil_layers .EQ. 9) THEN
  1001. zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 0.60, 1.00, 1.60 , 2.20, 3.00 /)
  1002. ! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
  1003. ENDIF
  1004. IF ( num_soil_layers .EQ. 4 .OR. num_soil_layers .EQ. 5 ) THEN
  1005. write (message, FMT='(A)') 'The RUC LSM uses 6, 9 or more levels. Change this in the namelist.'
  1006. CALL wrf_error_fatal ( message )
  1007. END IF
  1008. END SUBROUTINE init_soil_depth_3
  1009. SUBROUTINE init_soil_depth_7 ( zs , dzs , num_soil_layers )
  1010. IMPLICIT NONE
  1011. INTEGER, INTENT(IN) :: num_soil_layers
  1012. REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
  1013. INTEGER :: l
  1014. dzs = (/ 0.01 , 0.99 /)
  1015. IF ( num_soil_layers .NE. 2 ) THEN
  1016. PRINT '(A)','Usually, the PX LSM uses 2 layers. Change this in the namelist.'
  1017. CALL wrf_error_fatal ( 'PXLSM_uses_2_layers' )
  1018. END IF
  1019. zs(1) = 0.5 * dzs(1)
  1020. zs(2) = dzs(1) + 0.5 * dzs(2)
  1021. END SUBROUTINE init_soil_depth_7
  1022. SUBROUTINE init_soil_depth_8 ( zs , dzs , num_soil_layers )
  1023. IMPLICIT NONE
  1024. INTEGER, INTENT(IN) :: num_soil_layers
  1025. REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
  1026. INTEGER :: l
  1027. zs(1) = 0.00
  1028. zs(2) = 0.05
  1029. zs(3) = 1.05
  1030. END SUBROUTINE init_soil_depth_8
  1031. SUBROUTINE init_soil_1_real ( tsk , tmn , tslb , zs , dzs , &
  1032. num_soil_layers , real_data_init_type , &
  1033. landmask , sst , flag_sst , &
  1034. ids , ide , jds , jde , kds , kde , &
  1035. ims , ime , jms , jme , kms , kme , &
  1036. its , ite , jts , jte , kts , kte )
  1037. IMPLICIT NONE
  1038. INTEGER , INTENT(IN) :: num_soil_layers , real_data_init_type , &
  1039. ids , ide , jds , jde , kds , kde , &
  1040. ims , ime , jms , jme , kms , kme , &
  1041. its , ite , jts , jte , kts , kte
  1042. INTEGER , INTENT(IN) :: flag_sst
  1043. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
  1044. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tsk , tmn
  1045. REAL , DIMENSION(num_soil_layers) :: zs , dzs
  1046. REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb
  1047. INTEGER :: i , j , l
  1048. ! Soil temperature is linearly interpolated between the skin temperature (taken to be at a
  1049. ! depth of 0.5 cm) and the deep soil, annual temperature (taken to be at a depth of 23 cm).
  1050. ! The tslb(i,1,j) is the skin temperature, and the tslb(i,num_soil_layers,j) level is the
  1051. ! annual mean temperature.
  1052. DO j = jts , MIN(jde-1,jte)
  1053. DO i = its , MIN(ide-1,ite)
  1054. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1055. IF ( landmask(i,j) .GT. 0.5 ) THEN
  1056. DO l = 1 , num_soil_layers
  1057. tslb(i,l,j)= ( tsk(i,j) * ( zs(num_soil_layers) - zs(l) ) + &
  1058. tmn(i,j) * ( zs( l) - zs(1) ) ) / &
  1059. ( zs(num_soil_layers) - zs(1) )
  1060. END DO
  1061. ELSE
  1062. IF ( ( real_data_init_type .EQ. 1 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
  1063. DO l = 1 , num_soil_layers
  1064. tslb(i,l,j)= sst(i,j)
  1065. END DO
  1066. ELSE
  1067. DO l = 1 , num_soil_layers
  1068. tslb(i,l,j)= tsk(i,j)
  1069. END DO
  1070. END IF
  1071. END IF
  1072. END DO
  1073. END DO
  1074. END SUBROUTINE init_soil_1_real
  1075. SUBROUTINE init_soil_1_ideal(tsk,tmn,tslb,xland, &
  1076. ivgtyp,ZS,DZS,num_soil_layers, &
  1077. ids,ide, jds,jde, kds,kde, &
  1078. ims,ime, jms,jme, kms,kme, &
  1079. its,ite, jts,jte, kts,kte )
  1080. IMPLICIT NONE
  1081. INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
  1082. ims,ime, jms,jme, kms,kme, &
  1083. its,ite, jts,jte, kts,kte
  1084. INTEGER, INTENT(IN ) :: num_soil_layers
  1085. REAL, DIMENSION( ims:ime , num_soil_layers , jms:jme ), INTENT(OUT) :: tslb
  1086. REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: xland
  1087. INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: ivgtyp
  1088. REAL, DIMENSION(1:), INTENT(IN) :: dzs,zs
  1089. REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: tsk, tmn
  1090. ! Lcal variables.
  1091. INTEGER :: l,j,i,itf,jtf
  1092. itf=MIN(ite,ide-1)
  1093. jtf=MIN(jte,jde-1)
  1094. IF (num_soil_layers.NE.1)THEN
  1095. DO j=jts,jtf
  1096. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1097. DO l=1,num_soil_layers
  1098. DO i=its,itf
  1099. tslb(i,l,j)=( tsk(i,j)*(zs(num_soil_layers)-zs(l)) + tmn(i,j)*(zs(l)-zs(1)) ) / &
  1100. ( zs(num_soil_layers)-zs(1) )
  1101. ENDDO
  1102. ENDDO
  1103. ENDDO
  1104. ENDIF
  1105. ! DO j=jts,jtf
  1106. ! DO i=its,itf
  1107. ! xland(i,j) = 2
  1108. ! ivgtyp(i,j) = 7
  1109. ! ENDDO
  1110. ! ENDDO
  1111. END SUBROUTINE init_soil_1_ideal
  1112. SUBROUTINE init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
  1113. st_input , sm_input , sw_input , landmask , sst , &
  1114. zs , dzs , &
  1115. st_levels_input , sm_levels_input , sw_levels_input , &
  1116. num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  1117. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
  1118. flag_sst , flag_soil_layers , flag_soil_levels , &
  1119. ids , ide , jds , jde , kds , kde , &
  1120. ims , ime , jms , jme , kms , kme , &
  1121. its , ite , jts , jte , kts , kte )
  1122. IMPLICIT NONE
  1123. INTEGER , INTENT(IN) :: num_soil_layers , &
  1124. num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  1125. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
  1126. ids , ide , jds , jde , kds , kde , &
  1127. ims , ime , jms , jme , kms , kme , &
  1128. its , ite , jts , jte , kts , kte
  1129. INTEGER , INTENT(IN) :: flag_sst, flag_soil_layers, flag_soil_levels
  1130. INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
  1131. INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
  1132. INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
  1133. REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
  1134. REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
  1135. REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
  1136. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
  1137. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
  1138. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
  1139. REAL , DIMENSION(num_soil_layers) :: zs , dzs
  1140. REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
  1141. REAL , ALLOCATABLE , DIMENSION(:) :: zhave
  1142. INTEGER :: i , j , l , lout , lin , lwant , lhave , num
  1143. REAL :: temp
  1144. LOGICAL :: found_levels
  1145. CHARACTER (LEN=132) :: message
  1146. ! Are there any soil temp and moisture levels - ya know, they are mandatory.
  1147. num = num_st_levels_input * num_sm_levels_input
  1148. IF ( num .GE. 1 ) THEN
  1149. ! Ordered levels that we have data for.
  1150. !tgs add option to initialize from RUCLSM
  1151. IF ( flag_soil_levels == 1 ) THEN
  1152. write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
  1153. CALL wrf_message ( message )
  1154. ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) ) )
  1155. ELSE
  1156. write(message, FMT='(A)') ' Assume Noah LSM input'
  1157. CALL wrf_message ( message )
  1158. ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
  1159. END IF
  1160. ! Sort the levels for temperature.
  1161. outert : DO lout = 1 , num_st_levels_input-1
  1162. innert : DO lin = lout+1 , num_st_levels_input
  1163. IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
  1164. temp = st_levels_input(lout)
  1165. st_levels_input(lout) = st_levels_input(lin)
  1166. st_levels_input(lin) = NINT(temp)
  1167. DO j = jts , MIN(jde-1,jte)
  1168. DO i = its , MIN(ide-1,ite)
  1169. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1170. temp = st_input(i,lout+1,j)
  1171. st_input(i,lout+1,j) = st_input(i,lin+1,j)
  1172. st_input(i,lin+1,j) = temp
  1173. END DO
  1174. END DO
  1175. END IF
  1176. END DO innert
  1177. END DO outert
  1178. !tgs add IF
  1179. IF ( flag_soil_layers == 1 ) THEN
  1180. DO j = jts , MIN(jde-1,jte)
  1181. DO i = its , MIN(ide-1,ite)
  1182. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1183. st_input(i,1,j) = tsk(i,j)
  1184. st_input(i,num_st_levels_input+2,j) = tmn(i,j)
  1185. END DO
  1186. END DO
  1187. ENDIF
  1188. ! Sort the levels for moisture.
  1189. outerm: DO lout = 1 , num_sm_levels_input-1
  1190. innerm : DO lin = lout+1 , num_sm_levels_input
  1191. IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
  1192. temp = sm_levels_input(lout)
  1193. sm_levels_input(lout) = sm_levels_input(lin)
  1194. sm_levels_input(lin) = NINT(temp)
  1195. DO j = jts , MIN(jde-1,jte)
  1196. DO i = its , MIN(ide-1,ite)
  1197. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1198. temp = sm_input(i,lout+1,j)
  1199. sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
  1200. sm_input(i,lin+1,j) = temp
  1201. END DO
  1202. END DO
  1203. END IF
  1204. END DO innerm
  1205. END DO outerm
  1206. !tgs add IF
  1207. IF ( flag_soil_layers == 1 ) THEN
  1208. DO j = jts , MIN(jde-1,jte)
  1209. DO i = its , MIN(ide-1,ite)
  1210. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1211. sm_input(i,1,j) = sm_input(i,2,j)
  1212. sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
  1213. END DO
  1214. END DO
  1215. ENDIF
  1216. ! Sort the levels for liquid moisture.
  1217. outerw: DO lout = 1 , num_sw_levels_input-1
  1218. innerw : DO lin = lout+1 , num_sw_levels_input
  1219. IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
  1220. temp = sw_levels_input(lout)
  1221. sw_levels_input(lout) = sw_levels_input(lin)
  1222. sw_levels_input(lin) = NINT(temp)
  1223. DO j = jts , MIN(jde-1,jte)
  1224. DO i = its , MIN(ide-1,ite)
  1225. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1226. temp = sw_input(i,lout+1,j)
  1227. sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
  1228. sw_input(i,lin+1,j) = temp
  1229. END DO
  1230. END DO
  1231. END IF
  1232. END DO innerw
  1233. END DO outerw
  1234. IF ( num_sw_levels_input .GT. 1 ) THEN
  1235. DO j = jts , MIN(jde-1,jte)
  1236. DO i = its , MIN(ide-1,ite)
  1237. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1238. sw_input(i,1,j) = sw_input(i,2,j)
  1239. sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
  1240. END DO
  1241. END DO
  1242. END IF
  1243. found_levels = .TRUE.
  1244. ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
  1245. found_levels = .FALSE.
  1246. ELSE
  1247. CALL wrf_error_fatal ( &
  1248. 'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
  1249. END IF
  1250. ! Is it OK to continue?
  1251. IF ( found_levels ) THEN
  1252. !tgs: Here are the levels that we have from the input for temperature.
  1253. IF ( flag_soil_levels == 1 ) THEN
  1254. DO l = 1 , num_st_levels_input
  1255. zhave(l) = st_levels_input(l) / 100.
  1256. END DO
  1257. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  1258. z_wantt : DO lwant = 1 , num_soil_layers
  1259. z_havet : DO lhave = 1 , num_st_levels_input -1
  1260. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  1261. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  1262. DO j = jts , MIN(jde-1,jte)
  1263. DO i = its , MIN(ide-1,ite)
  1264. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1265. tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
  1266. st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  1267. ( zhave(lhave+1) - zhave(lhave) )
  1268. END DO
  1269. END DO
  1270. EXIT z_havet
  1271. END IF
  1272. END DO z_havet
  1273. END DO z_wantt
  1274. ELSE
  1275. ! Here are the levels that we have from the input for temperature. The input levels plus
  1276. ! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
  1277. zhave(1) = 0.
  1278. DO l = 1 , num_st_levels_input
  1279. zhave(l+1) = st_levels_input(l) / 100.
  1280. END DO
  1281. zhave(num_st_levels_input+2) = 300. / 100.
  1282. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  1283. z_wantt_2: DO lwant = 1 , num_soil_layers
  1284. z_havet_2 : DO lhave = 1 , num_st_levels_input +2 -1
  1285. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  1286. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  1287. DO j = jts , MIN(jde-1,jte)
  1288. DO i = its , MIN(ide-1,ite)
  1289. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1290. tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
  1291. st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  1292. ( zhave(lhave+1) - zhave(lhave) )
  1293. END DO
  1294. END DO
  1295. EXIT z_havet_2
  1296. END IF
  1297. END DO z_havet_2
  1298. END DO z_wantt_2
  1299. END IF
  1300. !tgs:
  1301. IF ( flag_soil_levels == 1 ) THEN
  1302. DO l = 1 , num_sm_levels_input
  1303. zhave(l) = sm_levels_input(l) / 100.
  1304. END DO
  1305. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  1306. z_wantm : DO lwant = 1 , num_soil_layers
  1307. z_havem : DO lhave = 1 , num_sm_levels_input -1
  1308. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  1309. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  1310. DO j = jts , MIN(jde-1,jte)
  1311. DO i = its , MIN(ide-1,ite)
  1312. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1313. smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
  1314. sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  1315. ( zhave(lhave+1) - zhave(lhave) )
  1316. END DO
  1317. END DO
  1318. EXIT z_havem
  1319. END IF
  1320. END DO z_havem
  1321. END DO z_wantm
  1322. ELSE
  1323. ! Here are the levels that we have from the input for moisture. The input levels plus
  1324. ! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
  1325. ! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
  1326. ! as the most deep layer's value.
  1327. zhave(1) = 0.
  1328. DO l = 1 , num_sm_levels_input
  1329. zhave(l+1) = sm_levels_input(l) / 100.
  1330. END DO
  1331. zhave(num_sm_levels_input+2) = 300. / 100.
  1332. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  1333. z_wantm_2 : DO lwant = 1 , num_soil_layers
  1334. z_havem_2 : DO lhave = 1 , num_sm_levels_input +2 -1
  1335. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  1336. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  1337. DO j = jts , MIN(jde-1,jte)
  1338. DO i = its , MIN(ide-1,ite)
  1339. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1340. smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
  1341. sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  1342. ( zhave(lhave+1) - zhave(lhave) )
  1343. END DO
  1344. END DO
  1345. EXIT z_havem_2
  1346. END IF
  1347. END DO z_havem_2
  1348. END DO z_wantm_2
  1349. ENDIF
  1350. ! Any liquid soil moisture to worry about?
  1351. IF ( num_sw_levels_input .GT. 1 ) THEN
  1352. zhave(1) = 0.
  1353. DO l = 1 , num_sw_levels_input
  1354. zhave(l+1) = sw_levels_input(l) / 100.
  1355. END DO
  1356. zhave(num_sw_levels_input+2) = 300. / 100.
  1357. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  1358. z_wantw : DO lwant = 1 , num_soil_layers
  1359. z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
  1360. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  1361. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  1362. DO j = jts , MIN(jde-1,jte)
  1363. DO i = its , MIN(ide-1,ite)
  1364. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1365. sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
  1366. sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  1367. ( zhave(lhave+1) - zhave(lhave) )
  1368. END DO
  1369. END DO
  1370. EXIT z_havew
  1371. END IF
  1372. END DO z_havew
  1373. END DO z_wantw
  1374. END IF
  1375. ! Over water, put in reasonable values for soil temperature and moisture. These won't be
  1376. ! used, but they will make a more continuous plot.
  1377. IF ( flag_sst .EQ. 1 ) THEN
  1378. DO j = jts , MIN(jde-1,jte)
  1379. DO i = its , MIN(ide-1,ite)
  1380. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1381. IF ( landmask(i,j) .LT. 0.5 ) THEN
  1382. DO l = 1 , num_soil_layers
  1383. tslb(i,l,j)= sst(i,j)
  1384. smois(i,l,j)= 1.0
  1385. sh2o (i,l,j)= 1.0
  1386. END DO
  1387. END IF
  1388. END DO
  1389. END DO
  1390. ELSE
  1391. DO j = jts , MIN(jde-1,jte)
  1392. DO i = its , MIN(ide-1,ite)
  1393. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1394. IF ( landmask(i,j) .LT. 0.5 ) THEN
  1395. DO l = 1 , num_soil_layers
  1396. tslb(i,l,j)= tsk(i,j)
  1397. smois(i,l,j)= 1.0
  1398. sh2o (i,l,j)= 1.0
  1399. END DO
  1400. END IF
  1401. END DO
  1402. END DO
  1403. END IF
  1404. DEALLOCATE (zhave)
  1405. END IF
  1406. END SUBROUTINE init_soil_2_real
  1407. SUBROUTINE init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
  1408. ivgtyp,isltyp,tslb,smois,tmn, &
  1409. num_soil_layers, &
  1410. ids,ide, jds,jde, kds,kde, &
  1411. ims,ime, jms,jme, kms,kme, &
  1412. its,ite, jts,jte, kts,kte )
  1413. IMPLICIT NONE
  1414. INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
  1415. ims,ime, jms,jme, kms,kme, &
  1416. its,ite, jts,jte, kts,kte
  1417. INTEGER, INTENT(IN) ::num_soil_layers
  1418. REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(OUT) :: smois, tslb
  1419. REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: xland
  1420. REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: snow, canwat, xice, vegfra, tmn
  1421. INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
  1422. INTEGER :: icm,jcm,itf,jtf
  1423. INTEGER :: i,j,l
  1424. itf=min0(ite,ide-1)
  1425. jtf=min0(jte,jde-1)
  1426. icm = ide/2
  1427. jcm = jde/2
  1428. DO j=jts,jtf
  1429. DO l=1,num_soil_layers
  1430. DO i=its,itf
  1431. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1432. if (xland(i,j) .lt. 1.5) then
  1433. smois(i,1,j)=0.30
  1434. smois(i,2,j)=0.30
  1435. smois(i,3,j)=0.30
  1436. smois(i,4,j)=0.30
  1437. tslb(i,1,j)=290.
  1438. tslb(i,2,j)=290.
  1439. tslb(i,3,j)=290.
  1440. tslb(i,4,j)=290.
  1441. endif
  1442. ENDDO
  1443. ENDDO
  1444. ENDDO
  1445. END SUBROUTINE init_soil_2_ideal
  1446. SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , &
  1447. st_input , sm_input , landmask, sst, &
  1448. zs , dzs , &
  1449. st_levels_input , sm_levels_input , &
  1450. num_soil_layers , num_st_levels_input , num_sm_levels_input , &
  1451. num_st_levels_alloc , num_sm_levels_alloc , &
  1452. flag_sst , flag_soil_layers , flag_soil_levels , &
  1453. ids , ide , jds , jde , kds , kde , &
  1454. ims , ime , jms , jme , kms , kme , &
  1455. its , ite , jts , jte , kts , kte )
  1456. IMPLICIT NONE
  1457. INTEGER , INTENT(IN) :: num_soil_layers , &
  1458. num_st_levels_input , num_sm_levels_input , &
  1459. num_st_levels_alloc , num_sm_levels_alloc , &
  1460. ids , ide , jds , jde , kds , kde , &
  1461. ims , ime , jms , jme , kms , kme , &
  1462. its , ite , jts , jte , kts , kte
  1463. INTEGER , INTENT(IN) :: flag_sst, flag_soil_layers, flag_soil_levels
  1464. INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
  1465. INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
  1466. REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
  1467. REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
  1468. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
  1469. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
  1470. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
  1471. REAL , DIMENSION(num_soil_layers) :: zs , dzs
  1472. REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois
  1473. REAL , ALLOCATABLE , DIMENSION(:) :: zhave
  1474. INTEGER :: i , j , l , lout , lin , lwant , lhave, k
  1475. REAL :: temp
  1476. CHARACTER (LEN=132) :: message
  1477. ! Allocate the soil layer array used for interpolating.
  1478. IF ( ( num_st_levels_input .LE. 0 ) .OR. &
  1479. ( num_sm_levels_input .LE. 0 ) ) THEN
  1480. write (message, FMT='(A)')&
  1481. 'No input soil level data (either temperature or moisture, or both are missing). Required for RUC LSM.'
  1482. CALL wrf_error_fatal ( message )
  1483. ELSE
  1484. IF ( flag_soil_levels == 1 ) THEN
  1485. write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
  1486. CALL wrf_message ( message )
  1487. ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input) ) )
  1488. ELSE
  1489. write(message, FMT='(A)') ' Assume non-RUC LSM input'
  1490. CALL wrf_message ( message )
  1491. ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers) ) )
  1492. END IF
  1493. END IF
  1494. ! Sort the levels for temperature.
  1495. outert : DO lout = 1 , num_st_levels_input-1
  1496. innert : DO lin = lout+1 , num_st_levels_input
  1497. IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
  1498. temp = st_levels_input(lout)
  1499. st_levels_input(lout) = st_levels_input(lin)
  1500. st_levels_input(lin) = NINT(temp)
  1501. DO j = jts , MIN(jde-1,jte)
  1502. DO i = its , MIN(ide-1,ite)
  1503. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1504. temp = st_input(i,lout,j)
  1505. st_input(i,lout,j) = st_input(i,lin,j)
  1506. st_input(i,lin,j) = temp
  1507. END DO
  1508. END DO
  1509. END IF
  1510. END DO innert
  1511. END DO outert
  1512. IF ( flag_soil_layers == 1 ) THEN
  1513. DO j = jts , MIN(jde-1,jte)
  1514. DO i = its , MIN(ide-1,ite)
  1515. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1516. st_input(i,1,j) = tsk(i,j)
  1517. st_input(i,num_st_levels_input+2,j) = tmn(i,j)
  1518. END DO
  1519. END DO
  1520. END IF
  1521. ! Sort the levels for moisture.
  1522. outerm: DO lout = 1 , num_sm_levels_input-1
  1523. innerm : DO lin = lout+1 , num_sm_levels_input
  1524. IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
  1525. temp = sm_levels_input(lout)
  1526. sm_levels_input(lout) = sm_levels_input(lin)
  1527. sm_levels_input(lin) = NINT(temp)
  1528. DO j = jts , MIN(jde-1,jte)
  1529. DO i = its , MIN(ide-1,ite)
  1530. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1531. temp = sm_input(i,lout,j)
  1532. sm_input(i,lout,j) = sm_input(i,lin,j)
  1533. sm_input(i,lin,j) = temp
  1534. END DO
  1535. END DO
  1536. END IF
  1537. END DO innerm
  1538. END DO outerm
  1539. IF ( flag_soil_layers == 1 ) THEN
  1540. DO j = jts , MIN(jde-1,jte)
  1541. DO i = its , MIN(ide-1,ite)
  1542. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1543. sm_input(i,1,j) = (sm_input(i,2,j)-sm_input(i,3,j))/ &
  1544. (st_levels_input(2)-st_levels_input(1))*st_levels_input(1)+ &
  1545. sm_input(i,2,j)
  1546. ! sm_input(i,1,j) = sm_input(i,2,j)
  1547. sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
  1548. END DO
  1549. END DO
  1550. END IF
  1551. ! Here are the levels that we have from the input for temperature.
  1552. IF ( flag_soil_levels == 1 ) THEN
  1553. DO l = 1 , num_st_levels_input
  1554. zhave(l) = st_levels_input(l) / 100.
  1555. END DO
  1556. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  1557. z_wantt : DO lwant = 1 , num_soil_layers
  1558. z_havet : DO lhave = 1 , num_st_levels_input -1
  1559. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  1560. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  1561. DO j = jts , MIN(jde-1,jte)
  1562. DO i = its , MIN(ide-1,ite)
  1563. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1564. tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
  1565. st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  1566. ( zhave(lhave+1) - zhave(lhave) )
  1567. END DO
  1568. END DO
  1569. EXIT z_havet
  1570. END IF
  1571. END DO z_havet
  1572. END DO z_wantt
  1573. ELSE
  1574. zhave(1) = 0.
  1575. DO l = 1 , num_st_levels_input
  1576. zhave(l+1) = st_levels_input(l) / 100.
  1577. END DO
  1578. zhave(num_st_levels_input+2) = 300. / 100.
  1579. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  1580. z_wantt_2 : DO lwant = 1 , num_soil_layers
  1581. z_havet_2 : DO lhave = 1 , num_st_levels_input +2
  1582. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  1583. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  1584. DO j = jts , MIN(jde-1,jte)
  1585. DO i = its , MIN(ide-1,ite)
  1586. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1587. tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
  1588. st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  1589. ( zhave(lhave+1) - zhave(lhave) )
  1590. END DO
  1591. END DO
  1592. EXIT z_havet_2
  1593. END IF
  1594. END DO z_havet_2
  1595. END DO z_wantt_2
  1596. END IF
  1597. ! Here are the levels that we have from the input for moisture.
  1598. IF ( flag_soil_levels .EQ. 1 ) THEN
  1599. DO l = 1 , num_sm_levels_input
  1600. zhave(l) = sm_levels_input(l) / 100.
  1601. END DO
  1602. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  1603. z_wantm : DO lwant = 1 , num_soil_layers
  1604. z_havem : DO lhave = 1 , num_sm_levels_input -1
  1605. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  1606. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  1607. DO j = jts , MIN(jde-1,jte)
  1608. DO i = its , MIN(ide-1,ite)
  1609. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1610. smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
  1611. sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  1612. ( zhave(lhave+1) - zhave(lhave) )
  1613. END DO
  1614. END DO
  1615. EXIT z_havem
  1616. END IF
  1617. END DO z_havem
  1618. END DO z_wantm
  1619. ELSE
  1620. zhave(1) = 0.
  1621. DO l = 1 , num_sm_levels_input
  1622. zhave(l+1) = sm_levels_input(l) / 100.
  1623. END DO
  1624. zhave(num_sm_levels_input+2) = 300. / 100.
  1625. z_wantm_2 : DO lwant = 1 , num_soil_layers
  1626. z_havem_2 : DO lhave = 1 , num_sm_levels_input +2
  1627. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  1628. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  1629. DO j = jts , MIN(jde-1,jte)
  1630. DO i = its , MIN(ide-1,ite)
  1631. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1632. smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
  1633. sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  1634. ( zhave(lhave+1) - zhave(lhave) )
  1635. END DO
  1636. END DO
  1637. EXIT z_havem_2
  1638. END IF
  1639. END DO z_havem_2
  1640. END DO z_wantm_2
  1641. END IF
  1642. ! Over water, put in reasonable values for soil temperature and moisture. These won't be
  1643. ! used, but they will make a more continuous plot.
  1644. IF ( flag_sst .EQ. 1 ) THEN
  1645. DO j = jts , MIN(jde-1,jte)
  1646. DO i = its , MIN(ide-1,ite)
  1647. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1648. IF ( landmask(i,j) .LT. 0.5 ) THEN
  1649. DO l = 1 , num_soil_layers
  1650. tslb(i,l,j) = sst(i,j)
  1651. tsk(i,j) = sst(i,j)
  1652. smois(i,l,j)= 1.0
  1653. END DO
  1654. END IF
  1655. END DO
  1656. END DO
  1657. ELSE
  1658. DO j = jts , MIN(jde-1,jte)
  1659. DO i = its , MIN(ide-1,ite)
  1660. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1661. IF ( landmask(i,j) .LT. 0.5 ) THEN
  1662. DO l = 1 , num_soil_layers
  1663. tslb(i,l,j)= tsk(i,j)
  1664. smois(i,l,j)= 1.0
  1665. END DO
  1666. END IF
  1667. END DO
  1668. END DO
  1669. END IF
  1670. DEALLOCATE (zhave)
  1671. END SUBROUTINE init_soil_3_real
  1672. SUBROUTINE init_soil_7_real ( tsk , tmn , smois , sh2o , tslb , &
  1673. st_input , sm_input , sw_input , landmask , sst , &
  1674. zs , dzs , &
  1675. st_levels_input , sm_levels_input , sw_levels_input , &
  1676. num_soil_layers , num_st_levels_input , &
  1677. num_sm_levels_input , num_sw_levels_input , &
  1678. num_st_levels_alloc , num_sm_levels_alloc , &
  1679. num_sw_levels_alloc , &
  1680. flag_sst , flag_soil_layers , flag_soil_levels , &
  1681. ids , ide , jds , jde , kds , kde , &
  1682. ims , ime , jms , jme , kms , kme , &
  1683. its , ite , jts , jte , kts , kte )
  1684. ! for soil temperature and moisture initialization for the PX LSM
  1685. IMPLICIT NONE
  1686. INTEGER , INTENT(IN) :: num_soil_layers , &
  1687. num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  1688. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
  1689. ids , ide , jds , jde , kds , kde , &
  1690. ims , ime , jms , jme , kms , kme , &
  1691. its , ite , jts , jte , kts , kte
  1692. INTEGER , INTENT(IN) :: flag_sst, flag_soil_layers, flag_soil_levels
  1693. INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
  1694. INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
  1695. INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
  1696. REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
  1697. REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
  1698. REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
  1699. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
  1700. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
  1701. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
  1702. REAL , DIMENSION(num_soil_layers) :: zs , dzs
  1703. REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
  1704. REAL , ALLOCATABLE , DIMENSION(:) :: zhave
  1705. INTEGER :: i , j , l , lout , lin , lwant , lhave , num
  1706. REAL :: temp
  1707. LOGICAL :: found_levels
  1708. ! Are there any soil temp and moisture levels - ya know, they are mandatory.
  1709. num = num_st_levels_input * num_sm_levels_input
  1710. IF ( num .GE. 1 ) THEN
  1711. ! Ordered levels that we have data for.
  1712. ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
  1713. ! Sort the levels for temperature.
  1714. outert : DO lout = 1 , num_st_levels_input-1
  1715. innert : DO lin = lout+1 , num_st_levels_input
  1716. IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
  1717. temp = st_levels_input(lout)
  1718. st_levels_input(lout) = st_levels_input(lin)
  1719. st_levels_input(lin) = NINT(temp)
  1720. DO j = jts , MIN(jde-1,jte)
  1721. DO i = its , MIN(ide-1,ite)
  1722. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1723. temp = st_input(i,lout+1,j)
  1724. st_input(i,lout+1,j) = st_input(i,lin+1,j)
  1725. st_input(i,lin+1,j) = temp
  1726. END DO
  1727. END DO
  1728. END IF
  1729. END DO innert
  1730. END DO outert
  1731. DO j = jts , MIN(jde-1,jte)
  1732. DO i = its , MIN(ide-1,ite)
  1733. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1734. st_input(i,1,j) = tsk(i,j)
  1735. st_input(i,num_st_levels_input+2,j) = tmn(i,j)
  1736. END DO
  1737. END DO
  1738. ! Sort the levels for moisture.
  1739. outerm: DO lout = 1 , num_sm_levels_input-1
  1740. innerm : DO lin = lout+1 , num_sm_levels_input
  1741. IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
  1742. temp = sm_levels_input(lout)
  1743. sm_levels_input(lout) = sm_levels_input(lin)
  1744. sm_levels_input(lin) = NINT(temp)
  1745. DO j = jts , MIN(jde-1,jte)
  1746. DO i = its , MIN(ide-1,ite)
  1747. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1748. temp = sm_input(i,lout+1,j)
  1749. sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
  1750. sm_input(i,lin+1,j) = temp
  1751. END DO
  1752. END DO
  1753. END IF
  1754. END DO innerm
  1755. END DO outerm
  1756. DO j = jts , MIN(jde-1,jte)
  1757. DO i = its , MIN(ide-1,ite)
  1758. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1759. sm_input(i,1,j) = sm_input(i,2,j)
  1760. sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
  1761. END DO
  1762. END DO
  1763. ! Sort the levels for liquid moisture.
  1764. outerw: DO lout = 1 , num_sw_levels_input-1
  1765. innerw : DO lin = lout+1 , num_sw_levels_input
  1766. IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
  1767. temp = sw_levels_input(lout)
  1768. sw_levels_input(lout) = sw_levels_input(lin)
  1769. sw_levels_input(lin) = NINT(temp)
  1770. DO j = jts , MIN(jde-1,jte)
  1771. DO i = its , MIN(ide-1,ite)
  1772. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1773. temp = sw_input(i,lout+1,j)
  1774. sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
  1775. sw_input(i,lin+1,j) = temp
  1776. END DO
  1777. END DO
  1778. END IF
  1779. END DO innerw
  1780. END DO outerw
  1781. IF ( num_sw_levels_input .GT. 1 ) THEN
  1782. DO j = jts , MIN(jde-1,jte)
  1783. DO i = its , MIN(ide-1,ite)
  1784. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1785. sw_input(i,1,j) = sw_input(i,2,j)
  1786. sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
  1787. END DO
  1788. END DO
  1789. END IF
  1790. found_levels = .TRUE.
  1791. ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
  1792. found_levels = .FALSE.
  1793. ELSE
  1794. CALL wrf_error_fatal ( &
  1795. 'No input soil level data (temperature, moisture or liquid, or all are missing). Required for PX LSM.' )
  1796. END IF
  1797. ! Is it OK to continue?
  1798. IF ( found_levels ) THEN
  1799. ! Here are the levels that we have from the input for temperature. The input levels plus
  1800. ! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
  1801. zhave(1) = 0.
  1802. DO l = 1 , num_st_levels_input
  1803. zhave(l+1) = st_levels_input(l) / 100.
  1804. END DO
  1805. zhave(num_st_levels_input+2) = 300. / 100.
  1806. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  1807. z_wantt : DO lwant = 1 , num_soil_layers
  1808. z_havet : DO lhave = 1 , num_st_levels_input +2 -1
  1809. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  1810. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  1811. DO j = jts , MIN(jde-1,jte)
  1812. DO i = its , MIN(ide-1,ite)
  1813. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1814. tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
  1815. st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  1816. ( zhave(lhave+1) - zhave(lhave) )
  1817. END DO
  1818. END DO
  1819. EXIT z_havet
  1820. END IF
  1821. END DO z_havet
  1822. END DO z_wantt
  1823. ! Here are the levels that we have from the input for moisture. The input levels plus
  1824. ! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
  1825. ! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
  1826. ! as the most deep layer's value.
  1827. zhave(1) = 0.
  1828. DO l = 1 , num_sm_levels_input
  1829. zhave(l+1) = sm_levels_input(l) / 100.
  1830. END DO
  1831. zhave(num_sm_levels_input+2) = 300. / 100.
  1832. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  1833. z_wantm : DO lwant = 1 , num_soil_layers
  1834. z_havem : DO lhave = 1 , num_sm_levels_input +2 -1
  1835. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  1836. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  1837. DO j = jts , MIN(jde-1,jte)
  1838. DO i = its , MIN(ide-1,ite)
  1839. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1840. smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
  1841. sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  1842. ( zhave(lhave+1) - zhave(lhave) )
  1843. END DO
  1844. END DO
  1845. EXIT z_havem
  1846. END IF
  1847. END DO z_havem
  1848. END DO z_wantm
  1849. ! Any liquid soil moisture to worry about?
  1850. IF ( num_sw_levels_input .GT. 1 ) THEN
  1851. zhave(1) = 0.
  1852. DO l = 1 , num_sw_levels_input
  1853. zhave(l+1) = sw_levels_input(l) / 100.
  1854. END DO
  1855. zhave(num_sw_levels_input+2) = 300. / 100.
  1856. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  1857. z_wantw : DO lwant = 1 , num_soil_layers
  1858. z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
  1859. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  1860. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  1861. DO j = jts , MIN(jde-1,jte)
  1862. DO i = its , MIN(ide-1,ite)
  1863. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1864. sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
  1865. sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  1866. ( zhave(lhave+1) - zhave(lhave) )
  1867. END DO
  1868. END DO
  1869. EXIT z_havew
  1870. END IF
  1871. END DO z_havew
  1872. END DO z_wantw
  1873. END IF
  1874. ! Over water, put in reasonable values for soil temperature and moisture. These won't be
  1875. ! used, but they will make a more continuous plot.
  1876. DO j = jts , MIN(jde-1,jte)
  1877. DO i = its , MIN(ide-1,ite)
  1878. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1879. tslb(i,1,j)= tsk(i,j)
  1880. tslb(i,2,j)= tmn(i,j)
  1881. END DO
  1882. END DO
  1883. IF ( flag_sst .EQ. 1 ) THEN
  1884. DO j = jts , MIN(jde-1,jte)
  1885. DO i = its , MIN(ide-1,ite)
  1886. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1887. IF ( landmask(i,j) .LT. 0.5 ) THEN
  1888. DO l = 1 , num_soil_layers
  1889. tslb(i,l,j)= sst(i,j)
  1890. smois(i,l,j)= 1.0
  1891. sh2o (i,l,j)= 1.0
  1892. END DO
  1893. END IF
  1894. END DO
  1895. END DO
  1896. ELSE
  1897. DO j = jts , MIN(jde-1,jte)
  1898. DO i = its , MIN(ide-1,ite)
  1899. IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
  1900. IF ( landmask(i,j) .LT. 0.5 ) THEN
  1901. DO l = 1 , num_soil_layers
  1902. tslb(i,l,j)= tsk(i,j)
  1903. smois(i,l,j)= 1.0
  1904. sh2o (i,l,j)= 1.0
  1905. END DO
  1906. END IF
  1907. END DO
  1908. END DO
  1909. END IF
  1910. DEALLOCATE (zhave)
  1911. END IF
  1912. END SUBROUTINE init_soil_7_real
  1913. !------------SSIB (fds 06/2010)-----------------------------
  1914. SUBROUTINE init_soil_8_real ( tsk , tmn , smois , sh2o , tslb , &
  1915. st_input , sm_input , sw_input , landmask , sst , &
  1916. zs , dzs , &
  1917. st_levels_input , sm_levels_input , sw_levels_input , &
  1918. num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  1919. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
  1920. flag_sst , flag_soil_layers , flag_soil_levels , &
  1921. st000010, st010200, sm000010, sm010200, &
  1922. st010040 , st040100 , st100200, &
  1923. sm010040 , sm040100 , sm100200, &
  1924. ids , ide , jds , jde , kds , kde , &
  1925. ims , ime , jms , jme , kms , kme , &
  1926. its , ite , jts , jte , kts , kte )
  1927. IMPLICIT NONE
  1928. INTEGER , INTENT(IN) :: num_soil_layers , &
  1929. num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  1930. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
  1931. ids , ide , jds , jde , kds , kde , &
  1932. ims , ime , jms , jme , kms , kme , &
  1933. its , ite , jts , jte , kts , kte
  1934. INTEGER , INTENT(IN) :: flag_sst, flag_soil_layers , flag_soil_levels
  1935. INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
  1936. INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
  1937. INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
  1938. REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
  1939. REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
  1940. REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
  1941. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
  1942. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: st000010, st010200, st010040 , st040100 , st100200
  1943. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: sm000010, sm010200, sm010040 , sm040100 , sm100200
  1944. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
  1945. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
  1946. REAL , DIMENSION(num_soil_layers) :: zs , dzs
  1947. REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
  1948. REAL , ALLOCATABLE , DIMENSION(:) :: zhave
  1949. INTEGER :: i , j , l , lout , lin , lwant , lhave , num
  1950. REAL :: temp
  1951. LOGICAL :: found_levels
  1952. found_levels = .false.
  1953. ! write(message, FMT='(A)') ' Prepare SSiB LSM soil fields'
  1954. !fds (06/2010)
  1955. !-- SSiB needs NCEP/NCAR reanalysis soil temperature and moisture at 0-10cm and 10-200cm --
  1956. IF ( flag_soil_layers == 1 ) THEN
  1957. DO j = jts , MIN(jde-1,jte)
  1958. DO i = its , MIN(ide-1,ite)
  1959. tslb(i,1,j) = st000010(i,j)
  1960. tslb(i,2,j) = st010040(i,j)
  1961. tslb(i,3,j) = st040100(i,j)
  1962. smois(i,1,j) = sm000010(i,j)
  1963. smois(i,2,j) = sm010040(i,j)
  1964. smois(i,3,j) = sm040100(i,j)
  1965. ! temporary fix (JD) for missing 4-layer input - use 2 instead
  1966. if(tslb(i,2,j) .lt. 10.)tslb(i,2,j) = st000010(i,j)
  1967. if(smois(i,2,j) .lt. 0.01)smois(i,2,j) = sm000010(i,j)
  1968. if(tslb(i,3,j) .lt. 10.)tslb(i,3,j) = st010200(i,j)
  1969. if(smois(i,3,j) .lt. 0.01)smois(i,3,j) = sm010200(i,j)
  1970. END DO
  1971. END DO
  1972. ELSE
  1973. CALL wrf_debug ( 0 , 'SSiB LSM needs 0-10 cm soil t')
  1974. CALL wrf_debug ( 0 , 'SSiB LSM needs 0-10 cm soil m')
  1975. CALL wrf_debug ( 0 , 'SSiB LSM needs 10-200 cm soil t')
  1976. CALL wrf_debug ( 0 , 'SSiB LSM needs 10-200 cm soil m')
  1977. ENDIF
  1978. ! Over water, put in reasonable values for soil temperature and moisture. These won't be
  1979. ! used, but they will make a more continuous plot.
  1980. IF ( flag_sst .EQ. 1 ) THEN
  1981. DO j = jts , MIN(jde-1,jte)
  1982. DO i = its , MIN(ide-1,ite)
  1983. IF ( landmask(i,j) .LT. 0.5 ) THEN
  1984. DO l = 1 , num_soil_layers
  1985. tslb(i,l,j)= sst(i,j)
  1986. smois(i,l,j)= 1.0
  1987. sh2o (i,l,j)= 1.0
  1988. END DO
  1989. END IF
  1990. END DO
  1991. END DO
  1992. ELSE
  1993. DO j = jts , MIN(jde-1,jte)
  1994. DO i = its , MIN(ide-1,ite)
  1995. IF ( landmask(i,j) .LT. 0.5 ) THEN
  1996. DO l = 1 , num_soil_layers
  1997. tslb(i,l,j)= tsk(i,j)
  1998. smois(i,l,j)= 1.0
  1999. sh2o (i,l,j)= 1.0
  2000. END DO
  2001. END IF
  2002. END DO
  2003. END DO
  2004. END IF
  2005. END SUBROUTINE init_soil_8_real
  2006. !--------------------------------------------------------
  2007. SUBROUTINE aggregate_categories_part1 ( landuse_frac , iswater , num_veg_cat , mminlu )
  2008. IMPLICIT NONE
  2009. INTEGER , INTENT(IN) :: iswater
  2010. INTEGER , INTENT(IN) :: num_veg_cat
  2011. CHARACTER (LEN=4) , INTENT(IN) :: mminlu
  2012. REAL , DIMENSION(1:num_veg_cat) , INTENT(INOUT):: landuse_frac
  2013. ! Local variables
  2014. INTEGER , PARAMETER :: num_special_bins = 3 ! grass, shrubs, trees; add a new line in the "cib" array if updating this
  2015. INTEGER , PARAMETER :: max_cats_per_bin = 8 ! larger than max num of cats inside each of the special bins
  2016. ! So, no more than 8 total tree categories that we need to consider.
  2017. INTEGER , PARAMETER :: fl = -1 ! Flag value so that we know this is not a valid category to consider.
  2018. INTEGER , DIMENSION ( max_cats_per_bin * num_special_bins ) :: cib_usgs = &
  2019. (/ 2 , 3 , 4 , 5 , 6 , 7 , 10 , fl , & ! grass
  2020. 8 , 9 , fl , fl , fl , fl , fl , fl , & ! shrubs
  2021. 11 , 12 , 13 , 14 , 15 , fl , fl , fl /) ! trees
  2022. INTEGER , DIMENSION ( max_cats_per_bin * num_special_bins ) :: cib_modis = &
  2023. (/ 9 , 10 , 12 , 14 , fl , fl , fl , fl , & ! grass
  2024. 6 , 7 , 8 , fl , fl , fl , fl , fl , & ! shrubs
  2025. 1 , 2 , 3 , 4 , 5 , fl , fl , fl /) ! trees
  2026. IF ( mminlu(1:4) .EQ. 'USGS' ) THEN
  2027. CALL aggregate_categories_part2 ( landuse_frac , iswater , num_veg_cat , &
  2028. num_special_bins , max_cats_per_bin , fl , cib_usgs )
  2029. ELSE IF ( mminlu(1:4) .EQ. 'MODI' ) THEN
  2030. CALL aggregate_categories_part2 ( landuse_frac , iswater , num_veg_cat , &
  2031. num_special_bins , max_cats_per_bin , fl , cib_modis )
  2032. END IF
  2033. END SUBROUTINE aggregate_categories_part1
  2034. SUBROUTINE aggregate_categories_part2 ( landuse_frac , iswater , num_veg_cat , &
  2035. num_special_bins , max_cats_per_bin , fl , cib )
  2036. IMPLICIT NONE
  2037. INTEGER , INTENT(IN) :: iswater
  2038. INTEGER , INTENT(IN) :: num_veg_cat
  2039. REAL , DIMENSION(1:num_veg_cat) , INTENT(INOUT):: landuse_frac
  2040. INTEGER , INTENT(IN) :: num_special_bins , max_cats_per_bin , fl
  2041. INTEGER , DIMENSION ( max_cats_per_bin * num_special_bins ) , INTENT(IN) :: cib
  2042. ! Local variables
  2043. REAL , DIMENSION(1:num_veg_cat) :: landuse_frac_work ! copy of the input array, allows us to be wreckless
  2044. INTEGER , DIMENSION ( max_cats_per_bin , num_special_bins ) :: cats_in_bin
  2045. INTEGER :: ib , ic ! indexes for the bin and the category loops
  2046. REAL , DIMENSION ( num_special_bins ) :: bin_max_val , bin_sum ! max category value in this bin, sum of all cats in this bin
  2047. INTEGER , DIMENSION ( num_special_bins ) :: bin_max_idx ! index of this maximum category in this bin
  2048. INTEGER :: bin_work , bin_orig ! the bin from whence the maximum hails, respectively for the aggregated and the original data
  2049. INTEGER :: max_cat_orig , max_cat_work
  2050. REAL :: max_val_orig , max_val_work
  2051. ! Find the max in the original. If it is >= 50%, no need to even be in here.
  2052. DO ic = 1 , num_veg_cat
  2053. IF ( landuse_frac(ic) .GE. 0.5 ) THEN
  2054. RETURN
  2055. END IF
  2056. END DO
  2057. ! Put the categories in the bin into a 2d array.
  2058. cats_in_bin = RESHAPE ( cib , (/ max_cats_per_bin , num_special_bins /) )
  2059. ! Make a copy of the incoming array so that we can eventually diff our working copy with
  2060. ! the original.
  2061. landuse_frac_work = landuse_frac
  2062. ! Loop over each of the special bins that we know about. Find the max values and the locations of such.
  2063. DO ib = 1 , num_special_bins
  2064. ! For this bin, we know about specific categories. We get the sum of all of those
  2065. ! categories that we know about for this bin, and we keep track of which category
  2066. ! is dominant.
  2067. bin_sum (ib) = 0
  2068. bin_max_val(ib) = fl
  2069. cat_loop_accumulate : DO ic = 1 , max_cats_per_bin
  2070. ! Have we run out of valid categories in this bin? For example, we would not necessarily have
  2071. ! the same number of tree categories as there are grass categories.
  2072. IF ( cats_in_bin(ic,ib) .EQ. fl ) THEN
  2073. EXIT cat_loop_accumulate
  2074. END IF
  2075. bin_sum(ib) = bin_sum(ib) + landuse_frac(cats_in_bin(ic,ib))
  2076. IF ( landuse_frac(cats_in_bin(ic,ib)) .GT. bin_max_val(ib) ) THEN
  2077. bin_max_val(ib) = landuse_frac(cats_in_bin(ic,ib))
  2078. bin_max_idx(ib) = cats_in_bin(ic,ib)
  2079. END IF
  2080. END DO cat_loop_accumulate
  2081. cat_loop_assign : DO ic = 1 , max_cats_per_bin
  2082. ! Plow through each cat in the bin. If we find the dominant one, he gets the total sum. If we land on
  2083. ! the other guys in the bin, they get set back to zero to maintain the original aggregate influence.
  2084. IF ( cats_in_bin(ic,ib) .EQ. fl ) THEN
  2085. EXIT cat_loop_assign
  2086. ELSE IF ( cats_in_bin(ic,ib) .EQ. bin_max_idx(ib) ) THEN
  2087. landuse_frac_work(cats_in_bin(ic,ib)) = bin_sum(ib)
  2088. ELSE
  2089. landuse_frac_work(cats_in_bin(ic,ib)) = 0
  2090. END IF
  2091. END DO cat_loop_assign
  2092. END DO
  2093. ! Now we loop through the categorical data, and get the max+location for the original input data and the
  2094. ! modified work data. Water is not allowed to be a "max" category unless it is greater than 50%. We hit
  2095. ! that test up at the top already, so we can toss out water willy nilly here.
  2096. max_cat_orig = fl
  2097. max_val_orig = 0
  2098. max_cat_work = fl
  2099. max_val_work = 0
  2100. DO ic = 1 , num_veg_cat
  2101. IF ( ic .EQ. iswater ) THEN
  2102. CYCLE
  2103. END IF
  2104. IF ( landuse_frac(ic) .GT. max_val_orig ) THEN
  2105. max_val_orig = landuse_frac(ic)
  2106. max_cat_orig = ic
  2107. END IF
  2108. IF ( landuse_frac_work(ic) .GT. max_val_work ) THEN
  2109. max_val_work = landuse_frac_work(ic)
  2110. max_cat_work = ic
  2111. END IF
  2112. END DO
  2113. ! Find the bin for the maximimum value of the original data.
  2114. bin_orig = -1
  2115. bin_loop_orig : DO ib = 1 , num_special_bins
  2116. cat_loop_orig : DO ic = 1 , max_cats_per_bin
  2117. IF ( cats_in_bin(ic,ib) .EQ. fl ) THEN
  2118. EXIT cat_loop_orig
  2119. ELSE IF ( cats_in_bin(ic,ib) .EQ. max_cat_orig ) THEN
  2120. bin_orig = ib
  2121. EXIT bin_loop_orig
  2122. END IF
  2123. END DO cat_loop_orig
  2124. END DO bin_loop_orig
  2125. ! Find the bin for the maximimum value of the aggregated data.
  2126. bin_work = -1
  2127. bin_loop_work : DO ib = 1 , num_special_bins
  2128. cat_loop_work : DO ic = 1 , max_cats_per_bin
  2129. IF ( cats_in_bin(ic,ib) .EQ. fl ) THEN
  2130. EXIT cat_loop_work
  2131. ELSE IF ( cats_in_bin(ic,ib) .EQ. max_cat_work ) THEN
  2132. bin_work = ib
  2133. EXIT bin_loop_work
  2134. END IF
  2135. END DO cat_loop_work
  2136. END DO bin_loop_work
  2137. ! So the big question is, did the aggregation change the bin? If the aggregation does not change the resulting
  2138. ! bin, then we leave everything alone. However, if there would be a change in the eventual bin chosen by
  2139. ! the simple dominant-category method, then we become pro-active interventionists and rectify this heinous
  2140. ! injustice foisted upon the lowly flora.
  2141. IF ( bin_work .EQ. bin_orig ) THEN
  2142. ! No op, we do nothing.
  2143. ELSE
  2144. DO ic = 1 , max_cats_per_bin
  2145. landuse_frac(cats_in_bin(ic,bin_work)) = landuse_frac_work(cats_in_bin(ic,bin_work))
  2146. END DO
  2147. END IF
  2148. END SUBROUTINE aggregate_categories_part2
  2149. END MODULE module_soil_pre
  2150. FUNCTION skip_middle_points_t ( ids , ide , jds , jde , &
  2151. i_in , j_in , width , &
  2152. subtleties_exist ) &
  2153. RESULT ( skip_it )
  2154. IMPLICIT NONE
  2155. INTEGER , INTENT(IN) :: ids , ide , jds , jde
  2156. INTEGER , INTENT(IN) :: i_in , j_in , width
  2157. LOGICAL , INTENT(IN) :: subtleties_exist
  2158. LOGICAL :: skip_it
  2159. INTEGER , PARAMETER :: slop = 0
  2160. IF ( ( subtleties_exist ) .OR. ( ide .EQ. 3 ) .OR. ( jde .EQ. 3 ) ) THEN
  2161. skip_it = .FALSE.
  2162. ELSE
  2163. IF ( ( i_in .GE. ids+width+slop ) .AND. ( i_in .LE. ide-1-width-slop ) .AND. &
  2164. ( j_in .GE. jds+width+slop ) .AND. ( j_in .LE. jde-1-width-slop ) ) THEN
  2165. skip_it = .TRUE.
  2166. ELSE
  2167. skip_it = .FALSE.
  2168. END IF
  2169. END IF
  2170. END FUNCTION skip_middle_points_t
  2171. #else
  2172. MODULE module_soil_pre
  2173. USE module_date_time
  2174. USE module_state_description
  2175. CONTAINS
  2176. SUBROUTINE adjust_for_seaice_pre ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
  2177. xland , landusef , isltyp , soilcat , soilctop , &
  2178. soilcbot , tmn , &
  2179. seaice_threshold , &
  2180. fractional_seaice, &
  2181. num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
  2182. iswater , isice , &
  2183. scheme , &
  2184. ids , ide , jds , jde , kds , kde , &
  2185. ims , ime , jms , jme , kms , kme , &
  2186. its , ite , jts , jte , kts , kte )
  2187. IMPLICIT NONE
  2188. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  2189. ims , ime , jms , jme , kms , kme , &
  2190. its , ite , jts , jte , kts , kte , &
  2191. iswater , isice
  2192. INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
  2193. REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
  2194. REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
  2195. REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
  2196. INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
  2197. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
  2198. vegcat, xland , soilcat , tmn
  2199. REAL , INTENT(IN) :: seaice_threshold
  2200. INTEGER :: i , j , num_seaice_changes , loop
  2201. CHARACTER (LEN=132) :: message
  2202. integer, intent(in) :: fractional_seaice
  2203. real :: xice_threshold
  2204. IF ( FRACTIONAL_SEAICE == 0 ) THEN
  2205. xice_threshold = 0.5
  2206. ELSEIF ( FRACTIONAL_SEAICE == 1 ) THEN
  2207. CALL WRF_ERROR_FATAL("NMM cannot use FRACTIONAL_SEAICE = 1")
  2208. xice_threshold = 0.02
  2209. ENDIF
  2210. fix_seaice : SELECT CASE ( scheme )
  2211. CASE ( SLABSCHEME )
  2212. DO j = jts , MIN(jde-1,jte)
  2213. DO i = its , MIN(ide-1,ite)
  2214. IF ( xice(i,j) .GT. 200.0 ) THEN
  2215. xice(i,j) = 0.
  2216. num_seaice_changes = num_seaice_changes + 1
  2217. END IF
  2218. END DO
  2219. END DO
  2220. IF ( num_seaice_changes .GT. 0 ) THEN
  2221. WRITE ( message , FMT='(A,I6)' ) &
  2222. 'Total pre number of sea ice locations removed (due to FLAG values) = ', &
  2223. num_seaice_changes
  2224. CALL wrf_debug ( 0 , message )
  2225. END IF
  2226. num_seaice_changes = 0
  2227. DO j = jts , MIN(jde-1,jte)
  2228. DO i = its , MIN(ide-1,ite)
  2229. IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
  2230. ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
  2231. IF ( FRACTIONAL_SEAICE == 0 ) THEN
  2232. xice(i,j) = 1.0
  2233. ENDIF
  2234. num_seaice_changes = num_seaice_changes + 1
  2235. tmn(i,j) = 271.4
  2236. vegcat(i,j)=isice
  2237. lu_index(i,j)=ivgtyp(i,j)
  2238. landmask(i,j)=1.
  2239. xland(i,j)=1.
  2240. DO loop=1,num_veg_cat
  2241. landusef(i,loop,j)=0.
  2242. END DO
  2243. landusef(i,ivgtyp(i,j),j)=1.
  2244. isltyp(i,j) = 16
  2245. soilcat(i,j)=isltyp(i,j)
  2246. DO loop=1,num_soil_top_cat
  2247. soilctop(i,loop,j)=0
  2248. END DO
  2249. DO loop=1,num_soil_bot_cat
  2250. soilcbot(i,loop,j)=0
  2251. END DO
  2252. soilctop(i,isltyp(i,j),j)=1.
  2253. soilcbot(i,isltyp(i,j),j)=1.
  2254. ELSE
  2255. xice(i,j) = 0.0
  2256. END IF
  2257. END DO
  2258. END DO
  2259. IF ( num_seaice_changes .GT. 0 ) THEN
  2260. WRITE ( message , FMT='(A,I6)' ) &
  2261. 'Total pre number of sea ice location changes (water to land) = ', num_seaice_changes
  2262. CALL wrf_debug ( 0 , message )
  2263. END IF
  2264. CASE ( LSMSCHEME )
  2265. CASE ( RUCLSMSCHEME )
  2266. END SELECT fix_seaice
  2267. END SUBROUTINE adjust_for_seaice_pre
  2268. SUBROUTINE adjust_for_seaice_post ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
  2269. xland , landusef , isltyp , soilcat , soilctop , &
  2270. soilcbot , tmn , &
  2271. tslb , smois , sh2o , &
  2272. seaice_threshold , &
  2273. fractional_seaice, &
  2274. num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
  2275. num_soil_layers , &
  2276. iswater , isice , &
  2277. scheme , &
  2278. ids , ide , jds , jde , kds , kde , &
  2279. ims , ime , jms , jme , kms , kme , &
  2280. its , ite , jts , jte , kts , kte )
  2281. IMPLICIT NONE
  2282. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  2283. ims , ime , jms , jme , kms , kme , &
  2284. its , ite , jts , jte , kts , kte , &
  2285. iswater , isice
  2286. INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
  2287. INTEGER , INTENT(IN) :: num_soil_layers
  2288. REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
  2289. REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
  2290. REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
  2291. REAL , DIMENSION(ims:ime,1:num_soil_layers,jms:jme) , INTENT(INOUT):: tslb , smois , sh2o
  2292. INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
  2293. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
  2294. vegcat, xland , soilcat , tmn
  2295. REAL , INTENT(IN) :: seaice_threshold
  2296. REAL :: total_depth , mid_point_depth
  2297. INTEGER :: i , j , num_seaice_changes , loop
  2298. CHARACTER (LEN=132) :: message
  2299. integer, intent(in) :: fractional_seaice
  2300. real :: xice_threshold
  2301. IF ( FRACTIONAL_SEAICE == 0 ) THEN
  2302. xice_threshold = 0.5
  2303. ELSEIF ( FRACTIONAL_SEAICE == 1 ) THEN
  2304. CALL WRF_ERROR_FATAL("NMM cannot use FRACTIONAL_SEAICE = 1")
  2305. xice_threshold = 0.02
  2306. ENDIF
  2307. fix_seaice : SELECT CASE ( scheme )
  2308. CASE ( SLABSCHEME )
  2309. !SBAO
  2310. CASE ( LSMSCHEME, GFDLSLAB )
  2311. DO j = jts , MIN(jde-1,jte)
  2312. DO i = its , MIN(ide-1,ite)
  2313. IF ( xice(i,j) .GT. 200.0 ) THEN
  2314. xice(i,j) = 0.
  2315. num_seaice_changes = num_seaice_changes + 1
  2316. END IF
  2317. END DO
  2318. END DO
  2319. IF ( num_seaice_changes .GT. 0 ) THEN
  2320. WRITE ( message , FMT='(A,I6)' ) &
  2321. 'Total post number of sea ice locations removed (due to FLAG values) = ', &
  2322. num_seaice_changes
  2323. CALL wrf_debug ( 0 , message )
  2324. END IF
  2325. num_seaice_changes = 0
  2326. DO j = jts , MIN(jde-1,jte)
  2327. DO i = its , MIN(ide-1,ite)
  2328. IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
  2329. ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
  2330. IF ( FRACTIONAL_SEAICE == 0 ) THEN
  2331. xice(i,j) = 1.0
  2332. ENDIF
  2333. num_seaice_changes = num_seaice_changes + 1
  2334. tmn(i,j) = 271.16
  2335. vegcat(i,j)=isice
  2336. lu_index(i,j)=ivgtyp(i,j)
  2337. landmask(i,j)=1.
  2338. xland(i,j)=1.
  2339. DO loop=1,num_veg_cat
  2340. landusef(i,loop,j)=0.
  2341. END DO
  2342. landusef(i,ivgtyp(i,j),j)=1.
  2343. isltyp(i,j) = 16
  2344. soilcat(i,j)=isltyp(i,j)
  2345. DO loop=1,num_soil_top_cat
  2346. soilctop(i,loop,j)=0
  2347. END DO
  2348. DO loop=1,num_soil_bot_cat
  2349. soilcbot(i,loop,j)=0
  2350. END DO
  2351. soilctop(i,isltyp(i,j),j)=1.
  2352. soilcbot(i,isltyp(i,j),j)=1.
  2353. total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
  2354. DO loop = 1,num_soil_layers
  2355. mid_point_depth=(total_depth/num_soil_layers)/2. + &
  2356. (loop-1)*(total_depth/num_soil_layers)
  2357. tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
  2358. mid_point_depth*tmn(i,j) ) / total_depth
  2359. END DO
  2360. DO loop=1,num_soil_layers
  2361. smois(i,loop,j) = 1.0
  2362. sh2o(i,loop,j) = 0.0
  2363. END DO
  2364. ELSE
  2365. xice(i,j) = 0.0
  2366. END IF
  2367. END DO
  2368. END DO
  2369. IF ( num_seaice_changes .GT. 0 ) THEN
  2370. WRITE ( message , FMT='(A,I6)' ) &
  2371. 'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
  2372. CALL wrf_debug ( 0 , message )
  2373. END IF
  2374. CASE ( RUCLSMSCHEME )
  2375. DO j = jts , MIN(jde-1,jte)
  2376. DO i = its , MIN(ide-1,ite)
  2377. IF ( xice(i,j) .GT. 200.0 ) THEN
  2378. xice(i,j) = 0.
  2379. num_seaice_changes = num_seaice_changes + 1
  2380. END IF
  2381. END DO
  2382. END DO
  2383. IF ( num_seaice_changes .GT. 0 ) THEN
  2384. WRITE ( message , FMT='(A,I6)' ) &
  2385. 'Total post number of sea ice locations removed (due to FLAG values) = ', &
  2386. num_seaice_changes
  2387. CALL wrf_debug ( 0 , message )
  2388. END IF
  2389. num_seaice_changes = 0
  2390. DO j = jts , MIN(jde-1,jte)
  2391. DO i = its , MIN(ide-1,ite)
  2392. IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
  2393. ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
  2394. IF ( FRACTIONAL_SEAICE == 0 ) THEN
  2395. xice(i,j) = 1.0
  2396. ELSE
  2397. xice(i,j)=max(0.25,xice(i,j))
  2398. ENDIF
  2399. num_seaice_changes = num_seaice_changes + 1
  2400. if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
  2401. vegcat(i,j)=isice
  2402. ivgtyp(i,j)=isice
  2403. lu_index(i,j)=isice
  2404. landmask(i,j)=1.
  2405. xland(i,j)=1.
  2406. DO loop=1,num_veg_cat
  2407. landusef(i,loop,j)=0.
  2408. END DO
  2409. landusef(i,ivgtyp(i,j),j)=1.
  2410. isltyp(i,j) = 16
  2411. soilcat(i,j)=isltyp(i,j)
  2412. DO loop=1,num_soil_top_cat
  2413. soilctop(i,loop,j)=0
  2414. END DO
  2415. DO loop=1,num_soil_bot_cat
  2416. soilcbot(i,loop,j)=0
  2417. END DO
  2418. soilctop(i,isltyp(i,j),j)=1.
  2419. soilcbot(i,isltyp(i,j),j)=1.
  2420. total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
  2421. tslb(i,1,j) = tsk(i,j)
  2422. tslb(i,num_soil_layers,j) = tmn(i,j)
  2423. DO loop = 2,num_soil_layers-1
  2424. mid_point_depth=(total_depth/num_soil_layers)/4. + &
  2425. (loop-2)*(total_depth/num_soil_layers)
  2426. tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
  2427. mid_point_depth*tmn(i,j) ) / total_depth
  2428. END DO
  2429. DO loop=1,num_soil_layers
  2430. smois(i,loop,j) = 1.0
  2431. sh2o(i,loop,j) = 0.0
  2432. END DO
  2433. ELSE IF ( xice(i,j) .LT. xice_threshold ) THEN
  2434. xice(i,j) = 0.
  2435. END IF
  2436. END DO
  2437. END DO
  2438. IF ( num_seaice_changes .GT. 0 ) THEN
  2439. WRITE ( message , FMT='(A,I6)' ) &
  2440. 'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
  2441. CALL wrf_debug ( 0 , message )
  2442. END IF
  2443. END SELECT fix_seaice
  2444. END SUBROUTINE adjust_for_seaice_post
  2445. SUBROUTINE process_percent_cat_new ( landmask , &
  2446. landuse_frac , soil_top_cat , soil_bot_cat , &
  2447. isltyp , ivgtyp , &
  2448. num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
  2449. ids , ide , jds , jde , kds , kde , &
  2450. ims , ime , jms , jme , kms , kme , &
  2451. its , ite , jts , jte , kts , kte , &
  2452. iswater )
  2453. IMPLICIT NONE
  2454. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  2455. ims , ime , jms , jme , kms , kme , &
  2456. its , ite , jts , jte , kts , kte , &
  2457. iswater
  2458. INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
  2459. REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landuse_frac
  2460. REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(IN):: soil_top_cat
  2461. REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(IN):: soil_bot_cat
  2462. INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
  2463. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask
  2464. INTEGER :: i , j , l , ll, dominant_index
  2465. REAL :: dominant_value
  2466. #ifdef WRF_CHEM
  2467. ! REAL :: lwthresh = .99
  2468. REAL :: lwthresh = .50
  2469. #else
  2470. REAL :: lwthresh = .50
  2471. #endif
  2472. INTEGER , PARAMETER :: iswater_soil = 14
  2473. INTEGER :: iforce
  2474. CHARACTER (LEN=1024) :: message
  2475. ! Sanity check on the 50/50 points
  2476. DO j = jts , MIN(jde-1,jte)
  2477. DO i = its , MIN(ide-1,ite)
  2478. dominant_value = landuse_frac(i,iswater,j)
  2479. IF ( dominant_value .EQ. lwthresh ) THEN
  2480. DO l = 1 , num_veg_cat
  2481. IF ( l .EQ. iswater ) CYCLE
  2482. IF ( landuse_frac(i,l,j) .EQ. lwthresh ) THEN
  2483. write(message, FMT='(I4,I4,A,I4,A,F12.4)') i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
  2484. call wrf_message ( message )
  2485. landuse_frac(i,l,j) = lwthresh - .01
  2486. landuse_frac(i,iswater,j) = 1. - (lwthresh + 0.01)
  2487. END IF
  2488. END DO
  2489. END IF
  2490. END DO
  2491. END DO
  2492. ! Compute the dominant VEGETATION INDEX.
  2493. DO j = jts , MIN(jde-1,jte)
  2494. DO i = its , MIN(ide-1,ite)
  2495. dominant_value = landuse_frac(i,1,j)
  2496. dominant_index = 1
  2497. DO l = 2 , num_veg_cat
  2498. IF ( l .EQ. iswater ) THEN
  2499. ! wait a bit
  2500. ELSE IF ( ( l .NE. iswater ) .AND. ( landuse_frac(i,l,j) .GT. dominant_value ) ) THEN
  2501. dominant_value = landuse_frac(i,l,j)
  2502. dominant_index = l
  2503. END IF
  2504. END DO
  2505. IF ( landuse_frac(i,iswater,j) .GT. lwthresh ) THEN
  2506. dominant_value = landuse_frac(i,iswater,j)
  2507. dominant_index = iswater
  2508. #if 0
  2509. si needs to beef up consistency checks before we can use this part
  2510. ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
  2511. ( dominant_value .EQ. lwthresh) ) THEN
  2512. ! no op
  2513. #else
  2514. ELSEIF((landuse_frac(i,iswater,j).EQ.lwthresh).AND.(dominant_value.EQ.lwthresh).and.(dominant_index.LT.iswater))THEN
  2515. write(message,*)'temporary SI landmask fix'
  2516. call wrf_message(trim(message))
  2517. ! no op
  2518. ELSEIF((landuse_frac(i,iswater,j).EQ.lwthresh).AND.(dominant_value.EQ.lwthresh).and.(dominant_index.GT.iswater))THEN
  2519. write(message,*)'temporary SI landmask fix'
  2520. call wrf_message(trim(message))
  2521. dominant_value = landuse_frac(i,iswater,j)
  2522. dominant_index = iswater
  2523. #endif
  2524. ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh ) .AND. &
  2525. ( dominant_value .LT. lwthresh ) ) THEN
  2526. dominant_value = landuse_frac(i,iswater,j)
  2527. dominant_index = iswater
  2528. END IF
  2529. IF ( dominant_index .EQ. iswater ) THEN
  2530. if(landmask(i,j).gt.lwthresh) then
  2531. write(message,*) 'changing to water at point ',i,j
  2532. call wrf_message(trim(message))
  2533. write(message,*) nint(landuse_frac(i,:,j)*100)
  2534. call wrf_message(trim(message))
  2535. endif
  2536. landmask(i,j) = 0
  2537. ELSE IF ( dominant_index .NE. iswater ) THEN
  2538. if(landmask(i,j).lt.lwthresh) then
  2539. write(message,*) 'changing to land at point ',i,j
  2540. call wrf_message(trim(message))
  2541. write(message,*) nint(landuse_frac(i,:,j)*100)
  2542. call wrf_message(trim(message))
  2543. endif
  2544. landmask(i,j) = 1
  2545. END IF
  2546. ivgtyp(i,j) = dominant_index
  2547. END DO
  2548. END DO
  2549. ! Compute the dominant SOIL TEXTURE INDEX, TOP.
  2550. iforce = 0
  2551. DO i = its , MIN(ide-1,ite)
  2552. DO j = jts , MIN(jde-1,jte)
  2553. dominant_value = soil_top_cat(i,1,j)
  2554. dominant_index = 1
  2555. IF ( landmask(i,j) .GT. lwthresh ) THEN
  2556. DO l = 2 , num_soil_top_cat
  2557. IF ( ( l .NE. iswater_soil ) .AND. ( soil_top_cat(i,l,j) .GT. dominant_value ) ) THEN
  2558. dominant_value = soil_top_cat(i,l,j)
  2559. dominant_index = l
  2560. END IF
  2561. END DO
  2562. IF ( dominant_value .LT. 0.01 ) THEN
  2563. iforce = iforce + 1
  2564. WRITE ( message , FMT = '(A,I4,I4)' ) &
  2565. 'based on landuse, changing soil to land at point ',i,j
  2566. CALL wrf_debug(1,message)
  2567. !atec WRITE ( message , FMT = '(16(i3,1x))' ) &
  2568. !atec 1, 2, 3, 4, 5, 6, 7, 8, 9, 10,11,12, 13, 14, 15, 16
  2569. !atec CALL wrf_debug(1,message)
  2570. !atec WRITE ( message , FMT = '(16(i3,1x))' ) &
  2571. !atec nint(soil_top_cat(i,:,j)*100)
  2572. !atec CALL wrf_debug(1,message)
  2573. dominant_index = 8
  2574. END IF
  2575. ELSE
  2576. dominant_index = iswater_soil
  2577. END IF
  2578. isltyp(i,j) = dominant_index
  2579. END DO
  2580. END DO
  2581. if(iforce.ne.0)then
  2582. WRITE(message,FMT='(A,I4,A,I6)' ) &
  2583. 'forcing artificial silty clay loam at ',iforce,' points, out of ',&
  2584. (MIN(ide-1,ite)-its+1)*(MIN(jde-1,jte)-jts+1)
  2585. CALL wrf_debug(0,message)
  2586. endif
  2587. END SUBROUTINE process_percent_cat_new
  2588. SUBROUTINE process_soil_real ( tsk , tmn , &
  2589. landmask , sst , &
  2590. st_input , sm_input , sw_input , st_levels_input , sm_levels_input , sw_levels_input , &
  2591. zs , dzs , tslb , smois , sh2o , &
  2592. flag_sst , flag_soilt000, flag_soilm000, &
  2593. ids , ide , jds , jde , kds , kde , &
  2594. ims , ime , jms , jme , kms , kme , &
  2595. its , ite , jts , jte , kts , kte , &
  2596. sf_surface_physics , num_soil_layers , real_data_init_type , &
  2597. num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  2598. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
  2599. IMPLICIT NONE
  2600. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  2601. ims , ime , jms , jme , kms , kme , &
  2602. its , ite , jts , jte , kts , kte , &
  2603. sf_surface_physics , num_soil_layers , real_data_init_type , &
  2604. num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  2605. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc
  2606. INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
  2607. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
  2608. INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
  2609. INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
  2610. INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
  2611. REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
  2612. REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
  2613. REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
  2614. REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
  2615. REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
  2616. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
  2617. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
  2618. INTEGER :: i , j , l , dominant_index , num_soil_cat , num_veg_cat
  2619. REAL :: dominant_value
  2620. ! Initialize the soil depth, and the soil temperature and moisture.
  2621. IF ( ( sf_surface_physics .EQ. SLABSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  2622. CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
  2623. CALL init_soil_1_real ( tsk , tmn , tslb , zs , dzs , num_soil_layers , real_data_init_type , &
  2624. landmask , sst , flag_sst , &
  2625. ids , ide , jds , jde , kds , kde , &
  2626. ims , ime , jms , jme , kms , kme , &
  2627. its , ite , jts , jte , kts , kte )
  2628. ELSE IF ( ( sf_surface_physics .EQ. LSMSCHEME .or. sf_surface_physics .eq. 88 ) .AND. &
  2629. ( num_soil_layers .GT. 1 ) ) THEN
  2630. CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
  2631. CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
  2632. st_input , sm_input , sw_input , landmask , sst , &
  2633. zs , dzs , &
  2634. st_levels_input , sm_levels_input , sw_levels_input , &
  2635. num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  2636. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
  2637. flag_sst , flag_soilt000 , flag_soilm000 , &
  2638. ids , ide , jds , jde , kds , kde , &
  2639. ims , ime , jms , jme , kms , kme , &
  2640. its , ite , jts , jte , kts , kte )
  2641. ! CALL init_soil_old ( tsk , tmn , &
  2642. ! smois , tslb , zs , dzs , num_soil_layers , &
  2643. ! st000010_input , st010040_input , st040100_input , st100200_input , &
  2644. ! st010200_input , &
  2645. ! sm000010_input , sm010040_input , sm040100_input , sm100200_input , &
  2646. ! sm010200_input , &
  2647. ! landmask_input , sst_input , &
  2648. ! ids , ide , jds , jde , kds , kde , &
  2649. ! ims , ime , jms , jme , kms , kme , &
  2650. ! its , ite , jts , jte , kts , kte )
  2651. ELSE IF ( ( sf_surface_physics .EQ. NOAHMPSCHEME .or. sf_surface_physics .eq. 88 ) .AND. &
  2652. ( num_soil_layers .GT. 1 ) ) THEN
  2653. CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
  2654. CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
  2655. st_input , sm_input , sw_input , landmask , sst , &
  2656. zs , dzs , &
  2657. st_levels_input , sm_levels_input , sw_levels_input , &
  2658. num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  2659. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
  2660. flag_sst , flag_soilt000 , flag_soilm000 , &
  2661. ids , ide , jds , jde , kds , kde , &
  2662. ims , ime , jms , jme , kms , kme , &
  2663. its , ite , jts , jte , kts , kte )
  2664. ! CALL init_soil_old ( tsk , tmn , &
  2665. ! smois , tslb , zs , dzs , num_soil_layers , &
  2666. ! st000010_input , st010040_input , st040100_input , st100200_input , &
  2667. ! st010200_input , &
  2668. ! sm000010_input , sm010040_input , sm040100_input , sm100200_input , &
  2669. ! sm010200_input , &
  2670. ! landmask_input , sst_input , &
  2671. ! ids , ide , jds , jde , kds , kde , &
  2672. ! ims , ime , jms , jme , kms , kme , &
  2673. ! its , ite , jts , jte , kts , kte )
  2674. ELSE IF ( ( sf_surface_physics .EQ. RUCLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  2675. CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
  2676. CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
  2677. st_input , sm_input , landmask , sst , &
  2678. zs , dzs , &
  2679. st_levels_input , sm_levels_input , &
  2680. num_soil_layers , num_st_levels_input , num_sm_levels_input , &
  2681. num_st_levels_alloc , num_sm_levels_alloc , &
  2682. flag_sst , flag_soilt000 , flag_soilm000 , &
  2683. ids , ide , jds , jde , kds , kde , &
  2684. ims , ime , jms , jme , kms , kme , &
  2685. its , ite , jts , jte , kts , kte )
  2686. END IF
  2687. END SUBROUTINE process_soil_real
  2688. SUBROUTINE process_soil_ideal ( xland,xice,vegfra,snow,canwat, &
  2689. ivgtyp,isltyp,tslb,smois, &
  2690. tsk,tmn,zs,dzs, &
  2691. num_soil_layers, &
  2692. sf_surface_physics , &
  2693. ids,ide, jds,jde, kds,kde,&
  2694. ims,ime, jms,jme, kms,kme,&
  2695. its,ite, jts,jte, kts,kte )
  2696. IMPLICIT NONE
  2697. INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
  2698. ims,ime, jms,jme, kms,kme, &
  2699. its,ite, jts,jte, kts,kte
  2700. INTEGER, INTENT(IN) :: num_soil_layers , sf_surface_physics
  2701. REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(INOUT) :: smois, tslb
  2702. REAL, DIMENSION(num_soil_layers), INTENT(OUT) :: dzs,zs
  2703. REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(INOUT) :: tsk, tmn
  2704. REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra
  2705. INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
  2706. ! Local variables.
  2707. INTEGER :: itf,jtf
  2708. itf=MIN(ite,ide-1)
  2709. jtf=MIN(jte,jde-1)
  2710. IF ( ( sf_surface_physics .EQ. SLABSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  2711. CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
  2712. CALL init_soil_1_ideal(tsk,tmn,tslb,xland, &
  2713. ivgtyp,zs,dzs,num_soil_layers, &
  2714. ids,ide, jds,jde, kds,kde, &
  2715. ims,ime, jms,jme, kms,kme, &
  2716. its,ite, jts,jte, kts,kte )
  2717. ELSE IF ( ( sf_surface_physics .EQ. LSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  2718. CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
  2719. CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
  2720. ivgtyp,isltyp,tslb,smois,tmn, &
  2721. num_soil_layers, &
  2722. ids,ide, jds,jde, kds,kde, &
  2723. ims,ime, jms,jme, kms,kme, &
  2724. its,ite, jts,jte, kts,kte )
  2725. ELSE IF ( ( sf_surface_physics .EQ. NOAHMPSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  2726. CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
  2727. CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
  2728. ivgtyp,isltyp,tslb,smois,tmn, &
  2729. num_soil_layers, &
  2730. ids,ide, jds,jde, kds,kde, &
  2731. ims,ime, jms,jme, kms,kme, &
  2732. its,ite, jts,jte, kts,kte )
  2733. ELSE IF ( ( sf_surface_physics .EQ. RUCLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
  2734. CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
  2735. END IF
  2736. END SUBROUTINE process_soil_ideal
  2737. SUBROUTINE adjust_soil_temp_new ( tmn , sf_surface_physics , &
  2738. tsk , ter , toposoil , landmask , flag_toposoil , &
  2739. st000010 , st010040 , st040100 , st100200 , st010200 , &
  2740. flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200 , &
  2741. soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300 , &
  2742. flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300 , &
  2743. ids , ide , jds , jde , kds , kde , &
  2744. ims , ime , jms , jme , kms , kme , &
  2745. its , ite , jts , jte , kts , kte )
  2746. IMPLICIT NONE
  2747. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  2748. ims , ime , jms , jme , kms , kme , &
  2749. its , ite , jts , jte , kts , kte
  2750. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: ter , toposoil , landmask
  2751. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tmn , tsk
  2752. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: st000010 , st010040 , st040100 , st100200 , st010200
  2753. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: soilt000 , soilt005 , soilt020 , soilt040 , soilt160 , soilt300
  2754. INTEGER , INTENT(IN) :: flag_st000010 , flag_st010040 , flag_st040100 , flag_st100200 , flag_st010200
  2755. INTEGER , INTENT(IN) :: flag_soilt000 , flag_soilt005 , flag_soilt020 , flag_soilt040 , flag_soilt160 , flag_soilt300
  2756. INTEGER , INTENT(IN) :: sf_surface_physics , flag_toposoil
  2757. INTEGER :: i , j
  2758. REAL :: soil_elev_min_val , soil_elev_max_val , soil_elev_min_dif , soil_elev_max_dif
  2759. ! Do we have a soil field with which to modify soil temperatures?
  2760. IF ( flag_toposoil .EQ. 1 ) THEN
  2761. DO j = jts , MIN(jde-1,jte)
  2762. DO i = its , MIN(ide-1,ite)
  2763. ! Is the toposoil field OK, or is it a subversive soil elevation field. We can tell
  2764. ! usually by looking at values. Anything less than -1000 m (lower than the Dead Sea) is
  2765. ! bad. Anything larger than 10 km (taller than Everest) is toast. Also, anything where
  2766. ! the difference between the soil elevation and the terrain is greater than 3 km means
  2767. ! that the soil data is either all zeros or that the data are inconsistent. Any of these
  2768. ! three conditions is grievous enough to induce a WRF fatality. However, if they are at
  2769. ! a water point, then we can safely ignore them.
  2770. soil_elev_min_val = toposoil(i,j)
  2771. soil_elev_max_val = toposoil(i,j)
  2772. soil_elev_min_dif = ter(i,j) - toposoil(i,j)
  2773. soil_elev_max_dif = ter(i,j) - toposoil(i,j)
  2774. IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
  2775. CYCLE
  2776. ELSE IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
  2777. !print *,'no soil temperature elevation adjustment, soil height too small = ',toposoil(i,j)
  2778. cycle
  2779. ! CALL wrf_error_fatal ( 'TOPOSOIL values have large negative values < -1000 m, unrealistic.' )
  2780. ENDIF
  2781. IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
  2782. CYCLE
  2783. ELSE IF ( ( soil_elev_max_val .GT. 10000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
  2784. print *,'no soil temperature elevation adjustment, soil height too high = ',toposoil(i,j)
  2785. cycle
  2786. CALL wrf_error_fatal ( 'TOPOSOIL values have large positive values > 10,000 m , unrealistic.' )
  2787. ENDIF
  2788. IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
  2789. ( landmask(i,j) .LT. 0.5 ) ) THEN
  2790. CYCLE
  2791. ELSE IF ( ( ( soil_elev_min_dif .LT. -3000 ) .OR. ( soil_elev_max_dif .GT. 3000 ) ) .AND. &
  2792. ( landmask(i,j) .GT. 0.5 ) ) THEN
  2793. print *,'no soil temperature elevation adjustment, diff of soil height and terrain = ',ter(i,j) - toposoil(i,j)
  2794. cycle
  2795. CALL wrf_error_fatal ( 'TOPOSOIL difference with terrain elevation differs by more than 3000 m, unrealistic' )
  2796. ENDIF
  2797. ! For each of the fields that we would like to modify, check to see if it came in from the SI.
  2798. ! If so, then use a -6.5 K/km lapse rate (based on the elevation diffs). We only adjust when we
  2799. ! are not at a water point.
  2800. IF (landmask(i,j) .GT. 0.5 ) THEN
  2801. IF ( sf_surface_physics .EQ. SLABSCHEME .OR. sf_surface_physics .EQ. PXLSMSCHEME ) THEN
  2802. tmn(i,j) = tmn(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2803. END IF
  2804. tsk(i,j) = tsk(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2805. IF ( flag_st000010 .EQ. 1 ) THEN
  2806. st000010(i,j) = st000010(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2807. END IF
  2808. IF ( flag_st010040 .EQ. 1 ) THEN
  2809. st010040(i,j) = st010040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2810. END IF
  2811. IF ( flag_st040100 .EQ. 1 ) THEN
  2812. st040100(i,j) = st040100(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2813. END IF
  2814. IF ( flag_st100200 .EQ. 1 ) THEN
  2815. st100200(i,j) = st100200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2816. END IF
  2817. IF ( flag_st010200 .EQ. 1 ) THEN
  2818. st010200(i,j) = st010200(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2819. END IF
  2820. IF ( flag_soilt000 .EQ. 1 ) THEN
  2821. soilt000(i,j) = soilt000(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2822. END IF
  2823. IF ( flag_soilt005 .EQ. 1 ) THEN
  2824. soilt005(i,j) = soilt005(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2825. END IF
  2826. IF ( flag_soilt020 .EQ. 1 ) THEN
  2827. soilt020(i,j) = soilt020(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2828. END IF
  2829. IF ( flag_soilt040 .EQ. 1 ) THEN
  2830. soilt040(i,j) = soilt040(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2831. END IF
  2832. IF ( flag_soilt160 .EQ. 1 ) THEN
  2833. soilt160(i,j) = soilt160(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2834. END IF
  2835. IF ( flag_soilt300 .EQ. 1 ) THEN
  2836. soilt300(i,j) = soilt300(i,j) - 0.0065 * ( ter(i,j) - toposoil(i,j) )
  2837. END IF
  2838. END IF
  2839. END DO
  2840. END DO
  2841. END IF
  2842. END SUBROUTINE adjust_soil_temp_new
  2843. SUBROUTINE init_soil_depth_1 ( zs , dzs , num_soil_layers )
  2844. IMPLICIT NONE
  2845. INTEGER, INTENT(IN) :: num_soil_layers
  2846. REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
  2847. INTEGER :: l
  2848. CHARACTER (LEN=132) :: message
  2849. ! Define layers (top layer = 0.01 m). Double the thicknesses at each step (dzs values).
  2850. ! The distance from the ground level to the midpoint of the layer is given by zs.
  2851. ! ------- Ground Level ---------- || || || ||
  2852. ! || || || || zs(1) = 0.005 m
  2853. ! -- -- -- -- -- -- -- -- -- || || || \/
  2854. ! || || ||
  2855. ! ----------------------------------- || || || \/ dzs(1) = 0.01 m
  2856. ! || || ||
  2857. ! || || || zs(2) = 0.02
  2858. ! -- -- -- -- -- -- -- -- -- || || \/
  2859. ! || ||
  2860. ! || ||
  2861. ! ----------------------------------- || || \/ dzs(2) = 0.02 m
  2862. ! || ||
  2863. ! || ||
  2864. ! || ||
  2865. ! || || zs(3) = 0.05
  2866. ! -- -- -- -- -- -- -- -- -- || \/
  2867. ! ||
  2868. ! ||
  2869. ! ||
  2870. ! ||
  2871. ! ----------------------------------- \/ dzs(3) = 0.04 m
  2872. IF ( num_soil_layers .NE. 5 ) THEN
  2873. write(message,FMT= '(A)') 'Usually, the 5-layer diffusion uses 5 layers. Change this in the namelist.'
  2874. CALL wrf_error_fatal ( message )
  2875. END IF
  2876. dzs(1)=.01
  2877. zs(1)=.5*dzs(1)
  2878. DO l=2,num_soil_layers
  2879. dzs(l)=2*dzs(l-1)
  2880. zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
  2881. ENDDO
  2882. END SUBROUTINE init_soil_depth_1
  2883. SUBROUTINE init_soil_depth_2 ( zs , dzs , num_soil_layers )
  2884. IMPLICIT NONE
  2885. INTEGER, INTENT(IN) :: num_soil_layers
  2886. REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
  2887. INTEGER :: l
  2888. CHARACTER (LEN=132) :: message
  2889. dzs = (/ 0.1 , 0.3 , 0.6 , 1.0 /)
  2890. IF ( num_soil_layers .NE. 4 ) THEN
  2891. write(message,FMT='(A)') 'Usually, the LSM uses 4 layers. Change this in the namelist.'
  2892. CALL wrf_error_fatal ( message )
  2893. END IF
  2894. zs(1)=.5*dzs(1)
  2895. DO l=2,num_soil_layers
  2896. zs(l)=zs(l-1)+.5*dzs(l-1)+.5*dzs(l)
  2897. ENDDO
  2898. END SUBROUTINE init_soil_depth_2
  2899. SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_layers )
  2900. IMPLICIT NONE
  2901. INTEGER, INTENT(IN) :: num_soil_layers
  2902. REAL, DIMENSION(1:num_soil_layers), INTENT(OUT) :: zs,dzs
  2903. INTEGER :: l
  2904. CHARACTER (LEN=132) :: message
  2905. ! in RUC LSM ZS - soil levels, and DZS - soil layer thicknesses, not used
  2906. ! ZS is specified in the namelist: num_soil_layers = 6 or 9.
  2907. ! Other options with number of levels are possible, but
  2908. ! WRF users should change consistently the namelist entry with the
  2909. ! ZS array in this subroutine.
  2910. IF ( num_soil_layers .EQ. 6) THEN
  2911. zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 1.60 , 3.00 /)
  2912. ! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
  2913. ELSEIF ( num_soil_layers .EQ. 9) THEN
  2914. zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 0.60, 1.00, 1.60 , 2.20, 3.00 /)
  2915. ! dzs = (/ 0.00 , 0.125, 0.175 , 0.70 , 1.30 , 1.40 /)
  2916. ENDIF
  2917. IF ( num_soil_layers .EQ. 4 .OR. num_soil_layers .EQ. 5 ) THEN
  2918. WRITE(message,FMT= '(A)')'Usually, the RUC LSM uses 6, 9 or more levels. Change this in the namelist.'
  2919. CALL wrf_error_fatal ( message )
  2920. END IF
  2921. END SUBROUTINE init_soil_depth_3
  2922. SUBROUTINE init_soil_1_real ( tsk , tmn , tslb , zs , dzs , &
  2923. num_soil_layers , real_data_init_type , &
  2924. landmask , sst , flag_sst , &
  2925. ids , ide , jds , jde , kds , kde , &
  2926. ims , ime , jms , jme , kms , kme , &
  2927. its , ite , jts , jte , kts , kte )
  2928. IMPLICIT NONE
  2929. INTEGER , INTENT(IN) :: num_soil_layers , real_data_init_type , &
  2930. ids , ide , jds , jde , kds , kde , &
  2931. ims , ime , jms , jme , kms , kme , &
  2932. its , ite , jts , jte , kts , kte
  2933. INTEGER , INTENT(IN) :: flag_sst
  2934. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
  2935. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tsk , tmn
  2936. REAL , DIMENSION(num_soil_layers) :: zs , dzs
  2937. REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb
  2938. INTEGER :: i , j , l
  2939. ! Soil temperature is linearly interpolated between the skin temperature (taken to be at a
  2940. ! depth of 0.5 cm) and the deep soil, annual temperature (taken to be at a depth of 23 cm).
  2941. ! The tslb(i,1,j) is the skin temperature, and the tslb(i,num_soil_layers,j) level is the
  2942. ! annual mean temperature.
  2943. DO j = jts , MIN(jde-1,jte)
  2944. DO i = its , MIN(ide-1,ite)
  2945. IF ( landmask(i,j) .GT. 0.5 ) THEN
  2946. DO l = 1 , num_soil_layers
  2947. tslb(i,l,j)= ( tsk(i,j) * ( zs(num_soil_layers) - zs(l) ) + &
  2948. tmn(i,j) * ( zs( l) - zs(1) ) ) / &
  2949. ( zs(num_soil_layers) - zs(1) )
  2950. END DO
  2951. ELSE
  2952. IF ( ( real_data_init_type .EQ. 1 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
  2953. DO l = 1 , num_soil_layers
  2954. tslb(i,l,j)= sst(i,j)
  2955. END DO
  2956. ELSE
  2957. DO l = 1 , num_soil_layers
  2958. tslb(i,l,j)= tsk(i,j)
  2959. END DO
  2960. END IF
  2961. END IF
  2962. END DO
  2963. END DO
  2964. END SUBROUTINE init_soil_1_real
  2965. SUBROUTINE init_soil_1_ideal(tsk,tmn,tslb,xland, &
  2966. ivgtyp,ZS,DZS,num_soil_layers, &
  2967. ids,ide, jds,jde, kds,kde, &
  2968. ims,ime, jms,jme, kms,kme, &
  2969. its,ite, jts,jte, kts,kte )
  2970. IMPLICIT NONE
  2971. INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
  2972. ims,ime, jms,jme, kms,kme, &
  2973. its,ite, jts,jte, kts,kte
  2974. INTEGER, INTENT(IN ) :: num_soil_layers
  2975. REAL, DIMENSION( ims:ime , 1 , jms:jme ), INTENT(OUT) :: tslb
  2976. REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: xland
  2977. INTEGER, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: ivgtyp
  2978. REAL, DIMENSION(1:), INTENT(IN) :: dzs,zs
  2979. REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN) :: tsk, tmn
  2980. ! Lcal variables.
  2981. INTEGER :: l,j,i,itf,jtf
  2982. itf=MIN(ite,ide-1)
  2983. jtf=MIN(jte,jde-1)
  2984. IF (num_soil_layers.NE.1)THEN
  2985. DO j=jts,jtf
  2986. DO l=1,num_soil_layers
  2987. DO i=its,itf
  2988. tslb(i,l,j)=( tsk(i,j)*(zs(num_soil_layers)-zs(l)) + tmn(i,j)*(zs(l)-zs(1)) ) / &
  2989. ( zs(num_soil_layers)-zs(1) )
  2990. ENDDO
  2991. ENDDO
  2992. ENDDO
  2993. ENDIF
  2994. DO j=jts,jtf
  2995. DO i=its,itf
  2996. xland(i,j) = 2
  2997. ivgtyp(i,j) = 7
  2998. ENDDO
  2999. ENDDO
  3000. END SUBROUTINE init_soil_1_ideal
  3001. SUBROUTINE init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
  3002. st_input , sm_input , sw_input , landmask , sst , &
  3003. zs , dzs , &
  3004. st_levels_input , sm_levels_input , sw_levels_input , &
  3005. num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  3006. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
  3007. flag_sst , flag_soilt000 , flag_soilm000 , &
  3008. ids , ide , jds , jde , kds , kde , &
  3009. ims , ime , jms , jme , kms , kme , &
  3010. its , ite , jts , jte , kts , kte )
  3011. IMPLICIT NONE
  3012. INTEGER , INTENT(IN) :: num_soil_layers , &
  3013. num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
  3014. num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
  3015. ids , ide , jds , jde , kds , kde , &
  3016. ims , ime , jms , jme , kms , kme , &
  3017. its , ite , jts , jte , kts , kte
  3018. INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
  3019. INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
  3020. INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
  3021. INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
  3022. REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
  3023. REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
  3024. REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
  3025. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
  3026. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
  3027. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
  3028. REAL , DIMENSION(num_soil_layers) :: zs , dzs
  3029. REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
  3030. REAL , ALLOCATABLE , DIMENSION(:) :: zhave
  3031. INTEGER :: i , j , l , lout , lin , lwant , lhave , num
  3032. REAL :: temp
  3033. LOGICAL :: found_levels
  3034. CHARACTER (LEN=132) :: message
  3035. ! Are there any soil temp and moisture levels - ya know, they are mandatory.
  3036. num = num_st_levels_input * num_sm_levels_input
  3037. IF ( num .GE. 1 ) THEN
  3038. ! Ordered levels that we have data for.
  3039. IF ( flag_soilt000 .eq. 1 ) THEN
  3040. write(message, FMT='(A)') ' Assume RUC LSM 6-level input'
  3041. CALL wrf_message ( message )
  3042. ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) ) )
  3043. ELSE
  3044. write(message, FMT='(A)') ' Assume Noah LSM input'
  3045. CALL wrf_message ( message )
  3046. ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input,num_sw_levels_input) +2) )
  3047. END IF
  3048. ! Sort the levels for temperature.
  3049. !print*,'num_st_levels_input, st_levels_input', num_st_levels_input, st_levels_input
  3050. !print*,'num_sm_levels_input,num_sw_levels_input',num_sm_levels_input,num_sw_levels_input
  3051. outert : DO lout = 1 , num_st_levels_input-1
  3052. innert : DO lin = lout+1 , num_st_levels_input
  3053. IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
  3054. temp = st_levels_input(lout)
  3055. st_levels_input(lout) = st_levels_input(lin)
  3056. st_levels_input(lin) = NINT(temp)
  3057. DO j = jts , MIN(jde-1,jte)
  3058. DO i = its , MIN(ide-1,ite)
  3059. temp = st_input(i,lout+1,j)
  3060. st_input(i,lout+1,j) = st_input(i,lin+1,j)
  3061. st_input(i,lin+1,j) = temp
  3062. END DO
  3063. END DO
  3064. END IF
  3065. END DO innert
  3066. END DO outert
  3067. !tgs add IF
  3068. IF ( flag_soilt000 .NE. 1 ) THEN
  3069. DO j = jts , MIN(jde-1,jte)
  3070. DO i = its , MIN(ide-1,ite)
  3071. st_input(i,1,j) = tsk(i,j)
  3072. st_input(i,num_st_levels_input+2,j) = tmn(i,j)
  3073. END DO
  3074. END DO
  3075. ENDIF
  3076. #if ( NMM_CORE == 1 )
  3077. !new
  3078. ! write(0,*) 'st_input(1) in init_soil_2_real'
  3079. ! DO J=MIN(jde-1,jte),JTS, -MIN(jde-1,jte)/15
  3080. ! write(0,616) (st_input(I,1,J),I=its , MIN(ide-1,ite),MIN(ide-1,ite)/10)
  3081. ! ENDDO
  3082. ! write(0,*) 'st_input(2) in init_soil_2_real'
  3083. ! DO J=MIN(jde-1,jte),JTS, -MIN(jde-1,jte)/15
  3084. ! write(0,616) (st_input(I,2,J),I=its , MIN(ide-1,ite),MIN(ide-1,ite)/10)
  3085. ! ENDDO
  3086. 616 format(20(f4.0,1x))
  3087. #endif
  3088. ! Sort the levels for moisture.
  3089. outerm: DO lout = 1 , num_sm_levels_input-1
  3090. innerm : DO lin = lout+1 , num_sm_levels_input
  3091. IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
  3092. temp = sm_levels_input(lout)
  3093. sm_levels_input(lout) = sm_levels_input(lin)
  3094. sm_levels_input(lin) = NINT(temp)
  3095. DO j = jts , MIN(jde-1,jte)
  3096. DO i = its , MIN(ide-1,ite)
  3097. temp = sm_input(i,lout+1,j)
  3098. sm_input(i,lout+1,j) = sm_input(i,lin+1,j)
  3099. sm_input(i,lin+1,j) = temp
  3100. END DO
  3101. END DO
  3102. END IF
  3103. END DO innerm
  3104. END DO outerm
  3105. !tgs add IF
  3106. IF ( flag_soilm000 .NE. 1 ) THEN
  3107. DO j = jts , MIN(jde-1,jte)
  3108. DO i = its , MIN(ide-1,ite)
  3109. sm_input(i,1,j) = sm_input(i,2,j)
  3110. sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
  3111. END DO
  3112. END DO
  3113. ENDIF
  3114. ! Sort the levels for liquid moisture.
  3115. outerw: DO lout = 1 , num_sw_levels_input-1
  3116. innerw : DO lin = lout+1 , num_sw_levels_input
  3117. IF ( sw_levels_input(lout) .GT. sw_levels_input(lin) ) THEN
  3118. temp = sw_levels_input(lout)
  3119. sw_levels_input(lout) = sw_levels_input(lin)
  3120. sw_levels_input(lin) = NINT(temp)
  3121. DO j = jts , MIN(jde-1,jte)
  3122. DO i = its , MIN(ide-1,ite)
  3123. temp = sw_input(i,lout+1,j)
  3124. sw_input(i,lout+1,j) = sw_input(i,lin+1,j)
  3125. sw_input(i,lin+1,j) = temp
  3126. END DO
  3127. END DO
  3128. END IF
  3129. END DO innerw
  3130. END DO outerw
  3131. IF ( num_sw_levels_input .GT. 1 ) THEN
  3132. DO j = jts , MIN(jde-1,jte)
  3133. DO i = its , MIN(ide-1,ite)
  3134. sw_input(i,1,j) = sw_input(i,2,j)
  3135. sw_input(i,num_sw_levels_input+2,j) = sw_input(i,num_sw_levels_input+1,j)
  3136. END DO
  3137. END DO
  3138. END IF
  3139. found_levels = .TRUE.
  3140. ELSE IF ( ( num .LE. 0 ) .AND. ( start_date .NE. current_date ) ) THEN
  3141. found_levels = .FALSE.
  3142. ELSE
  3143. CALL wrf_error_fatal ( &
  3144. 'No input soil level data (temperature, moisture or liquid, or all are missing). Required for LSM.' )
  3145. END IF
  3146. ! Is it OK to continue?
  3147. IF ( found_levels ) THEN
  3148. !tgs add IF
  3149. IF ( flag_soilt000 .NE.1) THEN
  3150. !tgs initialize from Noah data
  3151. ! Here are the levels that we have from the input for temperature. The input levels plus
  3152. ! two more: the skin temperature at 0 cm, and the annual mean temperature at 300 cm.
  3153. zhave(1) = 0.
  3154. DO l = 1 , num_st_levels_input
  3155. zhave(l+1) = st_levels_input(l) / 100.
  3156. END DO
  3157. zhave(num_st_levels_input+2) = 300. / 100.
  3158. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  3159. z_wantt : DO lwant = 1 , num_soil_layers
  3160. z_havet : DO lhave = 1 , num_st_levels_input +2 -1
  3161. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  3162. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  3163. DO j = jts , MIN(jde-1,jte)
  3164. DO i = its , MIN(ide-1,ite)
  3165. tslb(i,lwant,j)= ( st_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
  3166. st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  3167. ( zhave(lhave+1) - zhave(lhave) )
  3168. END DO
  3169. END DO
  3170. EXIT z_havet
  3171. END IF
  3172. END DO z_havet
  3173. END DO z_wantt
  3174. ELSE
  3175. !tgs initialize from RUCLSM data
  3176. DO l = 1 , num_st_levels_input
  3177. zhave(l) = st_levels_input(l) / 100.
  3178. END DO
  3179. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  3180. z_wantt_2 : DO lwant = 1 , num_soil_layers
  3181. z_havet_2 : DO lhave = 1 , num_st_levels_input -1
  3182. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  3183. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  3184. DO j = jts , MIN(jde-1,jte)
  3185. DO i = its , MIN(ide-1,ite)
  3186. tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
  3187. st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  3188. ( zhave(lhave+1) - zhave(lhave) )
  3189. END DO
  3190. END DO
  3191. EXIT z_havet_2
  3192. END IF
  3193. END DO z_havet_2
  3194. END DO z_wantt_2
  3195. ENDIF
  3196. IF ( flag_soilm000 .NE. 1 ) THEN
  3197. !tgs initialize from Noah
  3198. ! Here are the levels that we have from the input for moisture. The input levels plus
  3199. ! two more: a value at 0 cm and one at 300 cm. The 0 cm value is taken to be identical
  3200. ! to the most shallow layer's value. Similarly, the 300 cm value is taken to be the same
  3201. ! as the most deep layer's value.
  3202. zhave(1) = 0.
  3203. DO l = 1 , num_sm_levels_input
  3204. zhave(l+1) = sm_levels_input(l) / 100.
  3205. END DO
  3206. zhave(num_sm_levels_input+2) = 300. / 100.
  3207. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  3208. z_wantm : DO lwant = 1 , num_soil_layers
  3209. z_havem : DO lhave = 1 , num_sm_levels_input +2 -1
  3210. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  3211. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  3212. DO j = jts , MIN(jde-1,jte)
  3213. DO i = its , MIN(ide-1,ite)
  3214. smois(i,lwant,j)= ( sm_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
  3215. sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  3216. ( zhave(lhave+1) - zhave(lhave) )
  3217. END DO
  3218. END DO
  3219. EXIT z_havem
  3220. END IF
  3221. END DO z_havem
  3222. END DO z_wantm
  3223. ELSE
  3224. !tgs initialize from RUCLSM data
  3225. DO l = 1 , num_sm_levels_input
  3226. zhave(l) = sm_levels_input(l) / 100.
  3227. END DO
  3228. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  3229. z_wantm_2 : DO lwant = 1 , num_soil_layers
  3230. z_havem_2 : DO lhave = 1 , num_sm_levels_input -1
  3231. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  3232. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  3233. DO j = jts , MIN(jde-1,jte)
  3234. DO i = its , MIN(ide-1,ite)
  3235. smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
  3236. sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  3237. ( zhave(lhave+1) - zhave(lhave) )
  3238. END DO
  3239. END DO
  3240. EXIT z_havem_2
  3241. END IF
  3242. END DO z_havem_2
  3243. END DO z_wantm_2
  3244. ENDIF
  3245. ! Any liquid soil moisture to worry about?
  3246. IF ( num_sw_levels_input .GT. 1 ) THEN
  3247. zhave(1) = 0.
  3248. DO l = 1 , num_sw_levels_input
  3249. zhave(l+1) = sw_levels_input(l) / 100.
  3250. END DO
  3251. zhave(num_sw_levels_input+2) = 300. / 100.
  3252. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  3253. z_wantw : DO lwant = 1 , num_soil_layers
  3254. z_havew : DO lhave = 1 , num_sw_levels_input +2 -1
  3255. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  3256. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  3257. DO j = jts , MIN(jde-1,jte)
  3258. DO i = its , MIN(ide-1,ite)
  3259. sh2o(i,lwant,j)= ( sw_input(i,lhave ,j) * ( zhave(lhave+1) - zs (lwant) ) + &
  3260. sw_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  3261. ( zhave(lhave+1) - zhave(lhave) )
  3262. END DO
  3263. END DO
  3264. EXIT z_havew
  3265. END IF
  3266. END DO z_havew
  3267. END DO z_wantw
  3268. END IF
  3269. ! Over water, put in reasonable values for soil temperature and moisture. These won't be
  3270. ! used, but they will make a more continuous plot.
  3271. IF ( flag_sst .EQ. 1 ) THEN
  3272. DO j = jts , MIN(jde-1,jte)
  3273. DO i = its , MIN(ide-1,ite)
  3274. IF ( landmask(i,j) .LT. 0.5 ) THEN
  3275. DO l = 1 , num_soil_layers
  3276. tslb(i,l,j)= sst(i,j)
  3277. !tgs add line for tsk
  3278. tsk(i,j) = sst(i,j)
  3279. smois(i,l,j)= 1.0
  3280. sh2o (i,l,j)= 1.0
  3281. END DO
  3282. END IF
  3283. END DO
  3284. END DO
  3285. ELSE
  3286. DO j = jts , MIN(jde-1,jte)
  3287. DO i = its , MIN(ide-1,ite)
  3288. IF ( landmask(i,j) .LT. 0.5 ) THEN
  3289. DO l = 1 , num_soil_layers
  3290. tslb(i,l,j)= tsk(i,j)
  3291. smois(i,l,j)= 1.0
  3292. sh2o (i,l,j)= 1.0
  3293. END DO
  3294. END IF
  3295. END DO
  3296. END DO
  3297. END IF
  3298. DEALLOCATE (zhave)
  3299. END IF
  3300. END SUBROUTINE init_soil_2_real
  3301. SUBROUTINE init_soil_2_ideal ( xland,xice,vegfra,snow,canwat, &
  3302. ivgtyp,isltyp,tslb,smois,tmn, &
  3303. num_soil_layers, &
  3304. ids,ide, jds,jde, kds,kde, &
  3305. ims,ime, jms,jme, kms,kme, &
  3306. its,ite, jts,jte, kts,kte )
  3307. IMPLICIT NONE
  3308. INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde, &
  3309. ims,ime, jms,jme, kms,kme, &
  3310. its,ite, jts,jte, kts,kte
  3311. INTEGER, INTENT(IN) ::num_soil_layers
  3312. REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(OUT) :: smois, tslb
  3313. REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra, tmn
  3314. INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
  3315. INTEGER :: icm,jcm,itf,jtf
  3316. INTEGER :: i,j,l
  3317. itf=min0(ite,ide-1)
  3318. jtf=min0(jte,jde-1)
  3319. icm = ide/2
  3320. jcm = jde/2
  3321. DO j=jts,jtf
  3322. DO l=1,num_soil_layers
  3323. DO i=its,itf
  3324. smois(i,1,j)=0.10
  3325. smois(i,2,j)=0.10
  3326. smois(i,3,j)=0.10
  3327. smois(i,4,j)=0.10
  3328. tslb(i,1,j)=295.
  3329. tslb(i,2,j)=297.
  3330. tslb(i,3,j)=293.
  3331. tslb(i,4,j)=293.
  3332. ENDDO
  3333. ENDDO
  3334. ENDDO
  3335. DO j=jts,jtf
  3336. DO i=its,itf
  3337. xland(i,j) = 2
  3338. tmn(i,j) = 294.
  3339. xice(i,j) = 0.
  3340. vegfra(i,j) = 0.
  3341. snow(i,j) = 0.
  3342. canwat(i,j) = 0.
  3343. ivgtyp(i,j) = 7
  3344. isltyp(i,j) = 8
  3345. ENDDO
  3346. ENDDO
  3347. END SUBROUTINE init_soil_2_ideal
  3348. SUBROUTINE init_soil_3_real ( tsk , tmn , smois , tslb , &
  3349. st_input , sm_input , landmask, sst, &
  3350. zs , dzs , &
  3351. st_levels_input , sm_levels_input , &
  3352. num_soil_layers , num_st_levels_input , num_sm_levels_input , &
  3353. num_st_levels_alloc , num_sm_levels_alloc , &
  3354. flag_sst , flag_soilt000 , flag_soilm000 , &
  3355. ids , ide , jds , jde , kds , kde , &
  3356. ims , ime , jms , jme , kms , kme , &
  3357. its , ite , jts , jte , kts , kte )
  3358. IMPLICIT NONE
  3359. INTEGER , INTENT(IN) :: num_soil_layers , &
  3360. num_st_levels_input , num_sm_levels_input , &
  3361. num_st_levels_alloc , num_sm_levels_alloc , &
  3362. ids , ide , jds , jde , kds , kde , &
  3363. ims , ime , jms , jme , kms , kme , &
  3364. its , ite , jts , jte , kts , kte
  3365. INTEGER , INTENT(IN) :: flag_sst, flag_soilt000, flag_soilm000
  3366. INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
  3367. INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
  3368. REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
  3369. REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
  3370. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
  3371. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tmn
  3372. REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk
  3373. REAL , DIMENSION(num_soil_layers) :: zs , dzs
  3374. REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois
  3375. REAL , ALLOCATABLE , DIMENSION(:) :: zhave
  3376. INTEGER :: i , j , l , lout , lin , lwant , lhave
  3377. REAL :: temp
  3378. ! Allocate the soil layer array used for interpolating.
  3379. IF ( ( num_st_levels_input .LE. 0 ) .OR. &
  3380. ( num_sm_levels_input .LE. 0 ) ) THEN
  3381. PRINT '(A)','No input soil level data (either temperature or moisture, or both are missing). Required for RUC LSM.'
  3382. CALL wrf_error_fatal ( 'no soil data' )
  3383. ELSE
  3384. IF ( flag_soilt000 .eq. 1 ) THEN
  3385. PRINT '(A)',' Assume RUC LSM 6-level input'
  3386. ALLOCATE ( zhave( MAX(num_st_levels_input,num_sm_levels_input) ) )
  3387. ELSE
  3388. PRINT '(A)',' Assume non-RUC LSM input'
  3389. ALLOCATE ( zhave( MAX(num_st_levels_input,num_soil_layers) ) )
  3390. END IF
  3391. END IF
  3392. ! Sort the levels for temperature.
  3393. outert : DO lout = 1 , num_st_levels_input-1
  3394. innert : DO lin = lout+1 , num_st_levels_input
  3395. IF ( st_levels_input(lout) .GT. st_levels_input(lin) ) THEN
  3396. temp = st_levels_input(lout)
  3397. st_levels_input(lout) = st_levels_input(lin)
  3398. st_levels_input(lin) = NINT(temp)
  3399. DO j = jts , MIN(jde-1,jte)
  3400. DO i = its , MIN(ide-1,ite)
  3401. temp = st_input(i,lout,j)
  3402. st_input(i,lout,j) = st_input(i,lin,j)
  3403. st_input(i,lin,j) = temp
  3404. END DO
  3405. END DO
  3406. END IF
  3407. END DO innert
  3408. END DO outert
  3409. IF ( flag_soilt000 .NE. 1 ) THEN
  3410. DO j = jts , MIN(jde-1,jte)
  3411. DO i = its , MIN(ide-1,ite)
  3412. st_input(i,1,j) = tsk(i,j)
  3413. st_input(i,num_st_levels_input+2,j) = tmn(i,j)
  3414. END DO
  3415. END DO
  3416. END IF
  3417. ! Sort the levels for moisture.
  3418. outerm: DO lout = 1 , num_sm_levels_input-1
  3419. innerm : DO lin = lout+1 , num_sm_levels_input
  3420. IF ( sm_levels_input(lout) .GT. sm_levels_input(lin) ) THEN
  3421. temp = sm_levels_input(lout)
  3422. sm_levels_input(lout) = sm_levels_input(lin)
  3423. sm_levels_input(lin) = NINT(temp)
  3424. DO j = jts , MIN(jde-1,jte)
  3425. DO i = its , MIN(ide-1,ite)
  3426. temp = sm_input(i,lout,j)
  3427. sm_input(i,lout,j) = sm_input(i,lin,j)
  3428. sm_input(i,lin,j) = temp
  3429. END DO
  3430. END DO
  3431. END IF
  3432. END DO innerm
  3433. END DO outerm
  3434. IF ( flag_soilm000 .NE. 1 ) THEN
  3435. DO j = jts , MIN(jde-1,jte)
  3436. DO i = its , MIN(ide-1,ite)
  3437. sm_input(i,1,j) = sm_input(i,2,j)
  3438. sm_input(i,num_sm_levels_input+2,j) = sm_input(i,num_sm_levels_input+1,j)
  3439. END DO
  3440. END DO
  3441. END IF
  3442. ! Here are the levels that we have from the input for temperature.
  3443. IF ( flag_soilt000 .EQ. 1 ) THEN
  3444. DO l = 1 , num_st_levels_input
  3445. zhave(l) = st_levels_input(l) / 100.
  3446. END DO
  3447. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  3448. z_wantt : DO lwant = 1 , num_soil_layers
  3449. z_havet : DO lhave = 1 , num_st_levels_input -1
  3450. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  3451. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  3452. DO j = jts , MIN(jde-1,jte)
  3453. DO i = its , MIN(ide-1,ite)
  3454. tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
  3455. st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  3456. ( zhave(lhave+1) - zhave(lhave) )
  3457. END DO
  3458. END DO
  3459. EXIT z_havet
  3460. END IF
  3461. END DO z_havet
  3462. END DO z_wantt
  3463. ELSE
  3464. zhave(1) = 0.
  3465. DO l = 1 , num_st_levels_input
  3466. zhave(l+1) = st_levels_input(l) / 100.
  3467. END DO
  3468. zhave(num_st_levels_input+2) = 300. / 100.
  3469. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  3470. z_wantt_2 : DO lwant = 1 , num_soil_layers
  3471. z_havet_2 : DO lhave = 1 , num_st_levels_input +2
  3472. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  3473. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  3474. DO j = jts , MIN(jde-1,jte)
  3475. DO i = its , MIN(ide-1,ite)
  3476. tslb(i,lwant,j)= ( st_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
  3477. st_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  3478. ( zhave(lhave+1) - zhave(lhave) )
  3479. END DO
  3480. END DO
  3481. EXIT z_havet_2
  3482. END IF
  3483. END DO z_havet_2
  3484. END DO z_wantt_2
  3485. END IF
  3486. ! Here are the levels that we have from the input for moisture.
  3487. IF ( flag_soilm000 .EQ. 1 ) THEN
  3488. DO l = 1 , num_sm_levels_input
  3489. zhave(l) = sm_levels_input(l) / 100.
  3490. END DO
  3491. ! Interpolate between the layers we have (zhave) and those that we want (zs).
  3492. z_wantm : DO lwant = 1 , num_soil_layers
  3493. z_havem : DO lhave = 1 , num_sm_levels_input -1
  3494. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  3495. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  3496. DO j = jts , MIN(jde-1,jte)
  3497. DO i = its , MIN(ide-1,ite)
  3498. smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
  3499. sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  3500. ( zhave(lhave+1) - zhave(lhave) )
  3501. END DO
  3502. END DO
  3503. EXIT z_havem
  3504. END IF
  3505. END DO z_havem
  3506. END DO z_wantm
  3507. ELSE
  3508. zhave(1) = 0.
  3509. DO l = 1 , num_sm_levels_input
  3510. zhave(l+1) = sm_levels_input(l) / 100.
  3511. END DO
  3512. zhave(num_sm_levels_input+2) = 300. / 100.
  3513. z_wantm_2 : DO lwant = 1 , num_soil_layers
  3514. z_havem_2 : DO lhave = 1 , num_sm_levels_input +2
  3515. IF ( ( zs(lwant) .GE. zhave(lhave ) ) .AND. &
  3516. ( zs(lwant) .LE. zhave(lhave+1) ) ) THEN
  3517. DO j = jts , MIN(jde-1,jte)
  3518. DO i = its , MIN(ide-1,ite)
  3519. smois(i,lwant,j)= ( sm_input(i,lhave,j ) * ( zhave(lhave+1) - zs (lwant) ) + &
  3520. sm_input(i,lhave+1,j) * ( zs (lwant ) - zhave(lhave) ) ) / &
  3521. ( zhave(lhave+1) - zhave(lhave) )
  3522. END DO
  3523. END DO
  3524. EXIT z_havem_2
  3525. END IF
  3526. END DO z_havem_2
  3527. END DO z_wantm_2
  3528. END IF
  3529. ! Over water, put in reasonable values for soil temperature and moisture. These won't be
  3530. ! used, but they will make a more continuous plot.
  3531. IF ( flag_sst .EQ. 1 ) THEN
  3532. DO j = jts , MIN(jde-1,jte)
  3533. DO i = its , MIN(ide-1,ite)
  3534. IF ( landmask(i,j) .LT. 0.5 ) THEN
  3535. DO l = 1 , num_soil_layers
  3536. tslb(i,l,j) = sst(i,j)
  3537. tsk(i,j) = sst(i,j)
  3538. smois(i,l,j)= 1.0
  3539. END DO
  3540. END IF
  3541. END DO
  3542. END DO
  3543. ELSE
  3544. DO j = jts , MIN(jde-1,jte)
  3545. DO i = its , MIN(ide-1,ite)
  3546. IF ( landmask(i,j) .LT. 0.5 ) THEN
  3547. DO l = 1 , num_soil_layers
  3548. tslb(i,l,j)= tsk(i,j)
  3549. smois(i,l,j)= 1.0
  3550. END DO
  3551. END IF
  3552. END DO
  3553. END DO
  3554. END IF
  3555. DEALLOCATE (zhave)
  3556. END SUBROUTINE init_soil_3_real
  3557. END MODULE module_soil_pre
  3558. #endif