PageRenderTime 62ms CodeModel.GetById 17ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/dyn_nmm/module_initialize_tropical_cyclone.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 6633 lines | 4328 code | 1076 blank | 1229 comment | 142 complexity | 9d9f77a97607ac1b632315424e19711f 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. !REAL:MODEL_LAYER:INITIALIZATION
  2. ! This MODULE holds the routines which are used to perform various initializations
  3. ! for individual domains utilizing the NMM dynamical core.
  4. !-----------------------------------------------------------------------
  5. MODULE module_initialize_ideal
  6. USE module_bc
  7. USE module_configure
  8. USE module_domain
  9. USE module_io_domain
  10. USE module_model_constants
  11. ! USE module_si_io_nmm
  12. USE module_state_description
  13. USE module_timing
  14. USE module_soil_pre
  15. #ifdef DM_PARALLEL
  16. USE module_dm, ONLY : LOCAL_COMMUNICATOR &
  17. ,MYTASK,NTASKS,NTASKS_X &
  18. ,NTASKS_Y
  19. USE module_comm_dm
  20. USE module_ext_internal
  21. #endif
  22. INTEGER :: internal_time_loop
  23. INTEGER:: MPI_COMM_COMP,MYPE
  24. INTEGER:: loopinc, iloopinc
  25. CONTAINS
  26. !-------------------------------------------------------------------
  27. SUBROUTINE init_domain ( grid )
  28. IMPLICIT NONE
  29. ! Input space and data. No gridded meteorological data has been stored, though.
  30. ! TYPE (domain), POINTER :: grid
  31. TYPE (domain) :: grid
  32. ! Local data.
  33. INTEGER :: idum1, idum2
  34. CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
  35. CALL init_domain_nmm (grid &
  36. !
  37. #include <actual_new_args.inc>
  38. !
  39. )
  40. END SUBROUTINE init_domain
  41. !-------------------------------------------------------------------
  42. !---------------------------------------------------------------------
  43. SUBROUTINE init_domain_nmm ( grid &
  44. !
  45. # include <dummy_new_args.inc>
  46. !
  47. )
  48. USE module_optional_input
  49. IMPLICIT NONE
  50. ! Input space and data. No gridded meteorological data has been stored, though.
  51. ! TYPE (domain), POINTER :: grid
  52. TYPE (domain) :: grid
  53. # include <dummy_new_decl.inc>
  54. real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_metgrid_levels) :: ght_out
  55. real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_metgrid_levels) :: rh_out
  56. real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_metgrid_levels) :: t_out
  57. real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_metgrid_levels) :: u_out
  58. real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_metgrid_levels) :: v_out
  59. real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_land_cat) :: landusef_out
  60. real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_soil_cat) :: soilcbot_out
  61. real ,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,1:grid%num_soil_cat) :: soilctop_out
  62. TYPE (grid_config_rec_type) :: config_flags
  63. ! Local domain indices and counters.
  64. INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
  65. INTEGER :: num_veg_gc , num_soil_top_gc , num_soil_bot_gc
  66. INTEGER :: &
  67. ids, ide, jds, jde, kds, kde, &
  68. ims, ime, jms, jme, kms, kme, &
  69. its, ite, jts, jte, kts, kte, &
  70. ips, ipe, jps, jpe, kps, kpe, &
  71. i, j, k, NNXP, NNYP
  72. ! Local data
  73. CHARACTER(LEN=19):: start_date
  74. #ifdef DM_PARALLEL
  75. LOGICAL,EXTERNAL :: WRF_DM_ON_MONITOR
  76. logical :: test
  77. ! INTEGER :: DOMDESC
  78. REAL,ALLOCATABLE :: SICE_G(:,:), SM_G(:,:)
  79. INTEGER, ALLOCATABLE:: IHE_G(:),IHW_G(:)
  80. INTEGER, ALLOCATABLE:: ITEMP(:,:)
  81. INTEGER :: my_e,my_n,my_s,my_w,my_ne,my_nw,my_se,my_sw,myi,myj,npe
  82. INTEGER :: istat,inpes,jnpes
  83. #else
  84. integer, allocatable:: ihw(:),ihe(:)
  85. #endif
  86. CHARACTER (LEN=255) :: message
  87. INTEGER :: error
  88. REAL :: p_surf, p_level
  89. REAL :: cof1, cof2
  90. REAL :: qvf , qvf1 , qvf2 , pd_surf
  91. REAL :: p00 , t00 , a
  92. REAL :: hold_znw, rmin,rmax
  93. REAL :: p_top_requested , ptsgm
  94. INTEGER :: num_metgrid_levels, ICOUNT
  95. REAL , DIMENSION(max_eta) :: eta_levels
  96. LOGICAL :: stretch_grid, dry_sounding, debug, log_flag_sst, hyb_coor
  97. REAL, ALLOCATABLE,DIMENSION(:,:):: ADUM2D,SNOWC,HT,TG_ALT, &
  98. PDVP,PSFC_OUTV
  99. REAL, ALLOCATABLE,DIMENSION(:,:,:):: P3D_OUT,P3DV_OUT,P3DV_IN, &
  100. QTMP,QTMP2
  101. INTEGER, ALLOCATABLE, DIMENSION(:):: KHL2,KVL2,KHH2,KVH2, &
  102. KHLA,KHHA,KVLA,KVHA
  103. ! INTEGER, ALLOCATABLE, DIMENSION(:,:):: grid%lu_index
  104. REAL, ALLOCATABLE, DIMENSION(:):: DXJ,WPDARJ,CPGFUJ,CURVJ, &
  105. FCPJ,FDIVJ,EMJ,EMTJ,FADJ, &
  106. HDACJ,DDMPUJ,DDMPVJ
  107. REAL, ALLOCATABLE,DIMENSION(:),SAVE:: SG1,SG2,DSG1,DSG2, &
  108. SGML1,SGML2
  109. !-- Carsel and Parrish [1988]
  110. REAL , DIMENSION(100) :: lqmi
  111. integer iicount
  112. REAL:: TPH0D,TLM0D
  113. REAL:: TPH0,WB,SB,TDLM,TDPH
  114. REAL:: WBI,SBI,EBI,ANBI,STPH0,CTPH0
  115. REAL:: TSPH,DTAD,DTCF
  116. REAL:: ACDT,CDDAMP,DXP,FP
  117. REAL:: WBD,SBD
  118. REAL:: RSNOW,SNOFAC
  119. REAL, PARAMETER:: SALP=2.60
  120. REAL, PARAMETER:: SNUP=0.040
  121. REAL:: SMCSUM,STCSUM,SEAICESUM,FISX
  122. REAL:: cur_smc, aposs_smc
  123. REAL:: COAC, CODAMP
  124. INTEGER,PARAMETER:: DOUBLE=SELECTED_REAL_KIND(15,300)
  125. REAL(KIND=DOUBLE):: TERM1,APH,TLM,TPH,DLM,DPH,STPH,CTPH
  126. INTEGER:: KHH,KVH,JAM,JA, IHL, IHH, L
  127. INTEGER:: II,JJ,ISRCH,ISUM,ITER,Ilook,Jlook
  128. REAL, PARAMETER:: DTR=0.01745329
  129. REAL, PARAMETER:: W_NMM=0.08
  130. #if defined(HWRF)
  131. REAL, PARAMETER:: DDFC=1.0
  132. #else
  133. REAL, PARAMETER:: DDFC=8.0
  134. #endif
  135. REAL, PARAMETER:: TWOM=.00014584
  136. REAL, PARAMETER:: CP=1004.6
  137. REAL, PARAMETER:: DFC=1.0
  138. REAL, PARAMETER:: ROI=916.6
  139. REAL, PARAMETER:: R=287.04
  140. REAL, PARAMETER:: CI=2060.0
  141. REAL, PARAMETER:: ROS=1500.
  142. REAL, PARAMETER:: CS=1339.2
  143. REAL, PARAMETER:: DS=0.050
  144. REAL, PARAMETER:: AKS=.0000005
  145. REAL, PARAMETER:: DZG=2.85
  146. REAL, PARAMETER:: DI=.1000
  147. REAL, PARAMETER:: AKI=0.000001075
  148. REAL, PARAMETER:: DZI=2.0
  149. REAL, PARAMETER:: THL=210.
  150. REAL, PARAMETER:: PLQ=70000.
  151. REAL, PARAMETER:: ERAD=6371200.
  152. REAL, PARAMETER:: TG0=258.16
  153. REAL, PARAMETER:: TGA=30.0
  154. integer :: numzero,numexamined
  155. #ifdef HWRF
  156. !============================================================================
  157. ! gopal's doing for ocean coupling
  158. !============================================================================
  159. REAL, DIMENSION(:,:), ALLOCATABLE :: NHLAT,NHLON,NVLAT,NVLON,HRES_SM
  160. REAL :: NDLMD,NDPHD,NWBD,NSBD
  161. INTEGER :: NIDE,NJDE,ILOC,JLOC
  162. INTEGER fid, ierr, nprocs
  163. CHARACTER*255 f65name, SysString
  164. !============================================================================
  165. ! end gopal's doing for ocean coupling
  166. !============================================================================
  167. #endif
  168. if (ALLOCATED(ADUM2D)) DEALLOCATE(ADUM2D)
  169. if (ALLOCATED(TG_ALT)) DEALLOCATE(TG_ALT)
  170. !#define COPY_IN
  171. !#include <scalar_derefs.inc>
  172. #ifdef DM_PARALLEL
  173. # include <data_calls.inc>
  174. #endif
  175. SELECT CASE ( model_data_order )
  176. CASE ( DATA_ORDER_ZXY )
  177. kds = grid%sd31 ; kde = grid%ed31 ;
  178. ids = grid%sd32 ; ide = grid%ed32 ;
  179. jds = grid%sd33 ; jde = grid%ed33 ;
  180. kms = grid%sm31 ; kme = grid%em31 ;
  181. ims = grid%sm32 ; ime = grid%em32 ;
  182. jms = grid%sm33 ; jme = grid%em33 ;
  183. kts = grid%sp31 ; kte = grid%ep31 ; ! tile is entire patch
  184. its = grid%sp32 ; ite = grid%ep32 ; ! tile is entire patch
  185. jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
  186. CASE ( DATA_ORDER_XYZ )
  187. ids = grid%sd31 ; ide = grid%ed31 ;
  188. jds = grid%sd32 ; jde = grid%ed32 ;
  189. kds = grid%sd33 ; kde = grid%ed33 ;
  190. ims = grid%sm31 ; ime = grid%em31 ;
  191. jms = grid%sm32 ; jme = grid%em32 ;
  192. kms = grid%sm33 ; kme = grid%em33 ;
  193. its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
  194. jts = grid%sp32 ; jte = grid%ep32 ; ! tile is entire patch
  195. kts = grid%sp33 ; kte = grid%ep33 ; ! tile is entire patch
  196. CASE ( DATA_ORDER_XZY )
  197. ids = grid%sd31 ; ide = grid%ed31 ;
  198. kds = grid%sd32 ; kde = grid%ed32 ;
  199. jds = grid%sd33 ; jde = grid%ed33 ;
  200. ims = grid%sm31 ; ime = grid%em31 ;
  201. kms = grid%sm32 ; kme = grid%em32 ;
  202. jms = grid%sm33 ; jme = grid%em33 ;
  203. its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
  204. kts = grid%sp32 ; kte = grid%ep32 ; ! tile is entire patch
  205. jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
  206. END SELECT
  207. #ifdef DM_PARALLEL
  208. CALL WRF_GET_MYPROC(MYPE)
  209. CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP)
  210. call wrf_get_nprocx(inpes)
  211. call wrf_get_nprocy(jnpes)
  212. !
  213. allocate(itemp(inpes,jnpes),stat=istat)
  214. npe=0
  215. !
  216. do j=1,jnpes
  217. do i=1,inpes
  218. itemp(i,j)=npe
  219. if(npe==mype)then
  220. myi=i
  221. myj=j
  222. endif
  223. npe=npe+1
  224. enddo
  225. enddo
  226. !
  227. my_n=-1
  228. if(myj+1<=jnpes)my_n=itemp(myi,myj+1)
  229. !
  230. my_e=-1
  231. if(myi+1<=inpes)my_e=itemp(myi+1,myj)
  232. !
  233. my_s=-1
  234. if(myj-1>=1)my_s=itemp(myi,myj-1)
  235. !
  236. my_w=-1
  237. if(myi-1>=1)my_w=itemp(myi-1,myj)
  238. !
  239. my_ne=-1
  240. if((myi+1<=inpes).and.(myj+1<=jnpes)) &
  241. my_ne=itemp(myi+1,myj+1)
  242. !
  243. my_se=-1
  244. if((myi+1<=inpes).and.(myj-1>=1)) &
  245. my_se=itemp(myi+1,myj-1)
  246. !
  247. my_sw=-1
  248. if((myi-1>=1).and.(myj-1>=1)) &
  249. my_sw=itemp(myi-1,myj-1)
  250. !
  251. my_nw=-1
  252. if((myi-1>=1).and.(myj+1<=jnpes)) &
  253. my_nw=itemp(myi-1,myj+1)
  254. !
  255. ! my_neb(1)=my_n
  256. ! my_neb(2)=my_e
  257. ! my_neb(3)=my_s
  258. ! my_neb(4)=my_w
  259. ! my_neb(5)=my_ne
  260. ! my_neb(6)=my_se
  261. ! my_neb(7)=my_sw
  262. ! my_neb(8)=my_nw
  263. !
  264. deallocate(itemp)
  265. #endif
  266. grid%DT=float(grid%TIME_STEP)
  267. NNXP=min(ITE,IDE-1)
  268. NNYP=min(JTE,JDE-1)
  269. write(message,*) 'IDE, JDE: ', IDE,JDE
  270. write(message,*) 'NNXP, NNYP: ', NNXP,NNYP
  271. CALL wrf_message(message)
  272. JAM=6+2*(JDE-JDS-10)
  273. if (internal_time_loop .eq. 1) then
  274. ALLOCATE(ADUM2D(grid%sm31:grid%em31,jms:jme))
  275. ALLOCATE(KHL2(JTS:NNYP),KVL2(JTS:NNYP),KHH2(JTS:NNYP),KVH2(JTS:NNYP))
  276. ALLOCATE(DXJ(JTS:NNYP),WPDARJ(JTS:NNYP),CPGFUJ(JTS:NNYP),CURVJ(JTS:NNYP))
  277. ALLOCATE(FCPJ(JTS:NNYP),FDIVJ(JTS:NNYP),&
  278. FADJ(JTS:NNYP))
  279. ALLOCATE(HDACJ(JTS:NNYP),DDMPUJ(JTS:NNYP),DDMPVJ(JTS:NNYP))
  280. ALLOCATE(KHLA(JAM),KHHA(JAM))
  281. ALLOCATE(KVLA(JAM),KVHA(JAM))
  282. endif
  283. CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
  284. IF ( CONFIG_FLAGS%FRACTIONAL_SEAICE == 1 ) THEN
  285. CALL WRF_ERROR_FATAL("NMM cannot use FRACTIONAL_SEAICE = 1")
  286. ENDIF
  287. if ( config_flags%bl_pbl_physics == BOULACSCHEME ) then
  288. call wrf_error_fatal("Cannot use BOULAC PBL with NMM")
  289. endif
  290. write(message,*) 'cen_lat: ', config_flags%cen_lat
  291. CALL wrf_debug(100,message)
  292. write(message,*) 'cen_lon: ', config_flags%cen_lon
  293. CALL wrf_debug(100,message)
  294. write(message,*) 'dx: ', config_flags%dx
  295. CALL wrf_debug(100,message)
  296. write(message,*) 'dy: ', config_flags%dy
  297. CALL wrf_debug(100,message)
  298. write(message,*) 'config_flags%start_year: ', config_flags%start_year
  299. CALL wrf_debug(100,message)
  300. write(message,*) 'config_flags%start_month: ', config_flags%start_month
  301. CALL wrf_debug(100,message)
  302. write(message,*) 'config_flags%start_day: ', config_flags%start_day
  303. CALL wrf_debug(100,message)
  304. write(message,*) 'config_flags%start_hour: ', config_flags%start_hour
  305. CALL wrf_debug(100,message)
  306. write(start_date,435) config_flags%start_year, config_flags%start_month, &
  307. config_flags%start_day, config_flags%start_hour
  308. 435 format(I4,'-',I2.2,'-',I2.2,'_',I2.2,':00:00')
  309. grid%dlmd=config_flags%dx
  310. grid%dphd=config_flags%dy
  311. tph0d=config_flags%cen_lat
  312. tlm0d=config_flags%cen_lon
  313. !==========================================================================
  314. !!
  315. ! Check to see if the boundary conditions are set
  316. ! properly in the namelist file.
  317. ! This checks for sufficiency and redundancy.
  318. CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
  319. ! Some sort of "this is the first time" initialization. Who knows.
  320. grid%itimestep=0
  321. ! Pull in the info in the namelist to compare it to the input data.
  322. grid%real_data_init_type = model_config_rec%real_data_init_type
  323. write(message,*) 'what is flag_metgrid: ', flag_metgrid
  324. CALL wrf_message(message)
  325. IF ( flag_metgrid .EQ. 1 ) THEN ! <----- START OF VERTICAL INTERPOLATION PART ---->
  326. num_metgrid_levels = grid%num_metgrid_levels
  327. !---------------------------------------------------------------------
  328. !
  329. ! gopal's doing for ideal cases
  330. !
  331. num_veg_gc = SIZE ( grid%landusef_gc , DIM=3 )
  332. num_soil_top_gc = SIZE ( grid%soilctop_gc , DIM=3 )
  333. num_soil_bot_gc = SIZE ( grid%soilcbot_gc , DIM=3 )
  334. #ifdef DM_PARALLEL
  335. ips=its ; ipe=ite ; jps=jts ; jpe=jte ; kps=kts ; kpe=kte
  336. !JWB # include "HALO_NMM_IDEAL_1.inc"
  337. #endif
  338. !
  339. CALL vortex ( grid%ght_gc,grid%rh_gc,grid%t_gc,grid%u_gc,grid%v_gc,grid%p_gc &
  340. &, ght_out,rh_out,t_out,u_out,v_out &
  341. &, grid%ht_gc,grid%tsk_gc,grid%xice_gc &
  342. &, grid%hlat_gc,grid%hlon_gc,grid%vlat_gc,grid%vlon_gc &
  343. &, grid%greenfrac_gc,grid%albedo12m_gc,grid%landusef_gc &
  344. &, grid%soilctop_gc,grid%soilcbot_gc &
  345. &, landusef_out,soilctop_out,soilcbot_out &
  346. &, num_veg_gc,num_soil_top_gc,num_soil_bot_gc &
  347. &, config_flags%dx,internal_time_loop &
  348. &, 1,grid%num_metgrid_levels &
  349. &, ids,ide,jds,jde,kds,kde &
  350. &, ims,ime,jms,jme,kms,kme &
  351. &, its,ite,jts,jte,kts,kte )
  352. !----------------------------------------------------------------------
  353. IF (grid%ght_gc(its,jts,num_metgrid_levels/2) .lt. grid%ght_gc(its,jts,num_metgrid_levels/2+1)) THEN
  354. write(message,*) 'normal ground up file order'
  355. hyb_coor=.false.
  356. CALL wrf_message(message)
  357. ELSE
  358. hyb_coor=.true.
  359. write(message,*) 'reverse the order of coordinate'
  360. CALL wrf_message(message)
  361. CALL reverse_vert_coord(grid%ght_gc, 2, num_metgrid_levels &
  362. &, IDS,IDE,JDS,JDE,KDS,KDE &
  363. &, IMS,IME,JMS,JME,KMS,KME &
  364. &, ITS,ITE,JTS,JTE,KTS,KTE )
  365. #if defined(HWRF)
  366. if(.not. grid%use_prep_hybrid) then
  367. #endif
  368. CALL reverse_vert_coord(grid%p_gc, 2, num_metgrid_levels &
  369. &, IDS,IDE,JDS,JDE,KDS,KDE &
  370. &, IMS,IME,JMS,JME,KMS,KME &
  371. &, ITS,ITE,JTS,JTE,KTS,KTE )
  372. CALL reverse_vert_coord(grid%t_gc, 2, num_metgrid_levels &
  373. &, IDS,IDE,JDS,JDE,KDS,KDE &
  374. &, IMS,IME,JMS,JME,KMS,KME &
  375. &, ITS,ITE,JTS,JTE,KTS,KTE )
  376. CALL reverse_vert_coord(grid%u_gc, 2, num_metgrid_levels &
  377. &, IDS,IDE,JDS,JDE,KDS,KDE &
  378. &, IMS,IME,JMS,JME,KMS,KME &
  379. &, ITS,ITE,JTS,JTE,KTS,KTE )
  380. CALL reverse_vert_coord(grid%v_gc, 2, num_metgrid_levels &
  381. &, IDS,IDE,JDS,JDE,KDS,KDE &
  382. &, IMS,IME,JMS,JME,KMS,KME &
  383. &, ITS,ITE,JTS,JTE,KTS,KTE )
  384. CALL reverse_vert_coord(grid%rh_gc, 2, num_metgrid_levels &
  385. &, IDS,IDE,JDS,JDE,KDS,KDE &
  386. &, IMS,IME,JMS,JME,KMS,KME &
  387. &, ITS,ITE,JTS,JTE,KTS,KTE )
  388. #if defined(HWRF)
  389. endif
  390. #endif
  391. endif
  392. IF (hyb_coor) THEN
  393. ! limit extreme deviations from source model topography
  394. ! due to potential for nasty extrapolation/interpolation issues
  395. !
  396. write(message,*) 'min, max of grid%ht_gc before adjust: ', minval(grid%ht_gc), maxval(grid%ht_gc)
  397. CALL wrf_debug(100,message)
  398. ICOUNT=0
  399. DO J=JTS,min(JTE,JDE-1)
  400. DO I=ITS,min(ITE,IDE-1)
  401. IF ((grid%ht_gc(I,J) - grid%ght_gc(I,J,2)) .LT. -150.) THEN
  402. grid%ht_gc(I,J)=grid%ght_gc(I,J,2)-150.
  403. IF (ICOUNT .LT. 20) THEN
  404. write(message,*) 'increasing NMM topo toward RUC ', I,J
  405. CALL wrf_debug(100,message)
  406. ICOUNT=ICOUNT+1
  407. ENDIF
  408. ELSEIF ((grid%ht_gc(I,J) - grid%ght_gc(I,J,2)) .GT. 150.) THEN
  409. grid%ht_gc(I,J)=grid%ght_gc(I,J,2)+150.
  410. IF (ICOUNT .LT. 20) THEN
  411. write(message,*) 'decreasing NMM topo toward RUC ', I,J
  412. CALL wrf_debug(100,message)
  413. ICOUNT=ICOUNT+1
  414. ENDIF
  415. ENDIF
  416. END DO
  417. END DO
  418. write(message,*) 'min, max of ht_gc after correct: ', minval(grid%ht_gc), maxval(grid%ht_gc)
  419. CALL wrf_debug(100,message)
  420. ENDIF
  421. CALL boundary_smooth(grid%ht_gc,grid%landmask, grid, 12 , 12 &
  422. &, IDS,IDE,JDS,JDE,KDS,KDE &
  423. &, IMS,IME,JMS,JME,KMS,KME &
  424. &, ITS,ITE,JTS,JTE,KTS,KTE )
  425. DO j = jts, MIN(jte,jde-1)
  426. DO i = its, MIN(ite,ide-1)
  427. if (grid%landmask(I,J) .gt. 0.5) grid%sm(I,J)=0.
  428. if (grid%landmask(I,J) .le. 0.5) grid%sm(I,J)=1.
  429. if (grid%tsk_gc(I,J) .gt. 0.) then
  430. grid%nmm_tsk(I,J)=grid%tsk_gc(I,J)
  431. else
  432. #if defined(HWRF)
  433. if(grid%use_prep_hybrid) then
  434. if(grid%t(I,J,1)<100) then
  435. write(*,*) 'NO VALID SURFACE TEMPERATURE: I,J,TSK_GC(I,J),T(I,J,1) = ', &
  436. I,J,grid%TSK_GC(I,J),grid%T(I,J,1)
  437. else
  438. grid%nmm_tsk(I,J)=grid%t(I,J,1) ! stopgap measure
  439. end if
  440. else
  441. #endif
  442. grid%nmm_tsk(I,J)=grid%t_gc(I,J,1) ! stopgap measure
  443. #if defined(HWRF)
  444. endif
  445. #endif
  446. endif
  447. !
  448. grid%glat(I,J)=grid%hlat_gc(I,J)*DEGRAD
  449. grid%glon(I,J)=grid%hlon_gc(I,J)*DEGRAD
  450. grid%weasd(I,J)=grid%snow(I,J)
  451. grid%xice(I,J)=grid%xice_gc(I,J)
  452. ENDDO
  453. ENDDO
  454. ! First item is to define the target vertical coordinate
  455. num_metgrid_levels = grid%num_metgrid_levels
  456. eta_levels(1:kde) = model_config_rec%eta_levels(1:kde)
  457. ptsgm = model_config_rec%ptsgm
  458. p_top_requested = grid%p_top_requested
  459. grid%pt=p_top_requested
  460. if (internal_time_loop .eq. 1) then
  461. if (eta_levels(1) .ne. 1.0) then
  462. #if defined(HWRF)
  463. if(grid%use_prep_hybrid) then
  464. call wrf_error_fatal('PREP_HYBRID ERROR: eta_levels is not specified, but use_prep_hybrid=.true.')
  465. end if
  466. #endif
  467. write(message,*) '********************************************************************* '
  468. CALL wrf_message(message)
  469. write(message,*) '** eta_levels appears not to be specified in the namelist'
  470. CALL wrf_message(message)
  471. write(message,*) '** We will call compute_nmm_levels to define layer thicknesses.'
  472. CALL wrf_message(message)
  473. write(message,*) '** These levels should be reasonable for running the model, '
  474. CALL wrf_message(message)
  475. write(message,*) '** but may not be ideal for the simulation being made. Consider '
  476. CALL wrf_message(message)
  477. write(message,*) '** defining your own levels by specifying eta_levels in the model '
  478. CALL wrf_message(message)
  479. write(message,*) '** namelist. '
  480. CALL wrf_message(message)
  481. write(message,*) '********************************************************************** '
  482. CALL wrf_message(message)
  483. CALL compute_nmm_levels(KDE,p_top_requested,eta_levels)
  484. DO L=1,KDE
  485. write(message,*) 'L, eta_levels(L) returned :: ', L,eta_levels(L)
  486. CALL wrf_message(message)
  487. ENDDO
  488. endif
  489. write(message,*) 'KDE-1: ', KDE-1
  490. CALL wrf_debug(1,message)
  491. allocate(SG1(1:KDE-1))
  492. allocate(SG2(1:KDE-1))
  493. allocate(DSG1(1:KDE-1))
  494. allocate(DSG2(1:KDE-1))
  495. allocate(SGML1(1:KDE))
  496. allocate(SGML2(1:KDE))
  497. CALL define_nmm_vertical_coord (kde-1, ptsgm, grid%pt,grid%pdtop, eta_levels, &
  498. grid%eta1,grid%deta1,grid%aeta1, &
  499. grid%eta2,grid%deta2,grid%aeta2, grid%dfl, grid%dfrlg )
  500. DO L=KDS,KDE-1
  501. grid%deta(L)=eta_levels(L)-eta_levels(L+1)
  502. ENDDO
  503. endif
  504. write(message,*) 'num_metgrid_levels: ', num_metgrid_levels
  505. CALL wrf_message(message)
  506. DO j = jts, MIN(jte,jde-1)
  507. DO i = its, MIN(ite,ide-1)
  508. grid%fis(I,J)=grid%ht_gc(I,J)*g
  509. !
  510. ! IF ( grid%p_gc(I,J,1) .ne. 200100. .AND. (grid%ht_gc(I,J) .eq. grid%ght_gc(I,J,1)) .AND. grid%ht_gc(I,J) .ne. 0) THEN
  511. IF ( grid%p_gc(I,J,1) .ne. 200100. .AND. (abs(grid%ht_gc(I,J)-grid%ght_gc(I,J,1)) .lt. 0.01) .AND. grid%ht_gc(I,J) .ne. 0) THEN
  512. IF (mod(I,10) .eq. 0 .and. mod(J,10) .eq. 0) THEN
  513. write(message,*) 'grid%ht_gc and grid%toposoil to swap, flag_soilhgt ::: ', &
  514. I,J, grid%ht_gc(I,J),grid%toposoil(I,J),flag_soilhgt
  515. CALL wrf_debug(10,message)
  516. ENDIF
  517. IF ( ( flag_soilhgt.EQ. 1 ) ) THEN
  518. grid%ght_gc(I,J,1)=grid%toposoil(I,J)
  519. ENDIF
  520. ENDIF
  521. ENDDO
  522. ENDDO
  523. numzero=0
  524. numexamined=0
  525. DO j = jts, MIN(jte,jde-1)
  526. DO i = its, MIN(ite,ide-1)
  527. numexamined=numexamined+1
  528. if(grid%fis(i,j)<1e-5 .and. grid%fis(i,j)>-1e5 ) then
  529. numzero=numzero+1
  530. end if
  531. enddo
  532. enddo
  533. write(message,*) 'TOTAL NEAR-ZERO FIS POINTS: ',numzero,' OF ',numexamined
  534. call wrf_debug(10,message)
  535. #if defined(HWRF)
  536. interp_notph: if(.not. grid%use_prep_hybrid) then
  537. #endif
  538. if (.NOT. allocated(PDVP)) allocate(PDVP(IMS:IME,JMS:JME))
  539. if (.NOT. allocated(P3D_OUT)) allocate(P3D_OUT(IMS:IME,JMS:JME,KDS:KDE-1))
  540. if (.NOT. allocated(PSFC_OUTV)) allocate(PSFC_OUTV(IMS:IME,JMS:JME))
  541. if (.NOT. allocated(P3DV_OUT)) allocate(P3DV_OUT(IMS:IME,JMS:JME,KDS:KDE-1))
  542. if (.NOT. allocated(P3DV_IN)) allocate(P3DV_IN(IMS:IME,JMS:JME,num_metgrid_levels))
  543. CALL compute_nmm_surfacep (grid%ht_gc, grid%ght_gc, grid%p_gc , grid%t_gc &
  544. &, grid%psfc_out, num_metgrid_levels &
  545. &, IDS,IDE,JDS,JDE,KDS,KDE &
  546. &, IMS,IME,JMS,JME,KMS,KME &
  547. &, ITS,ITE,JTS,JTE,KTS,KTE ) ! H points
  548. if (internal_time_loop .eq. 1) then
  549. write(message,*) 'psfc points (final combined)'
  550. loopinc=max( (JTE-JTS)/20,1)
  551. iloopinc=max( (ITE-ITS)/10,1)
  552. CALL wrf_message(message)
  553. DO J=min(JTE,JDE-1),JTS,-loopinc
  554. write(message,633) (grid%psfc_out(I,J)/100.,I=its,min(ite,IDE-1),iloopinc)
  555. CALL wrf_message(message)
  556. ENDDO
  557. endif
  558. 633 format(35(f5.0,1x))
  559. CALL compute_3d_pressure (grid%psfc_out,grid%aeta1,grid%aeta2 &
  560. &, grid%pdtop,grid%pt,grid%pd,p3d_out &
  561. &, IDS,IDE,JDS,JDE,KDS,KDE &
  562. &, IMS,IME,JMS,JME,KMS,KME &
  563. &, ITS,ITE,JTS,JTE,KTS,KTE )
  564. #ifdef DM_PARALLEL
  565. ips=its ; ipe=ite ; jps=jts ; jpe=jte ; kps=kts ; kpe=kte
  566. # include "HALO_NMM_MG2.inc"
  567. #endif
  568. #ifdef DM_PARALLEL
  569. # include "HALO_NMM_MG3.inc"
  570. #endif
  571. do K=1,num_metgrid_levels
  572. do J=JTS,min(JTE,JDE-1)
  573. do I=ITS,min(ITE,IDE-1)
  574. IF (K .eq. KTS) THEN
  575. IF (J .eq. JDS .and. I .lt. IDE-1) THEN ! S boundary
  576. PDVP(I,J)=0.5*(grid%pd(I,J)+grid%pd(I+1,J))
  577. PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J)+grid%psfc_out(I+1,J))
  578. ELSEIF (J .eq. JDE-1 .and. I .lt. IDE-1) THEN ! N boundary
  579. PDVP(I,J)=0.5*(grid%pd(I,J)+grid%pd(I+1,J))
  580. PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J)+grid%psfc_out(I+1,J))
  581. ELSEIF (I .eq. IDS .and. mod(J,2) .eq. 0) THEN ! W boundary
  582. PDVP(I,J)=0.5*(grid%pd(I,J-1)+grid%pd(I,J+1))
  583. PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J-1)+grid%psfc_out(I,J+1))
  584. ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 0) THEN ! E boundary
  585. PDVP(I,J)=0.5*(grid%pd(I,J-1)+grid%pd(I,J+1))
  586. PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J-1)+grid%psfc_out(I,J+1))
  587. ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 1) THEN ! phantom E boundary
  588. PDVP(I,J)=grid%pd(I,J)
  589. PSFC_OUTV(I,J)=grid%psfc_out(I,J)
  590. ELSEIF (mod(J,2) .eq. 0) THEN ! interior even row
  591. PDVP(I,J)=0.25*(grid%pd(I,J)+grid%pd(I-1,J)+grid%pd(I,J+1)+grid%pd(I,J-1))
  592. PSFC_OUTV(I,J)=0.25*(grid%psfc_out(I,J)+grid%psfc_out(I-1,J)+ &
  593. grid%psfc_out(I,J+1)+grid%psfc_out(I,J-1))
  594. ELSE ! interior odd row
  595. PDVP(I,J)=0.25*(grid%pd(I,J)+grid%pd(I+1,J)+grid%pd(I,J+1)+grid%pd(I,J-1))
  596. PSFC_OUTV(I,J)=0.25*(grid%psfc_out(I,J)+grid%psfc_out(I+1,J)+ &
  597. grid%psfc_out(I,J+1)+grid%psfc_out(I,J-1))
  598. ENDIF
  599. ENDIF
  600. IF (J .eq. JDS .and. I .lt. IDE-1) THEN ! S boundary
  601. P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J,K)+grid%p_gc(I+1,J,K))
  602. ELSEIF (J .eq. JDE-1 .and. I .lt. IDE-1) THEN ! N boundary
  603. P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J,K)+grid%p_gc(I+1,J,K))
  604. ELSEIF (I .eq. IDS .and. mod(J,2) .eq. 0) THEN ! W boundary
  605. P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J-1,K)+grid%p_gc(I,J+1,K))
  606. ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 0) THEN ! E boundary
  607. P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J-1,K)+grid%p_gc(I,J+1,K))
  608. ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 1) THEN ! phantom E boundary
  609. P3DV_IN(I,J,K)=grid%p_gc(I,J,K)
  610. ELSEIF (mod(J,2) .eq. 0) THEN ! interior even row
  611. P3DV_IN(I,J,K)=0.25*(grid%p_gc(I,J,K)+grid%p_gc(I-1,J,K) + &
  612. grid%p_gc(I,J+1,K)+grid%p_gc(I,J-1,K))
  613. ELSE ! interior odd row
  614. P3DV_IN(I,J,K)=0.25*(grid%p_gc(I,J,K)+grid%p_gc(I+1,J,K) + &
  615. grid%p_gc(I,J+1,K)+grid%p_gc(I,J-1,K))
  616. ENDIF
  617. enddo
  618. enddo
  619. enddo
  620. CALL compute_3d_pressure (psfc_outv,grid%aeta1,grid%aeta2 &
  621. &, grid%pdtop,grid%pt,pdvp,p3dv_out &
  622. &, IDS,IDE,JDS,JDE,KDS,KDE &
  623. &, IMS,IME,JMS,JME,KMS,KME &
  624. &, ITS,ITE,JTS,JTE,KTS,KTE )
  625. CALL interp_press2press_lin(grid%p_gc, p3d_out &
  626. &, grid%t_gc, grid%t,num_metgrid_levels &
  627. &, .TRUE.,.TRUE.,.TRUE. & ! extrap, ignore_lowest, t_field
  628. &, IDS,IDE,JDS,JDE,KDS,KDE &
  629. &, IMS,IME,JMS,JME,KMS,KME &
  630. &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
  631. CALL interp_press2press_lin(p3dv_in, p3dv_out &
  632. &, grid%u_gc, grid%u,num_metgrid_levels &
  633. &, .FALSE.,.TRUE.,.FALSE. &
  634. &, IDS,IDE,JDS,JDE,KDS,KDE &
  635. &, IMS,IME,JMS,JME,KMS,KME &
  636. &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
  637. CALL interp_press2press_lin(p3dv_in, p3dv_out &
  638. &, grid%v_gc, grid%v,num_metgrid_levels &
  639. &, .FALSE.,.TRUE.,.FALSE. &
  640. &, IDS,IDE,JDS,JDE,KDS,KDE &
  641. &, IMS,IME,JMS,JME,KMS,KME &
  642. &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
  643. IF (hyb_coor) THEN
  644. CALL wind_adjust(p3dv_in,p3dv_out,grid%u_gc,grid%v_gc,grid%u,grid%v &
  645. &, num_metgrid_levels,5000. &
  646. &, IDS,IDE,JDS,JDE,KDS,KDE &
  647. &, IMS,IME,JMS,JME,KMS,KME &
  648. &, ITS,ITE,JTS,JTE,KTS,KTE )
  649. ENDIF
  650. ALLOCATE(qtmp(IMS:IME,JMS:JME,num_metgrid_levels))
  651. ALLOCATE(qtmp2(IMS:IME,JMS:JME,num_metgrid_levels))
  652. CALL rh_to_mxrat (grid%rh_gc, grid%t_gc, grid%p_gc, qtmp , .TRUE. , &
  653. ids , ide , jds , jde , 1 , num_metgrid_levels , &
  654. ims , ime , jms , jme , 1 , num_metgrid_levels , &
  655. its , ite , jts , jte , 1 , num_metgrid_levels )
  656. do K=1,num_metgrid_levels
  657. do J=JTS,min(JTE,JDE-1)
  658. do I=ITS,min(ITE,IDE-1)
  659. QTMP2(I,J,K)=QTMP(I,J,K)/(1.0+QTMP(I,J,K))
  660. end do
  661. end do
  662. end do
  663. CALL interp_press2press_log(grid%p_gc, p3d_out &
  664. &, QTMP2, grid%q,num_metgrid_levels &
  665. &, .FALSE.,.TRUE. &
  666. &, IDS,IDE,JDS,JDE,KDS,KDE &
  667. &, IMS,IME,JMS,JME,KMS,KME &
  668. &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
  669. IF (ALLOCATED(QTMP)) DEALLOCATE(QTMP)
  670. IF (ALLOCATED(QTMP)) DEALLOCATE(QTMP2)
  671. #if defined(HWRF)
  672. else ! we are using prep_hybrid
  673. ! Compute surface pressure:
  674. grid%psfc_out=grid%pdtop+grid%pd
  675. end if interp_notph
  676. #endif
  677. ! Get the monthly values interpolated to the current date
  678. ! for the traditional monthly
  679. ! fields of green-ness fraction and background grid%albedo.
  680. if (internal_time_loop .eq. 1 .or. config_flags%sst_update .eq. 1) then
  681. CALL monthly_interp_to_date ( grid%greenfrac_gc , current_date , grid%vegfra , &
  682. ids , ide , jds , jde , kds , kde , &
  683. ims , ime , jms , jme , kms , kme , &
  684. its , ite , jts , jte , kts , kte )
  685. CALL monthly_interp_to_date ( grid%albedo12m_gc , current_date , grid%albbck , &
  686. ids , ide , jds , jde , kds , kde , &
  687. ims , ime , jms , jme , kms , kme , &
  688. its , ite , jts , jte , kts , kte )
  689. ! Get the min/max of each i,j for the monthly green-ness fraction.
  690. CALL monthly_min_max ( grid%greenfrac_gc , grid%shdmin , grid%shdmax , &
  691. ids , ide , jds , jde , kds , kde , &
  692. ims , ime , jms , jme , kms , kme , &
  693. its , ite , jts , jte , kts , kte )
  694. ! The model expects the green-ness values in percent, not fraction.
  695. DO j = jts, MIN(jte,jde-1)
  696. DO i = its, MIN(ite,ide-1)
  697. !! grid%vegfra(i,j) = grid%vegfra(i,j) * 100.
  698. grid%shdmax(i,j) = grid%shdmax(i,j) * 100.
  699. grid%shdmin(i,j) = grid%shdmin(i,j) * 100.
  700. grid%vegfrc(I,J)=grid%vegfra(I,J)
  701. END DO
  702. END DO
  703. ! The model expects the albedo fields as
  704. ! a fraction, not a percent. Set the water values to 8%.
  705. DO j = jts, MIN(jte,jde-1)
  706. DO i = its, MIN(ite,ide-1)
  707. if (grid%albbck(i,j) .lt. 5.) then
  708. write(message,*) 'reset albedo to 8%... I,J,albbck:: ', I,J,grid%albbck(I,J)
  709. CALL wrf_debug(10,message)
  710. grid%albbck(I,J)=8.
  711. endif
  712. grid%albbck(i,j) = grid%albbck(i,j) / 100.
  713. grid%snoalb(i,j) = grid%snoalb(i,j) / 100.
  714. IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
  715. grid%albbck(i,j) = 0.08
  716. grid%snoalb(i,j) = 0.08
  717. END IF
  718. grid%albase(i,j)=grid%albbck(i,j)
  719. grid%mxsnal(i,j)=grid%snoalb(i,j)
  720. END DO
  721. END DO
  722. endif
  723. #if defined(HWRF)
  724. if(.not.grid%use_prep_hybrid) then
  725. #endif
  726. ! new deallocs
  727. DEALLOCATE(p3d_out,p3dv_out,p3dv_in)
  728. #if defined(HWRF)
  729. end if
  730. #endif
  731. END IF ! <----- END OF VERTICAL INTERPOLATION PART ---->
  732. !! compute SST at each time if updating SST
  733. if (config_flags%sst_update == 1) then
  734. DO j = jts, MIN(jde-1,jte)
  735. DO i = its, MIN(ide-1,ite)
  736. if (grid%SM(I,J) .lt. 0.5) then
  737. grid%SST(I,J)=0.
  738. endif
  739. if (grid%SM(I,J) .gt. 0.5) then
  740. grid%SST(I,J)=grid%NMM_TSK(I,J)
  741. grid%NMM_TSK(I,J)=0.
  742. endif
  743. IF ( (grid%NMM_TSK(I,J)+grid%SST(I,J)) .lt. 200. .or. &
  744. (grid%NMM_TSK(I,J)+grid%SST(I,J)) .gt. 350. ) THEN
  745. write(message,*) 'TSK, SST trouble at : ', I,J
  746. CALL wrf_message(message)
  747. write(message,*) 'SM, NMM_TSK,SST ', grid%SM(I,J),grid%NMM_TSK(I,J),grid%SST(I,J)
  748. CALL wrf_message(message)
  749. ENDIF
  750. ENDDO
  751. ENDDO
  752. endif ! sst_update test
  753. if (internal_time_loop .eq. 1) then
  754. !!! weasd has "snow water equivalent" in mm
  755. DO j = jts, MIN(jte,jde-1)
  756. DO i = its, MIN(ite,ide-1)
  757. IF(grid%sm(I,J).GT.0.9) THEN
  758. IF (grid%xice(I,J) .gt. 0) then
  759. grid%si(I,J)=1.0
  760. ENDIF
  761. ! SEA
  762. grid%epsr(I,J)=.97
  763. grid%embck(I,J)=.97
  764. grid%gffc(I,J)=0.
  765. grid%albedo(I,J)=.06
  766. grid%albase(I,J)=.06
  767. IF(grid%si (I,J).GT.0. ) THEN
  768. ! SEA-ICE
  769. grid%sm(I,J)=0.
  770. grid%si(I,J)=0.
  771. grid%sice(I,J)=1.
  772. grid%gffc(I,J)=0. ! just leave zero as irrelevant
  773. grid%albedo(I,J)=.60
  774. grid%albase(I,J)=.60
  775. ENDIF
  776. ELSE
  777. grid%si(I,J)=5.0*grid%weasd(I,J)/1000.
  778. ! LAND
  779. grid%epsr(I,J)=1.0
  780. grid%embck(I,J)=1.0
  781. grid%gffc(I,J)=0.0 ! just leave zero as irrelevant
  782. grid%sice(I,J)=0.
  783. grid%sno(I,J)=grid%si(I,J)*.20
  784. ENDIF
  785. ENDDO
  786. ENDDO
  787. ! DETERMINE grid%albedo OVER LAND
  788. DO j = jts, MIN(jte,jde-1)
  789. DO i = its, MIN(ite,ide-1)
  790. IF(grid%sm(I,J).LT.0.9.AND.grid%sice(I,J).LT.0.9) THEN
  791. ! SNOWFREE albedo
  792. IF ( (grid%sno(I,J) .EQ. 0.0) .OR. &
  793. (grid%albase(I,J) .GE. grid%mxsnal(I,J) ) ) THEN
  794. grid%albedo(I,J) = grid%albase(I,J)
  795. ELSE
  796. ! MODIFY albedo IF SNOWCOVER:
  797. ! BELOW SNOWDEPTH THRESHOLD...
  798. IF (grid%sno(I,J) .LT. SNUP) THEN
  799. RSNOW = grid%sno(I,J)/SNUP
  800. SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP))
  801. ! ABOVE SNOWDEPTH THRESHOLD...
  802. ELSE
  803. SNOFAC = 1.0
  804. ENDIF
  805. ! CALCULATE grid%albedo ACCOUNTING FOR SNOWDEPTH AND VGFRCK
  806. grid%albedo(I,J) = grid%albase(I,J) &
  807. + (1.0-grid%vegfra(I,J))*SNOFAC*(grid%mxsnal(I,J)-grid%albase(I,J))
  808. ENDIF
  809. END IF
  810. grid%si(I,J)=5.0*grid%weasd(I,J)
  811. grid%sno(I,J)=grid%weasd(I,J)
  812. !! convert vegfra
  813. grid%vegfra(I,J)=grid%vegfra(I,J)*100.
  814. !
  815. ENDDO
  816. ENDDO
  817. #ifdef DM_PARALLEL
  818. ALLOCATE(SM_G(IDS:IDE,JDS:JDE),SICE_G(IDS:IDE,JDS:JDE))
  819. CALL WRF_PATCH_TO_GLOBAL_REAL( grid%sice(IMS,JMS) &
  820. &, SICE_G,grid%DOMDESC &
  821. &, 'z','xy' &
  822. &, IDS,IDE-1,JDS,JDE-1,1,1 &
  823. &, IMS,IME,JMS,JME,1,1 &
  824. &, ITS,ITE,JTS,JTE,1,1 )
  825. CALL WRF_PATCH_TO_GLOBAL_REAL( grid%sm(IMS,JMS) &
  826. &, SM_G,grid%DOMDESC &
  827. &, 'z','xy' &
  828. &, IDS,IDE-1,JDS,JDE-1,1,1 &
  829. &, IMS,IME,JMS,JME,1,1 &
  830. &, ITS,ITE,JTS,JTE,1,1 )
  831. IF (WRF_DM_ON_MONITOR()) THEN
  832. 637 format(40(f3.0,1x))
  833. allocate(IHE_G(JDS:JDE-1),IHW_G(JDS:JDE-1))
  834. DO j = JDS, JDE-1
  835. IHE_G(J)=MOD(J+1,2)
  836. IHW_G(J)=IHE_G(J)-1
  837. ENDDO
  838. DO ITER=1,10
  839. DO j = jds+1, (jde-1)-1
  840. DO i = ids+1, (ide-1)-1
  841. ! any sea ice around point in question?
  842. IF (SM_G(I,J) .ge. 0.9) THEN
  843. SEAICESUM=SICE_G(I+IHE_G(J),J+1)+SICE_G(I+IHW_G(J),J+1)+ &
  844. SICE_G(I+IHE_G(J),J-1)+SICE_G(I+IHW_G(J),J-1)
  845. IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN
  846. IF ((SICE_G(I+IHE_G(J),J+1).lt.0.1 .and. SM_G(I+IHE_G(J),J+1).lt.0.1) .OR. &
  847. (SICE_G(I+IHW_G(J),J+1).lt.0.1 .and. SM_G(I+IHW_G(J),J+1).lt.0.1) .OR. &
  848. (SICE_G(I+IHE_G(J),J-1).lt.0.1 .and. SM_G(I+IHE_G(J),J-1).lt.0.1) .OR. &
  849. (SICE_G(I+IHW_G(J),J-1).lt.0.1 .and. SM_G(I+IHW_G(J),J-1).lt.0.1)) THEN
  850. ! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE
  851. write(message,*) 'making seaice (1): ', I,J
  852. CALL wrf_debug(100,message)
  853. SICE_G(I,J)=1.0
  854. SM_G(I,J)=0.
  855. ENDIF
  856. ELSEIF (SEAICESUM .ge. 3) THEN
  857. ! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE
  858. write(message,*) 'making seaice (2): ', I,J
  859. CALL wrf_debug(100,message)
  860. SICE_G(I,J)=1.0
  861. SM_G(I,J)=0.
  862. ENDIF
  863. ENDIF
  864. ENDDO
  865. ENDDO
  866. ENDDO
  867. ENDIF
  868. CALL WRF_GLOBAL_TO_PATCH_REAL( SICE_G, grid%sice &
  869. &, grid%DOMDESC &
  870. &, 'z','xy' &
  871. &, IDS,IDE-1,JDS,JDE-1,1,1 &
  872. &, IMS,IME,JMS,JME,1,1 &
  873. &, ITS,ITE,JTS,JTE,1,1 )
  874. CALL WRF_GLOBAL_TO_PATCH_REAL( SM_G,grid%sm &
  875. &, grid%DOMDESC &
  876. &, 'z','xy' &
  877. &, IDS,IDE-1,JDS,JDE-1,1,1 &
  878. &, IMS,IME,JMS,JME,1,1 &
  879. &, ITS,ITE,JTS,JTE,1,1 )
  880. IF (WRF_DM_ON_MONITOR()) THEN
  881. #if defined(HWRF)
  882. ! SM_G is still needed for the high-res grid
  883. #else
  884. DEALLOCATE(SM_G)
  885. #endif
  886. deallocate(SICE_G)
  887. DEALLOCATE(IHE_G,IHW_G)
  888. ENDIF
  889. ! write(message,*) 'revised sea ice on patch'
  890. ! CALL wrf_debug(100,message)
  891. ! DO J=JTE,JTS,-(((JTE-JTS)/25)+1)
  892. ! write(message,637) (grid%sice(I,J),I=ITS,ITE,ITE/20)
  893. ! CALL wrf_debug(100,message)
  894. ! END DO
  895. #else
  896. ! serial sea ice reprocessing
  897. allocate(IHE(JDS:JDE-1),IHW(JDS:JDE-1))
  898. DO j = jts, MIN(jte,jde-1)
  899. IHE(J)=MOD(J+1,2)
  900. IHW(J)=IHE(J)-1
  901. ENDDO
  902. DO ITER=1,10
  903. DO j = jts+1, MIN(jte,jde-1)-1
  904. DO i = its+1, MIN(ite,ide-1)-1
  905. ! any sea ice around point in question?
  906. IF (grid%sm(I,J) .gt. 0.9) THEN
  907. SEAICESUM=grid%sice(I+IHE(J),J+1)+grid%sice(I+IHW(J),J+1)+ &
  908. grid%sice(I+IHE(J),J-1)+grid%sice(I+IHW(J),J-1)
  909. IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN
  910. IF ((grid%sice(I+IHE(J),J+1).lt.0.1 .and. grid%sm(I+IHE(J),J+1).lt.0.1) .OR. &
  911. (grid%sice(I+IHW(J),J+1).lt.0.1 .and. grid%sm(I+IHW(J),J+1).lt.0.1) .OR. &
  912. (grid%sice(I+IHE(J),J-1).lt.0.1 .and. grid%sm(I+IHE(J),J-1).lt.0.1) .OR. &
  913. (grid%sice(I+IHW(J),J-1).lt.0.1 .and. grid%sm(I+IHW(J),J-1).lt.0.1)) THEN
  914. ! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE
  915. grid%sice(I,J)=1.0
  916. grid%sm(I,J)=0.
  917. ENDIF
  918. ELSEIF (SEAICESUM .ge. 3) THEN
  919. ! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE
  920. grid%sice(I,J)=1.0
  921. grid%sm(I,J)=0.
  922. ENDIF
  923. ENDIF
  924. ENDDO
  925. ENDDO
  926. ENDDO
  927. DEALLOCATE(IHE,IHW)
  928. #endif
  929. ! this block meant to guarantee land/sea agreement between sm and landmask
  930. DO j = jts, MIN(jte,jde-1)
  931. DO i = its, MIN(ite,ide-1)
  932. IF (grid%sm(I,J) .gt. 0.5) THEN
  933. grid%landmask(I,J)=0.0
  934. ELSEIF (grid%sm(I,J) .lt. 0.5 .and. grid%sice(I,J) .gt. 0.9) then
  935. grid%landmask(I,J)=0.0
  936. ELSEIF (grid%sm(I,J) .lt. 0.5 .and. grid%sice(I,J) .lt. 0.1) then
  937. grid%landmask(I,J)=1.0
  938. ELSE
  939. write(message,*) 'missed point in grid%landmask definition ' , I,J
  940. CALL wrf_message(message)
  941. grid%landmask(I,J)=0.0
  942. ENDIF
  943. !
  944. IF (grid%sice(I,J) .gt. 0.5 .and. grid%nmm_tsk(I,J) .lt. 0.1 .and. grid%sst(I,J) .gt. 0.) THEN
  945. write(message,*) 'set grid%nmm_tsk to: ', grid%sst(I,J)
  946. CALL wrf_message(message)
  947. grid%nmm_tsk(I,J)=grid%sst(I,J)
  948. grid%sst(I,J)=0.
  949. endif
  950. ENDDO
  951. ENDDO
  952. ! For sf_surface_physics = 1, we want to use close to a 10 cm value
  953. ! for the bottom level of the soil temps.
  954. IF ( ( model_config_rec%sf_surface_physics(grid%id) .EQ. 1 ) .AND. &
  955. ( flag_st000010 .EQ. 1 ) ) THEN
  956. DO j = jts , MIN(jde-1,jte)
  957. DO i = its , MIN(ide-1,ite)
  958. grid%soiltb(i,j) = grid%st000010(i,j)
  959. END DO
  960. END DO
  961. END IF
  962. ! Adjust the various soil temperature values depending on the difference in
  963. ! in elevation between the current model's elevation and the incoming data's
  964. ! orography.
  965. IF ( ( flag_toposoil .EQ. 1 ) ) THEN
  966. ALLOCATE(HT(ims:ime,jms:jme))
  967. DO J=jms,jme
  968. DO I=ims,ime
  969. HT(I,J)=grid%fis(I,J)/9.81
  970. END DO
  971. END DO
  972. ! if (maxval(grid%toposoil) .gt. 100.) then
  973. !
  974. ! Being avoided. Something to revisit eventually.
  975. !
  976. !1219 might be simply a matter of including toposoil
  977. !
  978. ! CODE NOT TESTED AT NCEP USING THIS FUNCTIONALITY,
  979. ! SO TO BE SAFE WILL AVOID FOR RETRO RUNS.
  980. !
  981. ! CALL adjust_soil_temp_new ( grid%soiltb , 2 , &
  982. ! grid%nmm_tsk , ht , grid%toposoil , grid%landmask, flag_toposoil , &
  983. ! grid%st000010 , st010040 , st040100 , st100200 , st010200 , &
  984. ! flag_st000010 , flag_st010040 , flag_st040100 , &
  985. ! flag_st100200 , flag_st010200 , &
  986. ! soilt000 , soilt005 , soilt020 , soilt040 , &
  987. ! soilt160 , soilt300 , &
  988. ! flag_soilt000 , flag_soilt005 , flag_soilt020 , &
  989. ! flag_soilt040 , flag_soilt160 , flag_soilt300 , &
  990. ! ids , ide , jds , jde , kds , kde , &
  991. ! ims , ime , jms , jme , kms , kme , &
  992. ! its , ite , jts , jte , kts , kte )
  993. ! endif
  994. END IF
  995. ! Process the LSM data.
  996. ! surface_input_source=1 => use data from static file
  997. ! (fractional category as input)
  998. ! surface_input_source=2 => use data from grib file
  999. ! (dominant category as input)
  1000. IF ( config_flags%surface_input_source .EQ. 1 ) THEN
  1001. grid%vegcat (its,jts) = 0
  1002. grid%soilcat(its,jts) = 0
  1003. END IF
  1004. ! Generate the vegetation and soil category information
  1005. ! from the fractional input
  1006. ! data, or use the existing dominant category fields if they exist.
  1007. IF ((grid%soilcat(its,jts) .LT. 0.5) .AND. (grid%vegcat(its,jts) .LT. 0.5)) THEN
  1008. num_veg_cat = SIZE ( grid%landusef_gc , DIM=3 )
  1009. num_soil_top_cat = SIZE ( grid%soilctop_gc , DIM=3 )
  1010. num_soil_bot_cat = SIZE ( grid%soilcbot_gc , DIM=3 )
  1011. do J=JMS,JME
  1012. do K=1,num_veg_cat
  1013. do I=IMS,IME
  1014. grid%landusef(I,K,J)=grid%landusef_gc(I,J,K)
  1015. enddo
  1016. enddo
  1017. enddo
  1018. do J=JMS,JME
  1019. do K=1,num_soil_top_cat
  1020. do I=IMS,IME
  1021. grid%soilctop(I,K,J)=grid%soilctop_gc(I,J,K)
  1022. enddo
  1023. enddo
  1024. enddo
  1025. do J=JMS,JME
  1026. do K=1,num_soil_bot_cat
  1027. do I=IMS,IME
  1028. grid%soilcbot(I,K,J)=grid%soilcbot_gc(I,J,K)
  1029. enddo
  1030. enddo
  1031. enddo
  1032. ! grid%sm (1=water, 0=land)
  1033. ! grid%landmask(0=water, 1=land)
  1034. write(message,*) 'landmask into process_percent_cat_new'
  1035. CALL wrf_debug(1,message)
  1036. do J=JTE,JTS,-(((JTE-JTS)/20)+1)
  1037. write(message,641) (grid%landmask(I,J),I=ITS,min(ITE,IDE-1),((ITE-ITS)/15)+1)
  1038. CALL wrf_debug(1,message)
  1039. enddo
  1040. 641 format(25(f3.0,1x))
  1041. CALL process_percent_cat_new ( grid%landmask , &
  1042. grid%landusef , grid%soilctop , grid%soilcbot , &
  1043. grid%isltyp , grid%ivgtyp , &
  1044. num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
  1045. ids , ide , jds , jde , kds , kde , &
  1046. ims , ime , jms , jme , kms , kme , &
  1047. its , ite , jts , jte , kts , kte , &
  1048. model_config_rec%iswater(grid%id) )
  1049. DO j = jts , MIN(jde-1,jte)
  1050. DO i = its , MIN(ide-1,ite)
  1051. grid%vegcat(i,j) = grid%ivgtyp(i,j)
  1052. grid%soilcat(i,j) = grid%isltyp(i,j)
  1053. END DO
  1054. END DO
  1055. ELSE
  1056. ! Do we have dominant soil and veg data from the input already?
  1057. IF ( grid%soilcat(its,jts) .GT. 0.5 ) THEN
  1058. DO j = jts, MIN(jde-1,jte)
  1059. DO i = its, MIN(ide-1,ite)
  1060. grid%isltyp(i,j) = NINT( grid%soilcat(i,j) )
  1061. END DO
  1062. END DO
  1063. END IF
  1064. IF ( grid%vegcat(its,jts) .GT. 0.5 ) THEN
  1065. DO j = jts, MIN(jde-1,jte)
  1066. DO i = its, MIN(ide-1,ite)
  1067. grid%ivgtyp(i,j) = NINT( grid%vegcat(i,j) )
  1068. END DO
  1069. END DO
  1070. END IF
  1071. ENDIF
  1072. DO j = jts, MIN(jde-1,jte)
  1073. DO i = its, MIN(ide-1,ite)
  1074. IF (grid%sice(I,J) .lt. 0.1) THEN
  1075. IF (grid%landmask(I,J) .gt. 0.5 .and. grid%sm(I,J) .gt. 0.5) THEN
  1076. write(message,*) 'land mask and grid%sm both > 0.5: ', &
  1077. I,J,grid%landmask(I,J),grid%sm(I,J)
  1078. CALL wrf_message(message)
  1079. grid%sm(I,J)=0.
  1080. ELSEIF (grid%landmask(I,J) .lt. 0.5 .and. grid%sm(I,J) .lt. 0.5) THEN
  1081. write(message,*) 'land mask and grid%sm both < 0.5: ', &
  1082. I,J, grid%landmask(I,J),grid%sm(I,J)
  1083. CALL wrf_message(message)
  1084. grid%sm(I,J)=1.
  1085. ENDIF
  1086. ELSE
  1087. IF (grid%landmask(I,J) .gt. 0.5 .and. grid%sm(I,J)+grid%sice(I,J) .gt. 0.9) then
  1088. write(message,*) 'landmask says LAND, sm/sice say SEAICE: ', I,J
  1089. ENDIF
  1090. ENDIF
  1091. ENDDO
  1092. ENDDO
  1093. DO j = jts, MIN(jde-1,jte)
  1094. DO i = its, MIN(ide-1,ite)
  1095. if (grid%sice(I,J) .gt. 0.9) then
  1096. grid%isltyp(I,J)=16
  1097. grid%ivgtyp(I,J)=24
  1098. endif
  1099. ENDDO
  1100. ENDDO
  1101. DO j = jts, MIN(jde-1,jte)
  1102. DO i = its, MIN(ide-1,ite)
  1103. if (grid%sm(I,J) .lt. 0.5) then
  1104. grid%sst(I,J)=0.
  1105. endif
  1106. if (grid%sm(I,J) .gt. 0.5) then
  1107. if (grid%sst(I,J) .lt. 0.1) then
  1108. grid%sst(I,J)=grid%nmm_tsk(I,J)
  1109. endif
  1110. grid%nmm_tsk(I,J)=0.
  1111. endif
  1112. IF ( (grid%nmm_tsk(I,J)+grid%sst(I,J)) .lt. 200. .or. &
  1113. (grid%nmm_tsk(I,J)+grid%sst(I,J)) .gt. 350. ) THEN
  1114. write(message,*) 'TSK, sst trouble at : ', I,J
  1115. CALL wrf_message(message)
  1116. write(message,*) 'sm, nmm_tsk,sst ', grid%sm(I,J),grid%nmm_tsk(I,J),grid%sst(I,J)
  1117. CALL wrf_message(message)

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