PageRenderTime 57ms CodeModel.GetById 13ms 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

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

    

  1. #if ( ! NMM_CORE == 1 )
    2. 

  2. 

  3. MODULE module_soil_pre
    4. 

  4. 

  5.    USE module_date_time
    6. 

  6.    USE module_state_description
    7. 

  7. 

  8.    CHARACTER (LEN=3) :: num_cat_count
    9. 

  9.    INTEGER , PARAMETER , DIMENSION(0:300) :: ints = &
    10. 

  10.    (/    0,    1,    2,    3,    4,    5,    6,    7,    8,    9, &
    11. 

  11.         10,   11,   12,   13,   14,   15,   16,   17,   18,   19, &
    12. 

  12.         20,   21,   22,   23,   24,   25,   26,   27,   28,   29, &
    13. 

  13.         30,   31,   32,   33,   34,   35,   36,   37,   38,   39, &
    14. 

  14.         40,   41,   42,   43,   44,   45,   46,   47,   48,   49, &
    15. 

  15.         50,   51,   52,   53,   54,   55,   56,   57,   58,   59, &
    16. 

  16.         60,   61,   62,   63,   64,   65,   66,   67,   68,   69, &
    17. 

  17.         70,   71,   72,   73,   74,   75,   76,   77,   78,   79, &
    18. 

  18.         80,   81,   82,   83,   84,   85,   86,   87,   88,   89, &
    19. 

  19.         90,   91,   92,   93,   94,   95,   96,   97,   98,   99, &
    20. 

  20.        100,  101,  102,  103,  104,  105,  106,  107,  108,  109, &
    21. 

  21.        110,  111,  112,  113,  114,  115,  116,  117,  118,  119, &
    22. 

  22.        120,  121,  122,  123,  124,  125,  126,  127,  128,  129, &
    23. 

  23.        130,  131,  132,  133,  134,  135,  136,  137,  138,  139, &
    24. 

  24.        140,  141,  142,  143,  144,  145,  146,  147,  148,  149, &
    25. 

  25.        150,  151,  152,  153,  154,  155,  156,  157,  158,  159, &
    26. 

  26.        160,  161,  162,  163,  164,  165,  166,  167,  168,  169, &
    27. 

  27.        170,  171,  172,  173,  174,  175,  176,  177,  178,  179, &
    28. 

  28.        180,  181,  182,  183,  184,  185,  186,  187,  188,  189, &
    29. 

  29.        190,  191,  192,  193,  194,  195,  196,  197,  198,  199, &
    30. 

  30.        200,  201,  202,  203,  204,  205,  206,  207,  208,  209, &
    31. 

  31.        210,  211,  212,  213,  214,  215,  216,  217,  218,  219, &
    32. 

  32.        220,  221,  222,  223,  224,  225,  226,  227,  228,  229, &
    33. 

  33.        230,  231,  232,  233,  234,  235,  236,  237,  238,  239, &
    34. 

  34.        240,  241,  242,  243,  244,  245,  246,  247,  248,  249, &
    35. 

  35.        250,  251,  252,  253,  254,  255,  256,  257,  258,  259, &
    36. 

  36.        260,  261,  262,  263,  264,  265,  266,  267,  268,  269, &
    37. 

  37.        270,  271,  272,  273,  274,  275,  276,  277,  278,  279, &
    38. 

  38.        280,  281,  282,  283,  284,  285,  286,  287,  288,  289, &
    39. 

  39.        290,  291,  292,  293,  294,  295,  296,  297,  298,  299, 300 /)
    40. 

  40. 

  41.    !  Excluded middle processing
    42. 

  42. 

  43.    LOGICAL , SAVE :: hold_ups
    44. 

  44.    INTEGER , SAVE :: em_width
    45. 

  45. 

  46.    LOGICAL , EXTERNAL :: skip_middle_points_t
    47. 

  47. 

  48. CONTAINS
    49. 

  49. 

  50.    SUBROUTINE adjust_for_seaice_pre ( xice , landmask , tsk , ivgtyp , vegcat , lu_index , &
    51. 

  51.                                       xland , landusef , isltyp , soilcat , soilctop , &
    52. 

  52.                                       soilcbot , tmn , &
    53. 

  53.                                       seaice_threshold , &
    54. 

  54.                                       fractional_seaice, &
    55. 

  55.                                       num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
    56. 

  56.                                       iswater , isice , &
    57. 

  57.                                       scheme , &
    58. 

  58.                                       ids , ide , jds , jde , kds , kde , &
    59. 

  59.                                       ims , ime , jms , jme , kms , kme , &
    60. 

  60.                                       its , ite , jts , jte , kts , kte )
    61. 

  61. 

  62.       IMPLICIT NONE
    63. 

  63. 

  64.       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
    65. 

  65.                               ims , ime , jms , jme , kms , kme , &
    66. 

  66.                               its , ite , jts , jte , kts , kte , &
    67. 

  67.                               iswater , isice
    68. 

  68.       INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
    69. 

  69. 

  70.       REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
    71. 

  71.       REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
    72. 

  72.       REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
    73. 

  73.       INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
    74. 

  74.       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
    75. 

  75.                                                            vegcat, xland , soilcat , tmn
    76. 

  76.       REAL , INTENT(IN) :: seaice_threshold
    77. 

  77. 

  78.       INTEGER :: i , j , num_seaice_changes , loop
    79. 

  79.       CHARACTER (LEN=132) :: message
    80. 

  80. 

  81.       INTEGER, INTENT(IN) :: fractional_seaice
    82. 

  82.       REAL :: XICE_THRESHOLD
    83. 

  83. 

  84.       IF ( FRACTIONAL_SEAICE == 0 ) THEN
    85. 

  85.          xice_threshold = 0.5
    86. 

  86.       ELSEIF ( FRACTIONAL_SEAICE == 1 ) THEN
    87. 

  87.          xice_threshold = 0.02
    88. 

  88.       ENDIF
    89. 

  89. 

  90.       num_seaice_changes = 0
    91. 

  91.       fix_seaice : SELECT CASE ( scheme )
    92. 

  92. 

  93.          CASE ( SLABSCHEME )
    94. 

  94.             DO j = jts , MIN(jde-1,jte)
    95. 

  95.                DO i = its , MIN(ide-1,ite)
    96. 

  96.                   IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    97. 

  97.                   IF ( xice(i,j) .GT. 200.0 ) THEN
    98. 

  98.                      xice(i,j) = 0.
    99. 

  99.                      num_seaice_changes = num_seaice_changes + 1
    100. 

  100.                   END IF
    101. 

  101.                END DO
    102. 

  102.             END DO
    103. 

  103.             IF ( num_seaice_changes .GT. 0 ) THEN
    104. 

  104.                WRITE ( message , FMT='(A,I6)' ) &
    105. 

  105.                'Total pre number of sea ice locations removed (due to FLAG values) = ', &
    106. 

  106.                num_seaice_changes
    107. 

  107.                CALL wrf_debug ( 0 , message )
    108. 

  108.             END IF
    109. 

  109.             num_seaice_changes = 0
    110. 

  110.             DO j = jts , MIN(jde-1,jte)
    111. 

  111.                DO i = its , MIN(ide-1,ite)
    112. 

  112.                   IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    113. 

  113.                   IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
    114. 

  114.                        ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
    115. 

  115.                      IF ( FRACTIONAL_SEAICE == 0 ) THEN
    116. 

  116.                         xice(i,j) = 1.0
    117. 

  117.                      ENDIF
    118. 

  118.                      num_seaice_changes = num_seaice_changes + 1
    119. 

  119.                      if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
    120. 

  120.                      vegcat(i,j)=isice
    121. 

  121.                      ivgtyp(i,j)=isice
    122. 

  122.                      lu_index(i,j)=isice
    123. 

  123.                      landmask(i,j)=1.
    124. 

  124.                      xland(i,j)=1.
    125. 

  125.                      DO loop=1,num_veg_cat
    126. 

  126.                         landusef(i,loop,j)=0.
    127. 

  127.                      END DO
    128. 

  128.                      landusef(i,ivgtyp(i,j),j)=1.
    129. 

  129. 

  130.                      isltyp(i,j) = 16
    131. 

  131.                      soilcat(i,j)=isltyp(i,j)
    132. 

  132.                      DO loop=1,num_soil_top_cat
    133. 

  133.                         soilctop(i,loop,j)=0
    134. 

  134.                      END DO
    135. 

  135.                      DO loop=1,num_soil_bot_cat
    136. 

  136.                         soilcbot(i,loop,j)=0
    137. 

  137.                      END DO
    138. 

  138.                      soilctop(i,isltyp(i,j),j)=1.
    139. 

  139.                      soilcbot(i,isltyp(i,j),j)=1.
    140. 

  140.                   ELSE
    141. 

  141.                      xice(i,j) = 0.0
    142. 

  142.                   END IF
    143. 

  143.                END DO
    144. 

  144.             END DO
    145. 

  145.             IF ( num_seaice_changes .GT. 0 ) THEN
    146. 

  146.                WRITE ( message , FMT='(A,I6)' ) &
    147. 

  147.                'Total pre number of sea ice location changes (water to land) = ', num_seaice_changes
    148. 

  148.                CALL wrf_debug ( 0 , message )
    149. 

  149.             END IF
    150. 

  150. 

  151.          CASE ( LSMSCHEME , NOAHMPSCHEME , RUCLSMSCHEME , SSIBSCHEME)   !mchen add for ssib
    152. 

  152.             num_seaice_changes = 0
    153. 

  153.             DO j = jts , MIN(jde-1,jte)
    154. 

  154.                DO i = its , MIN(ide-1,ite)
    155. 

  155.                   IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    156. 

  156.                   IF ( landmask(i,j) .GT. 0.5 ) THEN
    157. 

  157.                      if (xice(i,j).gt.0) num_seaice_changes = num_seaice_changes + 1
    158. 

  158.                      xice(i,j) = 0.
    159. 

  159.                   END IF
    160. 

  160.                END DO
    161. 

  161.             END DO
    162. 

  162.             IF ( num_seaice_changes .GT. 0 ) THEN
    163. 

  163.                WRITE ( message , FMT='(A,I6)' ) &
    164. 

  164.                'Total pre number of land location changes (seaice set to zero) = ', num_seaice_changes
    165. 

  165.                CALL wrf_debug ( 0 , message )
    166. 

  166.             END IF
    167. 

  167. 

  168.       END SELECT fix_seaice
    169. 

  169. 

  170.    END SUBROUTINE adjust_for_seaice_pre
    171. 

  171. 

  172.    SUBROUTINE adjust_for_seaice_post ( xice , landmask , tsk_old , tsk , ivgtyp , vegcat , lu_index , &
    173. 

  173.                                       xland , landusef , isltyp , soilcat , soilctop , &
    174. 

  174.                                       soilcbot , tmn , vegfra , &
    175. 

  175.                                       tslb , smois , sh2o , &
    176. 

  176.                                       seaice_threshold , &
    177. 

  177.                                       sst , flag_sst , &
    178. 

  178.                                       fractional_seaice, &
    179. 

  179.                                       num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
    180. 

  180.                                       num_soil_layers , &
    181. 

  181.                                       iswater , isice , &
    182. 

  182.                                       scheme , &
    183. 

  183.                                       ids , ide , jds , jde , kds , kde , &
    184. 

  184.                                       ims , ime , jms , jme , kms , kme , &
    185. 

  185.                                       its , ite , jts , jte , kts , kte )
    186. 

  186. 

  187.       IMPLICIT NONE
    188. 

  188. 

  189.       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
    190. 

  190.                               ims , ime , jms , jme , kms , kme , &
    191. 

  191.                               its , ite , jts , jte , kts , kte , &
    192. 

  192.                               iswater , isice
    193. 

  193.       INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat , scheme
    194. 

  194.       INTEGER , INTENT(IN) :: num_soil_layers
    195. 

  195. 

  196.       REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landusef
    197. 

  197.       REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(INOUT):: soilctop
    198. 

  198.       REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(INOUT):: soilcbot
    199. 

  199.       REAL , DIMENSION(ims:ime,1:num_soil_layers,jms:jme) , INTENT(INOUT):: tslb , smois , sh2o
    200. 

  200.       REAL , DIMENSION(ims:ime,jms:jme)                   , INTENT(IN):: sst
    201. 

  201.       INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
    202. 

  202.       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask , xice , tsk , lu_index , &
    203. 

  203.                                                            vegcat, xland , soilcat , tmn , &
    204. 

  204.                                                            tsk_old , vegfra
    205. 

  205.       INTEGER , INTENT(IN) :: flag_sst
    206. 

  206.       REAL , INTENT(IN) :: seaice_threshold
    207. 

  207.       REAL :: total_depth , mid_point_depth
    208. 

  208. 

  209.       INTEGER :: i , j , num_seaice_changes , loop
    210. 

  210.       CHARACTER (LEN=132) :: message
    211. 

  211. 

  212. 

  213.       INTEGER, INTENT(IN) :: fractional_seaice
    214. 

  214.       real :: xice_threshold
    215. 

  215. 

  216.       IF ( FRACTIONAL_SEAICE == 0 ) THEN
    217. 

  217.          xice_threshold = 0.5
    218. 

  218.       ELSEIF ( FRACTIONAL_SEAICE == 1 ) THEN
    219. 

  219.          xice_threshold = 0.02
    220. 

  220.       ENDIF
    221. 

  221.       num_seaice_changes = 0
    222. 

  222.       fix_seaice : SELECT CASE ( scheme )
    223. 

  223. 

  224.          CASE ( SLABSCHEME )
    225. 

  225. 

  226.          CASE ( LSMSCHEME , NOAHMPSCHEME , SSIBSCHEME )   !mchen add for ssib
    227. 

  227.             DO j = jts , MIN(jde-1,jte)
    228. 

  228.                DO i = its , MIN(ide-1,ite)
    229. 

  229.                   IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    230. 

  230.                   IF ( xice(i,j) .GT. 200.0 ) THEN
    231. 

  231.                      xice(i,j) = 0.
    232. 

  232.                      num_seaice_changes = num_seaice_changes + 1
    233. 

  233.                   END IF
    234. 

  234.                END DO
    235. 

  235.             END DO
    236. 

  236.             IF ( num_seaice_changes .GT. 0 ) THEN
    237. 

  237.                WRITE ( message , FMT='(A,I6)' ) &
    238. 

  238.                'Total post number of sea ice locations removed (due to FLAG values) = ', &
    239. 

  239.                num_seaice_changes
    240. 

  240.                CALL wrf_debug ( 0 , message )
    241. 

  241.             END IF
    242. 

  242.             num_seaice_changes = 0
    243. 

  243.             DO j = jts , MIN(jde-1,jte)
    244. 

  244.                DO i = its , MIN(ide-1,ite)
    245. 

  245.                   IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    246. 

  246.                   IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
    247. 

  247.                        ( ( tsk_old(i,j) .GT. 170 ) .AND. ( tsk_old(i,j) .LT. 400 ) ) )THEN
    248. 

  248.                      tsk(i,j) = tsk_old(i,j)
    249. 

  249.                   END IF
    250. 

  250.                   IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
    251. 

  251.                        ( ( tsk_old(i,j) .LT. 170 ) .OR. ( tsk_old(i,j) .GT. 400 ) ) )THEN
    252. 

  252.                      print *,'TSK woes in seaice post, i,j=',i,j,'  tsk = ',tsk(i,j), tsk_old(i,j)
    253. 

  253.                      CALL wrf_error_fatal ( 'TSK is unrealistic, problems for seaice post')
    254. 

  254.                   ELSE IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
    255. 

  255.                        ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
    256. 

  256.                      IF ( FRACTIONAL_SEAICE == 0 ) THEN
    257. 

  257.                         xice(i,j) = 1.0
    258. 

  258.                      ENDIF
    259. 

  259.                      num_seaice_changes = num_seaice_changes + 1
    260. 

  260.                      if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
    261. 

  261.                      vegcat(i,j)=isice
    262. 

  262.                      ivgtyp(i,j)=isice
    263. 

  263.                      lu_index(i,j)=isice
    264. 

  264.                      landmask(i,j)=1.
    265. 

  265.                      xland(i,j)=1.
    266. 

  266.                      vegfra(i,j)=0.
    267. 

  267.                      DO loop=1,num_veg_cat
    268. 

  268.                         landusef(i,loop,j)=0.
    269. 

  269.                      END DO
    270. 

  270.                      landusef(i,ivgtyp(i,j),j)=1.
    271. 

  271. 

  272.                      tsk_old(i,j) = tsk(i,j)
    273. 

  273. 

  274.                      isltyp(i,j) = 16
    275. 

  275.                      soilcat(i,j)=isltyp(i,j)
    276. 

  276.                      DO loop=1,num_soil_top_cat
    277. 

  277.                         soilctop(i,loop,j)=0
    278. 

  278.                      END DO
    279. 

  279.                      DO loop=1,num_soil_bot_cat
    280. 

  280.                         soilcbot(i,loop,j)=0
    281. 

  281.                      END DO
    282. 

  282.                      soilctop(i,isltyp(i,j),j)=1.
    283. 

  283.                      soilcbot(i,isltyp(i,j),j)=1.
    284. 

  284. 

  285.                      total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
    286. 

  286.                      DO loop = 1,num_soil_layers
    287. 

  287.                         mid_point_depth=(total_depth/num_soil_layers)/2. + &
    288. 

  288.                                         (loop-1)*(total_depth/num_soil_layers)
    289. 

  289.                         tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
    290. 

  290.                                             mid_point_depth*tmn(i,j) ) / total_depth
    291. 

  291.                      END DO
    292. 

  292. 

  293.                      DO loop=1,num_soil_layers
    294. 

  294.                         smois(i,loop,j) = 1.0
    295. 

  295.                         sh2o(i,loop,j)  = 0.0
    296. 

  296.                      END DO
    297. 

  297.                   ELSE IF ( xice(i,j) .LT. xice_threshold ) THEN
    298. 

  298.                      xice(i,j) = 0.
    299. 

  299.                   END IF
    300. 

  300.                END DO
    301. 

  301.             END DO
    302. 

  302.             IF ( num_seaice_changes .GT. 0 ) THEN
    303. 

  303.                WRITE ( message , FMT='(A,I6)' ) &
    304. 

  304.                'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
    305. 

  305.                CALL wrf_debug ( 0 , message )
    306. 

  306.             END IF
    307. 

  307. 

  308.         CASE ( RUCLSMSCHEME )
    309. 

  309.             DO j = jts , MIN(jde-1,jte)
    310. 

  310.                DO i = its , MIN(ide-1,ite)
    311. 

  311.                   IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    312. 

  312.                   IF ( xice(i,j) .GT. 200.0 ) THEN
    313. 

  313.                      xice(i,j) = 0.
    314. 

  314.                      num_seaice_changes = num_seaice_changes + 1
    315. 

  315.                   END IF
    316. 

  316.                END DO
    317. 

  317.             END DO
    318. 

  318.             IF ( num_seaice_changes .GT. 0 ) THEN
    319. 

  319.                WRITE ( message , FMT='(A,I6)' ) &
    320. 

  320.                'Total post number of sea ice locations removed (due to FLAG values) = ', &
    321. 

  321.                num_seaice_changes
    322. 

  322.                CALL wrf_debug ( 0 , message )
    323. 

  323.             END IF
    324. 

  324.             num_seaice_changes = 0
    325. 

  325.             DO j = jts , MIN(jde-1,jte)
    326. 

  326.                DO i = its , MIN(ide-1,ite)
    327. 

  327.                   IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    328. 

  328.                   IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
    329. 

  329.                        ( ( tsk_old(i,j) .GT. 170 ) .AND. ( tsk_old(i,j) .LT. 400 ) ) )THEN
    330. 

  330.                      tsk(i,j) = tsk_old(i,j)
    331. 

  331.                   END IF
    332. 

  332.                   IF ( ( ( tsk(i,j) .LT. 170 ) .OR. ( tsk(i,j) .GT. 400 ) ) .AND. &
    333. 

  333.                        ( ( tsk_old(i,j) .LT. 170 ) .OR. ( tsk_old(i,j) .GT. 400 ) ) )THEN
    334. 

  334.                      print *,'TSK woes in seaice post, i,j=',i,j,'  tsk = ',tsk(i,j), tsk_old(i,j)
    335. 

  335.                      CALL wrf_error_fatal ( 'TSK is unrealistic, problems for seaice post')
    336. 

  336.                   ELSE IF ( ( xice(i,j) .GE. xice_threshold ) .OR. &
    337. 

  337.                        ( ( landmask(i,j) .LT. 0.5 ) .AND. ( tsk(i,j) .LT. seaice_threshold ) ) ) THEN
    338. 

  338.                        IF ( FRACTIONAL_SEAICE == 0 ) THEN
    339. 

  339.                           xice(i,j) = 1.0
    340. 

  340.                        ELSE
    341. 

  341.                           xice(i,j)=max(0.25,xice(i,j))
    342. 

  342.                        ENDIF
    343. 

  343.                      num_seaice_changes = num_seaice_changes + 1
    344. 

  344.                      if(landmask(i,j) .LT. 0.5 )tmn(i,j) = 271.4
    345. 

  345.                      vegcat(i,j)=isice
    346. 

  346.                      ivgtyp(i,j)=isice
    347. 

  347.                      lu_index(i,j)=isice
    348. 

  348.                      landmask(i,j)=1.
    349. 

  349.                      xland(i,j)=1.
    350. 

  350.                      vegfra(i,j)=0.
    351. 

  351.                      DO loop=1,num_veg_cat
    352. 

  352.                         landusef(i,loop,j)=0.
    353. 

  353.                      END DO
    354. 

  354.                      landusef(i,ivgtyp(i,j),j)=1.
    355. 

  355. 

  356. !tgs - compute blended sea ice/water skin temperature
    357. 

  357.                    if(flag_sst.eq.1) then
    358. 

  358.                      tsk(i,j) = xice(i,j)*(min(seaice_threshold,tsk(i,j)))  &
    359. 

  359.                                 +(1-xice(i,j))*sst(i,j)
    360. 

  360.                    else
    361. 

  361.                      tsk(i,j) = xice(i,j)*(min(seaice_threshold,tsk(i,j)))  &
    362. 

  362.                                 +(1-xice(i,j))*tsk(i,j)
    363. 

  363.                    endif
    364. 

  364.                      tsk_old(i,j) = tsk(i,j)
    365. 

  365. 

  366.                      isltyp(i,j) = 16
    367. 

  367.                      soilcat(i,j)=isltyp(i,j)
    368. 

  368.                      DO loop=1,num_soil_top_cat
    369. 

  369.                         soilctop(i,loop,j)=0
    370. 

  370.                      END DO
    371. 

  371.                      DO loop=1,num_soil_bot_cat
    372. 

  372.                         soilcbot(i,loop,j)=0
    373. 

  373.                      END DO
    374. 

  374.                      soilctop(i,isltyp(i,j),j)=1.
    375. 

  375.                      soilcbot(i,isltyp(i,j),j)=1.
    376. 

  376. 

  377.                  total_depth = 3. ! ice is 3 m deep, num_soil_layers equispaced layers
    378. 

  378.                        tslb(i,1,j) = tsk(i,j)
    379. 

  379.                        tslb(i,num_soil_layers,j) = tmn(i,j)
    380. 

  380.                      DO loop = 2,num_soil_layers-1
    381. 

  381.                         mid_point_depth=(total_depth/num_soil_layers)/4. + &
    382. 

  382.                                         (loop-2)*(total_depth/num_soil_layers)
    383. 

  383.                         tslb(i,loop,j) = ( (total_depth-mid_point_depth)*tsk(i,j) + &
    384. 

  384.                                             mid_point_depth*tmn(i,j) ) / total_depth
    385. 

  385.                      END DO
    386. 

  386. 

  387.                      DO loop=1,num_soil_layers
    388. 

  388.                         smois(i,loop,j) = 1.0
    389. 

  389.                         sh2o(i,loop,j)  = 0.0
    390. 

  390.                      END DO
    391. 

  391.                   ELSE IF ( xice(i,j) .LT. xice_threshold ) THEN
    392. 

  392.                      xice(i,j) = 0.
    393. 

  393.                   END IF
    394. 

  394.                END DO
    395. 

  395.             END DO
    396. 

  396.             IF ( num_seaice_changes .GT. 0 ) THEN
    397. 

  397.                WRITE ( message , FMT='(A,I6)' ) &
    398. 

  398.                'Total post number of sea ice location changes (water to land) = ', num_seaice_changes
    399. 

  399.                CALL wrf_debug ( 0 , message )
    400. 

  400.             END IF
    401. 

  401. 

  402. 

  403.       END SELECT fix_seaice
    404. 

  404. 

  405.    END SUBROUTINE adjust_for_seaice_post
    406. 

  406. 

  407.    SUBROUTINE process_percent_cat_new ( landmask ,  &
    408. 

  408.                                 landuse_frac , soil_top_cat , soil_bot_cat , &
    409. 

  409.                                 isltyp , ivgtyp , &
    410. 

  410.                                 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
    411. 

  411.                                 ids , ide , jds , jde , kds , kde , &
    412. 

  412.                                 ims , ime , jms , jme , kms , kme , &
    413. 

  413.                                 its , ite , jts , jte , kts , kte , &
    414. 

  414.                                 iswater )
    415. 

  415. 

  416.       IMPLICIT NONE
    417. 

  417. 

  418.       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
    419. 

  419.                               ims , ime , jms , jme , kms , kme , &
    420. 

  420.                               its , ite , jts , jte , kts , kte , &
    421. 

  421.                               iswater
    422. 

  422.       INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
    423. 

  423.       REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landuse_frac
    424. 

  424.       REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(IN):: soil_top_cat
    425. 

  425.       REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(IN):: soil_bot_cat
    426. 

  426.       INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
    427. 

  427.       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: landmask
    428. 

  428. 

  429.       INTEGER :: i , j , l , ll, dominant_index
    430. 

  430.       REAL :: dominant_value
    431. 

  431. 

  432. #ifdef WRF_CHEM
    433. 

  433. !      REAL :: lwthresh = .99
    434. 

  434.       REAL :: lwthresh = .50
    435. 

  435. #else
    436. 

  436.       REAL :: lwthresh = .50
    437. 

  437. #endif
    438. 

  438. 

  439.       INTEGER , PARAMETER :: iswater_soil = 14
    440. 

  440.       INTEGER :: iforce
    441. 

  441.       CHARACTER (LEN=132) :: message
    442. 

  442.       CHARACTER(LEN=256) :: mminlu
    443. 

  443.       LOGICAL :: aggregate_lu
    444. 

  444. integer :: change_water , change_land
    445. 

  445. change_water = 0
    446. 

  446. change_land = 0
    447. 

  447. 

  448.       !  Sanity check on the 50/50 points
    449. 

  449. 

  450.       DO j = jts , MIN(jde-1,jte)
    451. 

  451.          DO i = its , MIN(ide-1,ite)
    452. 

  452.             IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    453. 

  453.             dominant_value = landuse_frac(i,iswater,j)
    454. 

  454.             IF ( dominant_value .EQ. lwthresh ) THEN
    455. 

  455.                DO l = 1 , num_veg_cat
    456. 

  456.                   IF ( l .EQ. iswater ) CYCLE
    457. 

  457.                   IF ( ( landuse_frac(i,l,j) .EQ. lwthresh ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
    458. 

  458.                      PRINT *,i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
    459. 

  459.                      landuse_frac(i,l,j) = lwthresh - .01
    460. 

  460.                      landuse_frac(i,iswater,j) = lwthresh + 0.01
    461. 

  461.                   ELSE IF ( ( landuse_frac(i,l,j) .EQ. lwthresh ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
    462. 

  462.                      PRINT *,i,j,' water and category ',l,' both at 50%, landmask is ',landmask(i,j)
    463. 

  463.                      landuse_frac(i,l,j) = lwthresh + .01
    464. 

  464.                      landuse_frac(i,iswater,j) = lwthresh - 0.01
    465. 

  465.                   END IF
    466. 

  466.                END DO
    467. 

  467.             END IF
    468. 

  468.          END DO
    469. 

  469.       END DO
    470. 

  470. 

  471.       !  Compute the aggregate of the vegetation/land use categories.  Lump all of the grasses together,
    472. 

  472.       !  all of the shrubs, all of the trees, etc.  Choose the correct set of available land use
    473. 

  473.       !  categories.  Also, make sure that the user wants to actually avail themselves of aforementioned
    474. 

  474.       !  opportunity, as mayhaps they don't.
    475. 

  475.  
    476. 

  476.       CALL nl_get_mminlu       ( 1 , mminlu       )
    477. 

  477.       CALL nl_get_aggregate_lu ( 1 , aggregate_lu )
    478. 

  478.       IF ( aggregate_lu ) THEN
    479. 

  479.          DO j = jts , MIN(jde-1,jte)
    480. 

  480.             DO i = its , MIN(ide-1,ite)
    481. 

  481.                IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    482. 

  482.                CALL aggregate_categories_part1 ( landuse_frac , iswater , num_veg_cat , mminlu(1:4) )
    483. 

  483.             END DO
    484. 

  484.          END DO
    485. 

  485.       END IF
    486. 

  486. 

  487.       !  Compute the dominant VEGETATION INDEX.
    488. 

  488. 

  489.       DO j = jts , MIN(jde-1,jte)
    490. 

  490.          DO i = its , MIN(ide-1,ite)
    491. 

  491.             IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    492. 

  492.             dominant_value = landuse_frac(i,1,j)
    493. 

  493.             dominant_index = 1
    494. 

  494.             DO l = 2 , num_veg_cat
    495. 

  495.                IF        ( l .EQ. iswater ) THEN
    496. 

  496.                   ! wait a bit
    497. 

  497.                ELSE IF ( ( l .NE. iswater ) .AND. ( landuse_frac(i,l,j) .GT. dominant_value ) ) THEN
    498. 

  498.                   dominant_value = landuse_frac(i,l,j)
    499. 

  499.                   dominant_index = l
    500. 

  500.                END IF
    501. 

  501.             END DO
    502. 

  502.             IF ( landuse_frac(i,iswater,j) .GT. lwthresh ) THEN
    503. 

  503.                dominant_value = landuse_frac(i,iswater,j)
    504. 

  504.                dominant_index = iswater
    505. 

  505.             ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
    506. 

  506.                       ( landmask(i,j) .LT. 0.5) .AND. &
    507. 

  507.                       ( dominant_value .EQ. lwthresh) ) THEN
    508. 

  508.                dominant_value = landuse_frac(i,iswater,j)
    509. 

  509.                dominant_index = iswater
    510. 

  510.             ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh) .AND. &
    511. 

  511.                       ( landmask(i,j) .GT. 0.5) .AND. &
    512. 

  512.                       ( dominant_value .EQ. lwthresh) ) THEN
    513. 

  513.                !no op
    514. 

  514.             ELSE IF ( ( landuse_frac(i,iswater,j) .EQ. lwthresh ) .AND. &
    515. 

  515.                       ( dominant_value .LT. lwthresh ) ) THEN
    516. 

  516.                dominant_value = landuse_frac(i,iswater,j)
    517. 

  517.                dominant_index = iswater
    518. 

  518.             END IF
    519. 

  519.             IF      ( dominant_index .EQ. iswater ) THEN
    520. 

  520. if(landmask(i,j).gt.lwthresh) then
    521. 

  521. !print *,'changing to water at point ',i,j
    522. 

  522. !WRITE ( num_cat_count , FMT = '(I3)' ) num_veg_cat
    523. 

  523. !WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) ints(1:num_veg_cat)
    524. 

  524. !CALL wrf_debug(1,message)
    525. 

  525. !WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) nint(landuse_frac(i,:,j)*100)
    526. 

  526. !CALL wrf_debug(1,message)
    527. 

  527. change_water=change_water+1
    528. 

  528. endif
    529. 

  529.                landmask(i,j) = 0
    530. 

  530.             ELSE IF ( dominant_index .NE. iswater ) THEN
    531. 

  531. if(landmask(i,j).lt.lwthresh) then
    532. 

  532. !print *,'changing to land at point ',i,j
    533. 

  533. !WRITE ( num_cat_count , FMT = '(I3)' ) num_veg_cat
    534. 

  534. !WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) ints(1:num_veg_cat)
    535. 

  535. !CALL wrf_debug(1,message)
    536. 

  536. !WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) nint(landuse_frac(i,:,j)*100)
    537. 

  537. !CALL wrf_debug(1,message)
    538. 

  538. change_land=change_land+1
    539. 

  539. endif
    540. 

  540.                landmask(i,j) = 1
    541. 

  541.             END IF
    542. 

  542.             ivgtyp(i,j) = dominant_index
    543. 

  543.          END DO
    544. 

  544.       END DO
    545. 

  545. 

  546.       !  Compute the dominant SOIL TEXTURE INDEX, TOP.
    547. 

  547. 

  548.       iforce = 0
    549. 

  549.       DO i = its , MIN(ide-1,ite)
    550. 

  550.          DO j = jts , MIN(jde-1,jte)
    551. 

  551.             IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    552. 

  552.             dominant_value = soil_top_cat(i,1,j)
    553. 

  553.             dominant_index = 1
    554. 

  554.             IF ( landmask(i,j) .GT. lwthresh ) THEN
    555. 

  555.                DO l = 2 , num_soil_top_cat
    556. 

  556.                   IF ( ( l .NE. iswater_soil ) .AND. ( soil_top_cat(i,l,j) .GT. dominant_value ) ) THEN
    557. 

  557.                      dominant_value = soil_top_cat(i,l,j)
    558. 

  558.                      dominant_index = l
    559. 

  559.                   END IF
    560. 

  560.                END DO
    561. 

  561.                IF ( dominant_value .LT. 0.01 ) THEN
    562. 

  562.                   iforce = iforce + 1
    563. 

  563.                   WRITE ( message , FMT = '(A,I4,I4)' ) &
    564. 

  564.                   'based on landuse, changing soil to land at point ',i,j
    565. 

  565.                   CALL wrf_debug(1,message)
    566. 

  566.                   WRITE ( num_cat_count , FMT = '(I3)' ) num_soil_top_cat
    567. 

  567.                   WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) (ints(l),l=1,num_soil_top_cat)
    568. 

  568.                   CALL wrf_debug(1,message)
    569. 

  569.                   WRITE ( message , FMT = '('//num_cat_count//'(i3,1x))' ) &
    570. 

  570.                     ((nint(soil_top_cat(i,ints(l),j)*100)), l=1,num_soil_top_cat)
    571. 

  571.                   CALL wrf_debug(1,message)
    572. 

  572.                   dominant_index = 8
    573. 

  573.                END IF
    574. 

  574.             ELSE
    575. 

  575.                dominant_index = iswater_soil
    576. 

  576.             END IF
    577. 

  577.             isltyp(i,j) = dominant_index
    578. 

  578.          END DO
    579. 

  579.       END DO
    580. 

  580. 

  581. if(iforce.ne.0)then
    582. 

  582. WRITE(message,FMT='(A,I4,A,I6)' ) &
    583. 

  583. 'forcing artificial silty clay loam at ',iforce,' points, out of ',&
    584. 

  584. (MIN(ide-1,ite)-its+1)*(MIN(jde-1,jte)-jts+1)
    585. 

  585. CALL wrf_debug(0,message)
    586. 

  586. endif
    587. 

  587. print *,'LAND  CHANGE = ',change_land
    588. 

  588. print *,'WATER CHANGE = ',change_water
    589. 

  589. 

  590.    END SUBROUTINE process_percent_cat_new
    591. 

  591. 

  592.    SUBROUTINE process_soil_real ( tsk , tmn , tavgsfc, &
    593. 

  593.                                 landmask , sst , ht, toposoil, &
    594. 

  594.                                 st_input , sm_input , sw_input , &
    595. 

  595.                                 st_levels_input , sm_levels_input , sw_levels_input , &
    596. 

  596.                                 zs , dzs , tslb , smois , sh2o , &
    597. 

  597.                                 flag_sst , flag_tavgsfc, flag_soilhgt, &
    598. 

  598.                                 flag_soil_layers, flag_soil_levels, &
    599. 

  599.                                 ids , ide , jds , jde , kds , kde , &
    600. 

  600.                                 ims , ime , jms , jme , kms , kme , &
    601. 

  601.                                 its , ite , jts , jte , kts , kte , &
    602. 

  602.                                 sf_surface_physics , num_soil_layers , real_data_init_type , &
    603. 

  603.                                 num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
    604. 

  604.                                 num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc )
    605. 

  605. 

  606.       IMPLICIT NONE
    607. 

  607. 

  608.       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
    609. 

  609.                               ims , ime , jms , jme , kms , kme , &
    610. 

  610.                               its , ite , jts , jte , kts , kte , &
    611. 

  611.                               sf_surface_physics , num_soil_layers , real_data_init_type , &
    612. 

  612.                               num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
    613. 

  613.                               num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc
    614. 

  614. 

  615.       INTEGER , INTENT(IN) :: flag_sst, flag_tavgsfc
    616. 

  616.       INTEGER , INTENT(IN) :: flag_soil_layers, flag_soil_levels, flag_soilhgt
    617. 

  617. 

  618.       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: landmask , sst
    619. 

  619. 

  620.       INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(INOUT) :: st_levels_input
    621. 

  621.       INTEGER , DIMENSION(1:num_sm_levels_input) , INTENT(INOUT) :: sm_levels_input
    622. 

  622.       INTEGER , DIMENSION(1:num_sw_levels_input) , INTENT(INOUT) :: sw_levels_input
    623. 

  623.       REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
    624. 

  624.       REAL , DIMENSION(ims:ime,1:num_sm_levels_alloc,jms:jme) , INTENT(INOUT) :: sm_input
    625. 

  625.       REAL , DIMENSION(ims:ime,1:num_sw_levels_alloc,jms:jme) , INTENT(INOUT) :: sw_input
    626. 

  626. 

  627.       REAL, DIMENSION(1:num_soil_layers), INTENT(OUT)  ::  zs,dzs
    628. 

  628.       REAL , DIMENSION(ims:ime,num_soil_layers,jms:jme) , INTENT(OUT) :: tslb , smois , sh2o
    629. 

  629.       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN) :: tavgsfc, ht, toposoil
    630. 

  630.       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tsk, tmn
    631. 

  631. 

  632.       INTEGER :: i , j , k, l , dominant_index , num_soil_cat , num_veg_cat, closest_layer
    633. 

  633.       REAL :: dominant_value, closest_depth, diff_cm
    634. 

  634.       REAL , ALLOCATABLE , DIMENSION(:) :: depth     ! Soil layer thicknesses (cm)
    635. 

  635.       REAL, PARAMETER :: get_temp_closest_to = 30.   ! use soil temperature closest to this depth (cm)
    636. 

  636.       REAL, PARAMETER :: something_big = 1.e6        ! Initialize closest depth as something big (cm)
    637. 

  637.       INTEGER :: something_far = 1000                ! Soil array index far away
    638. 

  638.       CHARACTER (LEN=132) :: message
    639. 

  639. 

  640.       !  Case statement for tmn initialization
    641. 

  641.       !  Need to have a reasonable default value for annual mean deeeeep soil temperature
    642. 

  642.       !  For sf_surface_physics = 1, we want to use close to a 30 cm value
    643. 

  643.       !  for the bottom level of the soil temps.
    644. 

  644.       ! NOTE:  We are assuming that soil_layers are the same for each grid point
    645. 

  645. 

  646.       fix_bottom_level_for_temp : SELECT CASE ( sf_surface_physics )
    647. 

  647.          CASE (SLABSCHEME)
    648. 

  648.             IF ( flag_tavgsfc  .EQ. 1 ) THEN
    649. 

  649.                CALL wrf_debug ( 0 , 'Using average surface temperature for tmn')
    650. 

  650.                DO j = jts , MIN(jde-1,jte)
    651. 

  651.                   DO i = its , MIN(ide-1,ite)
    652. 

  652.                      IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    653. 

  653.                      tmn(i,j) = tavgsfc(i,j)
    654. 

  654.                   END DO
    655. 

  655.                END DO
    656. 

  656.             ELSE
    657. 

  657.             ! Look for soil temp close to 30 cm
    658. 

  658.                closest_layer = something_far
    659. 

  659.                closest_depth = something_big
    660. 

  660.                DO k = 1, num_st_levels_input
    661. 

  661.                   diff_cm = abs( st_levels_input(k) - get_temp_closest_to )
    662. 

  662.                   IF ( diff_cm < closest_depth ) THEN
    663. 

  663.                      closest_depth = diff_cm
    664. 

  664.                      closest_layer = k
    665. 

  665.                   END IF
    666. 

  666.                END DO
    667. 

  667.                IF ( closest_layer == something_far ) THEN
    668. 

  668.                   CALL wrf_debug ( 0 , 'No soil temperature data for grid%tmn near 30 cm')
    669. 

  669.                   CALL wrf_debug ( 0 , 'Using 1 degree static annual mean temps' )
    670. 

  670.                ELSE
    671. 

  671.                   write(message, FMT='(A,F7.2,A,I3)')&
    672. 

  672.                      'Soil temperature closest to ',get_temp_closest_to, &
    673. 

  673.                      ' at level ',st_levels_input(closest_layer)
    674. 

  674.                   CALL wrf_debug ( 0 , message )
    675. 

  675.                   DO j = jts , MIN(jde-1,jte)
    676. 

  676.                      DO i = its , MIN(ide-1,ite)
    677. 

  677.                         IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    678. 

  678.                         tmn(i,j) = st_input(i,closest_layer+1,j)
    679. 

  679.                      END DO
    680. 

  680.                   END DO
    681. 

  681.                END IF
    682. 

  682.             END IF
    683. 

  683. 

  684.          CASE (LSMSCHEME)
    685. 

  685. 

  686.          CASE (NOAHMPSCHEME)
    687. 

  687. 

  688.          CASE (RUCLSMSCHEME)
    689. 

  689. 

  690.          CASE (PXLSMSCHEME)
    691. 

  691. ! When the input data from met_em is in layers, there is an additional level added to the beginning
    692. 

  692. ! of the array to define the surface, which is why we add the extra value (closest_layer+1)
    693. 

  693.             IF ( flag_tavgsfc  .EQ. 1 ) THEN
    694. 

  694.                CALL wrf_debug ( 0 , 'Using average surface temperature for tmn')
    695. 

  695.                DO j = jts , MIN(jde-1,jte)
    696. 

  696.                   DO i = its , MIN(ide-1,ite)
    697. 

  697.                      IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    698. 

  698.                      tmn(i,j) = tavgsfc(i,j)
    699. 

  699.                   END DO
    700. 

  700.                END DO
    701. 

  701.             ELSE
    702. 

  702.             ! Look for soil temp close to 30 cm
    703. 

  703.                closest_layer = num_st_levels_input+1
    704. 

  704.                closest_depth = something_big
    705. 

  705.                DO k = 1, num_st_levels_input
    706. 

  706.                   diff_cm = abs( st_levels_input(k) - get_temp_closest_to )
    707. 

  707.                   IF ( diff_cm < closest_depth ) THEN
    708. 

  708.                      closest_depth = diff_cm
    709. 

  709.                      closest_layer = k
    710. 

  710.                   END IF
    711. 

  711.                END DO
    712. 

  712.                IF ( closest_layer == num_st_levels_input + 1 ) THEN
    713. 

  713.                   CALL wrf_debug ( 0 , 'No soil temperature data for grid%tmn near 30 cm')
    714. 

  714.                   CALL wrf_debug ( 0 , 'Using 1 degree static annual mean temps' )
    715. 

  715.                ELSE
    716. 

  716.                   write(message, FMT='(A,F7.2,A,I3)')&
    717. 

  717.                      'Soil temperature closest to ',get_temp_closest_to, &
    718. 

  718.                      ' at level ',st_levels_input(closest_layer)
    719. 

  719.                   CALL wrf_debug ( 0 , message )
    720. 

  720.                   DO j = jts , MIN(jde-1,jte)
    721. 

  721.                      DO i = its , MIN(ide-1,ite)
    722. 

  722.                         IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    723. 

  723.                         tmn(i,j) = st_input(i,closest_layer+1,j)
    724. 

  724.                      END DO
    725. 

  725.                   END DO
    726. 

  726.                END IF
    727. 

  727.             END IF
    728. 

  728. 

  729. #if 0
    730. 

  730.             ! Loop over layers and do a weighted mean
    731. 

  731.             IF ( ALLOCATED ( depth ) ) DEALLOCATE ( depth )
    732. 

  732.             ALLOCATE ( depth(num_st_levels_input) )
    733. 

  733.             IF ( flag_soil_layers == 1 ) THEN
    734. 

  734.                DO k = num_st_levels_input, 2, -1
    735. 

  735.                   depth(k) = st_levels_input(k) - st_levels_input(k-1)
    736. 

  736.                END DO
    737. 

  737.                depth(1) = st_levels_input(1)
    738. 

  738.             ELSE IF ( flag_soil_levels == 1 ) THEN
    739. 

  739.                DO k = 2, num_st_levels_input
    740. 

  740.                   depth(k) = st_levels_input(k) - st_levels_input(k-1)
    741. 

  741.                END DO
    742. 

  742.                depth(1) = 0.
    743. 

  743.             END IF
    744. 

  744. 

  745.             DO j = jts , MIN(jde-1,jte)
    746. 

  746.                DO i = its , MIN(ide-1,ite)
    747. 

  747.                   IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    748. 

  748.                   tmn(i,j) = 0.
    749. 

  749.                   DO k = 1, num_st_levels_input
    750. 

  750.                      tmn(i,j) = tmn(i,j) + depth(k) * st_input(i,k,j)
    751. 

  751.                   END DO
    752. 

  752.                END DO
    753. 

  753.             END DO
    754. 

  754.             DEALLOCATE ( depth )
    755. 

  755. 

  756. #endif
    757. 

  757. 

  758.       END SELECT fix_bottom_level_for_temp
    759. 

  759. 

  760.       !  Adjust the various soil temperature values depending on the difference in
    761. 

  761.       !  elevation between the current model's elevation and the incoming data's
    762. 

  762.       !  orography.
    763. 

  763. 

  764.       adjust_soil : SELECT CASE ( sf_surface_physics )
    765. 

  765. 

  766.          CASE ( SLABSCHEME , LSMSCHEME , NOAHMPSCHEME , RUCLSMSCHEME, PXLSMSCHEME, SSIBSCHEME )
    767. 

  767.             CALL adjust_soil_temp_new ( tmn , sf_surface_physics , tsk , ht ,            &
    768. 

  768.                                         toposoil , landmask , st_input, st_levels_input, &
    769. 

  769.                                         flag_soilhgt , flag_tavgsfc ,                    &
    770. 

  770.                                         flag_soil_layers , flag_soil_levels,             &
    771. 

  771.                                         num_st_levels_input, num_st_levels_alloc,        &
    772. 

  772.                                         ids , ide , jds , jde , kds , kde ,              &
    773. 

  773.                                         ims , ime , jms , jme , kms , kme ,              &
    774. 

  774.                                         its , ite , jts , jte , kts , kte )
    775. 

  775. 

  776.       END SELECT adjust_soil
    777. 

  777. 

  778.       !  Initialize the soil depth, and the soil temperature and moisture.
    779. 

  779.    
    780. 

  780.       IF      ( ( sf_surface_physics .EQ. SLABSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
    781. 

  781.          CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
    782. 

  782.          CALL init_soil_1_real ( tsk , tmn , tslb , zs , dzs , num_soil_layers , real_data_init_type , &
    783. 

  783.                                  landmask , sst , flag_sst , &
    784. 

  784.                                  ids , ide , jds , jde , kds , kde , &
    785. 

  785.                                  ims , ime , jms , jme , kms , kme , &
    786. 

  786.                                  its , ite , jts , jte , kts , kte )
    787. 

  787. 

  788.       ELSE IF ( ( sf_surface_physics .EQ. LSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
    789. 

  789.          CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
    790. 

  790.          CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
    791. 

  791.                                  st_input , sm_input , sw_input , landmask , sst , &
    792. 

  792.                                  zs , dzs , &
    793. 

  793.                                  st_levels_input , sm_levels_input , sw_levels_input , &
    794. 

  794.                                  num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
    795. 

  795.                                  num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
    796. 

  796.                                  flag_sst , flag_soil_layers , flag_soil_levels , &
    797. 

  797.                                  ids , ide , jds , jde , kds , kde , &
    798. 

  798.                                  ims , ime , jms , jme , kms , kme , &
    799. 

  799.                                  its , ite , jts , jte , kts , kte )
    800. 

  800.       ELSE IF ( ( sf_surface_physics .EQ. NOAHMPSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
    801. 

  801.          CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
    802. 

  802.          CALL init_soil_2_real ( tsk , tmn , smois , sh2o , tslb , &
    803. 

  803.                                  st_input , sm_input , sw_input , landmask , sst , &
    804. 

  804.                                  zs , dzs , &
    805. 

  805.                                  st_levels_input , sm_levels_input , sw_levels_input , &
    806. 

  806.                                  num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
    807. 

  807.                                  num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
    808. 

  808.                                  flag_sst , flag_soil_layers , flag_soil_levels , &
    809. 

  809.                                  ids , ide , jds , jde , kds , kde , &
    810. 

  810.                                  ims , ime , jms , jme , kms , kme , &
    811. 

  811.                                  its , ite , jts , jte , kts , kte )
    812. 

  812.       ELSE IF ( ( sf_surface_physics .EQ. RUCLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
    813. 

  813.          CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
    814. 

  814.          CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
    815. 

  815.                                  st_input , sm_input , landmask , sst , &
    816. 

  816.                                  zs , dzs , &
    817. 

  817.                                  st_levels_input , sm_levels_input , &
    818. 

  818.                                  num_soil_layers , num_st_levels_input , num_sm_levels_input , &
    819. 

  819.                                  num_st_levels_alloc , num_sm_levels_alloc , &
    820. 

  820.                                  flag_sst , flag_soil_layers , flag_soil_levels , &
    821. 

  821.                                  ids , ide , jds , jde , kds , kde , &
    822. 

  822.                                  ims , ime , jms , jme , kms , kme , &
    823. 

  823.                                  its , ite , jts , jte , kts , kte )
    824. 

  824.       ELSE IF ( ( sf_surface_physics .EQ. PXLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
    825. 

  825.          CALL init_soil_depth_7 ( zs , dzs , num_soil_layers )
    826. 

  826.          CALL init_soil_7_real ( tsk , tmn , smois , sh2o, tslb , &
    827. 

  827.                                  st_input , sm_input , sw_input, landmask , sst , &
    828. 

  828.                                  zs , dzs , &
    829. 

  829.                                  st_levels_input , sm_levels_input , sw_levels_input, &
    830. 

  830.                                  num_soil_layers , num_st_levels_input , num_sm_levels_input , num_sw_levels_input , &
    831. 

  831.                                  num_st_levels_alloc , num_sm_levels_alloc , num_sw_levels_alloc , &
    832. 

  832.                                  flag_sst , flag_soil_layers , flag_soil_levels , &
    833. 

  833.                                  ids , ide , jds , jde , kds , kde , &
    834. 

  834.                                  ims , ime , jms , jme , kms , kme , &
    835. 

  835.                                  its , ite , jts , jte , kts , kte )
    836. 

  836.       ELSE IF ( ( sf_surface_physics .EQ. 8 ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
    837. 

  837.          CALL init_soil_depth_8 ( zs , dzs , num_soil_layers )
    838. 

  838.          CALL init_soil_3_real ( tsk , tmn , smois , tslb , &
    839. 

  839.                                  st_input , sm_input , landmask , sst , &
    840. 

  840.                                  zs , dzs , &
    841. 

  841.                                  st_levels_input , sm_levels_input , &
    842. 

  842.                                  num_soil_layers , num_st_levels_input , num_sm_levels_input , &
    843. 

  843.                                  num_st_levels_alloc , num_sm_levels_alloc , &
    844. 

  844.                                  flag_sst , flag_soil_layers , flag_soil_levels , &
    845. 

  845.                                  ids , ide , jds , jde , kds , kde , &
    846. 

  846.                                  ims , ime , jms , jme , kms , kme , &
    847. 

  847.                                  its , ite , jts , jte , kts , kte )
    848. 

  848.       END IF
    849. 

  849. 

  850.    END SUBROUTINE process_soil_real
    851. 

  851. 

  852.    SUBROUTINE process_soil_ideal ( xland,xice,vegfra,snow,canwat,  &
    853. 

  853.                                    ivgtyp,isltyp,tslb,smois, &
    854. 

  854.                                    tsk,tmn,zs,dzs,           &
    855. 

  855.                                    num_soil_layers,          &
    856. 

  856.                                    sf_surface_physics ,      &
    857. 

  857.                                    ids,ide, jds,jde, kds,kde,&
    858. 

  858.                                    ims,ime, jms,jme, kms,kme,&
    859. 

  859.                                    its,ite, jts,jte, kts,kte )
    860. 

  860. 

  861.       IMPLICIT NONE
    862. 

  862. 

  863.       INTEGER, INTENT(IN) ::ids,ide, jds,jde, kds,kde,  &
    864. 

  864.                             ims,ime, jms,jme, kms,kme,  &
    865. 

  865.                             its,ite, jts,jte, kts,kte
    866. 

  866. 

  867.       INTEGER, INTENT(IN) :: num_soil_layers , sf_surface_physics
    868. 

  868. 

  869.       REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , INTENT(INOUT) :: smois, tslb
    870. 

  870. 

  871.       REAL, DIMENSION(num_soil_layers), INTENT(OUT) :: dzs,zs
    872. 

  872. 

  873.       REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(INOUT) :: tsk, tmn
    874. 

  874.       REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: xland, snow, canwat, xice, vegfra
    875. 

  875.       INTEGER, DIMENSION( ims:ime, jms:jme ) , INTENT(OUT) :: ivgtyp, isltyp
    876. 

  876. 

  877.       !  Local variables.
    878. 

  878. 

  879.       INTEGER :: itf,jtf
    880. 

  880. 

  881.       itf=MIN(ite,ide-1)
    882. 

  882.       jtf=MIN(jte,jde-1)
    883. 

  883. 

  884.       IF      ( ( sf_surface_physics .EQ. SLABSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
    885. 

  885.          CALL init_soil_depth_1 ( zs , dzs , num_soil_layers )
    886. 

  886.          CALL init_soil_1_ideal(tsk,tmn,tslb,xland,                      &
    887. 

  887.                                 ivgtyp,zs,dzs,num_soil_layers,           &
    888. 

  888.                                 ids,ide, jds,jde, kds,kde,               &
    889. 

  889.                                 ims,ime, jms,jme, kms,kme,               &
    890. 

  890.                                 its,ite, jts,jte, kts,kte                )
    891. 

  891.       ELSE IF ( ( sf_surface_physics .EQ. LSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
    892. 

  892.          CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
    893. 

  893.          CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat,         &
    894. 

  894.                                   ivgtyp,isltyp,tslb,smois,tmn,          &
    895. 

  895.                                   num_soil_layers,                       &
    896. 

  896.                                   ids,ide, jds,jde, kds,kde,             &
    897. 

  897.                                   ims,ime, jms,jme, kms,kme,             &
    898. 

  898.                                   its,ite, jts,jte, kts,kte              )
    899. 

  899.       ELSE IF ( ( sf_surface_physics .EQ. NOAHMPSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
    900. 

  900.          CALL init_soil_depth_2 ( zs , dzs , num_soil_layers )
    901. 

  901.          CALL init_soil_2_ideal ( xland,xice,vegfra,snow,canwat,         &
    902. 

  902.                                   ivgtyp,isltyp,tslb,smois,tmn,          &
    903. 

  903.                                   num_soil_layers,                       &
    904. 

  904.                                   ids,ide, jds,jde, kds,kde,             &
    905. 

  905.                                   ims,ime, jms,jme, kms,kme,             &
    906. 

  906.                                   its,ite, jts,jte, kts,kte              )
    907. 

  907.       ELSE IF ( ( sf_surface_physics .EQ. RUCLSMSCHEME ) .AND. ( num_soil_layers .GT. 1 ) ) THEN
    908. 

  908.          CALL init_soil_depth_3 ( zs , dzs , num_soil_layers )
    909. 

  909. 

  910.       END IF
    911. 

  911. 

  912.    END SUBROUTINE process_soil_ideal
    913. 

  913. 

  914.    SUBROUTINE adjust_soil_temp_new ( tmn , sf_surface_physics , tsk , ter ,            &
    915. 

  915.                                      toposoil , landmask , st_input , st_levels_input, &
    916. 

  916.                                      flag_toposoil , flag_tavgsfc ,                    &
    917. 

  917.                                      flag_soil_layers , flag_soil_levels,              &
    918. 

  918.                                      num_st_levels_input, num_st_levels_alloc,         &
    919. 

  919.                                      ids , ide , jds , jde , kds , kde ,               &
    920. 

  920.                                      ims , ime , jms , jme , kms , kme ,               &
    921. 

  921.                                      its , ite , jts , jte , kts , kte )
    922. 

  922. 

  923.       IMPLICIT NONE
    924. 

  924. 

  925.       INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
    926. 

  926.                               ims , ime , jms , jme , kms , kme , &
    927. 

  927.                               its , ite , jts , jte , kts , kte 
    928. 

  928.       INTEGER , INTENT(IN) :: num_st_levels_input, num_st_levels_alloc
    929. 

  929. 

  930.       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(IN)    :: ter , toposoil , landmask
    931. 

  931.       REAL , DIMENSION(ims:ime,jms:jme) , INTENT(INOUT) :: tmn , tsk
    932. 

  932.       REAL , DIMENSION(ims:ime,1:num_st_levels_alloc,jms:jme) , INTENT(INOUT) :: st_input
    933. 

  933.       INTEGER , DIMENSION(1:num_st_levels_input) , INTENT(IN) :: st_levels_input
    934. 

  934. 

  935.       INTEGER , INTENT(IN) :: sf_surface_physics , flag_toposoil , flag_tavgsfc
    936. 

  936.       INTEGER , INTENT(IN) :: flag_soil_layers , flag_soil_levels
    937. 

  937.  
    938. 

  938.       INTEGER :: i , j, k , st_near_sfc
    939. 

  939. 

  940.       REAL :: soil_elev_min_val ,  soil_elev_max_val , soil_elev_min_dif , soil_elev_max_dif
    941. 

  941. 

  942.       !  Adjust the annual mean temperature as if it is based on from a sea-level elevation
    943. 

  943.       !  if the value used is from the actual annula mean data set.  If the input field to
    944. 

  944.       !  be used for tmn is one of the first-guess input temp fields, need to do an adjustment
    945. 

  945.       !  only on the diff in topo from the model terrain and the first-guess terrain.
    946. 

  946. 

  947.        SELECT CASE ( sf_surface_physics )
    948. 

  948. 

  949.          CASE ( LSMSCHEME , NOAHMPSCHEME )
    950. 

  950.             DO j = jts , MIN(jde-1,jte)
    951. 

  951.                DO i = its , MIN(ide-1,ite)
    952. 

  952.                   IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    953. 

  953.                   IF (landmask(i,j) .GT. 0.5 ) THEN
    954. 

  954.                      tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
    955. 

  955.                   END IF
    956. 

  956.                END DO
    957. 

  957.             END DO
    958. 

  958. 

  959.          CASE (RUCLSMSCHEME)
    960. 

  960.             DO j = jts , MIN(jde-1,jte)
    961. 

  961.                DO i = its , MIN(ide-1,ite)
    962. 

  962.                   IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    963. 

  963.                   IF (landmask(i,j) .GT. 0.5 ) THEN
    964. 

  964.                      tmn(i,j) = tmn(i,j) - 0.0065 * ter(i,j)
    965. 

  965.                   END IF
    966. 

  966.                END DO
    967. 

  967.             END DO
    968. 

  968. 

  969.       END SELECT
    970. 

  970. 

  971. 

  972.       !  Do we have a soil field with which to modify soil temperatures?
    973. 

  973. 

  974.       IF ( flag_toposoil .EQ. 1 ) THEN
    975. 

  975. 

  976.          DO j = jts , MIN(jde-1,jte)
    977. 

  977.             DO i = its , MIN(ide-1,ite)
    978. 

  978. 

  979.                IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE 
    980. 

  980.                !  Is the toposoil field OK, or is it a subversive soil elevation field.  We can tell
    981. 

  981.                !  usually by looking at values.  Anything less than -1000 m (lower than the Dead Sea) is
    982. 

  982.                !  bad.  Anything larger than 10 km (taller than Everest) is toast.  Also, anything where
    983. 

  983.                !  the difference between the soil elevation and the terrain is greater than 3 km means
    984. 

  984.                !  that the soil data is either all zeros or that the data are inconsistent.  Any of these
    985. 

  985.                !  three conditions is grievous enough to induce a WRF fatality.  However, if they are at
    986. 

  986.                !  a water point, then we can safely ignore them.
    987. 

  987. 

  988.                soil_elev_min_val = toposoil(i,j)
    989. 

  989.                soil_elev_max_val = toposoil(i,j)
    990. 

  990.                soil_elev_min_dif = ter(i,j) - toposoil(i,j)
    991. 

  991.                soil_elev_max_dif = ter(i,j) - toposoil(i,j)
    992. 

  992. 

  993.                IF      ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .LT. 0.5 ) ) THEN
    994. 

  994.                   CYCLE
    995. 

  995.                ELSE IF ( ( soil_elev_min_val .LT. -1000 ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
    996. 

  996. !print *,'no soil temperature elevation adjustment, soil height too small = ',toposoil(i,j)
    997. 

  997. cycle
    998. 

  998. !                 CALL wrf_error_f…
    999.

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