PageRenderTime 86ms CodeModel.GetById 22ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/dyn_nmm/module_initialize_real.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 5238 lines | 3350 code | 960 blank | 928 comment | 103 complexity | 132c80bdb132b364e3d5f3f4a821dd15 MD5 | raw file
Possible License(s): AGPL-1.0
  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_real
  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. TYPE (grid_config_rec_type) :: config_flags
  55. ! Local domain indices and counters.
  56. INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
  57. INTEGER :: &
  58. ids, ide, jds, jde, kds, kde, &
  59. ims, ime, jms, jme, kms, kme, &
  60. its, ite, jts, jte, kts, kte, &
  61. ips, ipe, jps, jpe, kps, kpe, &
  62. i, j, k, NNXP, NNYP
  63. ! Local data
  64. CHARACTER(LEN=19):: start_date
  65. #ifdef DM_PARALLEL
  66. LOGICAL,EXTERNAL :: WRF_DM_ON_MONITOR
  67. logical :: test
  68. ! INTEGER :: DOMDESC
  69. REAL,ALLOCATABLE :: SICE_G(:,:), SM_G(:,:)
  70. INTEGER, ALLOCATABLE:: IHE_G(:),IHW_G(:)
  71. INTEGER, ALLOCATABLE:: ITEMP(:,:)
  72. INTEGER :: my_e,my_n,my_s,my_w,my_ne,my_nw,my_se,my_sw,myi,myj,npe
  73. INTEGER :: istat,inpes,jnpes
  74. #else
  75. integer, allocatable:: ihw(:),ihe(:)
  76. #endif
  77. CHARACTER (LEN=255) :: message
  78. INTEGER :: error
  79. REAL :: p_surf, p_level
  80. REAL :: cof1, cof2
  81. REAL :: qvf , qvf1 , qvf2 , pd_surf
  82. REAL :: p00 , t00 , a
  83. REAL :: hold_znw, rmin,rmax
  84. REAL :: p_top_requested , ptsgm
  85. INTEGER :: num_metgrid_levels, ICOUNT
  86. REAL , DIMENSION(max_eta) :: eta_levels
  87. LOGICAL :: stretch_grid, dry_sounding, debug, log_flag_sst, hyb_coor
  88. REAL, ALLOCATABLE,DIMENSION(:,:):: ADUM2D,SNOWC,HT,TG_ALT, &
  89. PDVP,PSFC_OUTV
  90. REAL, ALLOCATABLE,DIMENSION(:,:,:):: P3D_OUT,P3DV_OUT,P3DV_IN, &
  91. QTMP,QTMP2
  92. INTEGER, ALLOCATABLE, DIMENSION(:):: KHL2,KVL2,KHH2,KVH2, &
  93. KHLA,KHHA,KVLA,KVHA
  94. ! INTEGER, ALLOCATABLE, DIMENSION(:,:):: grid%lu_index
  95. REAL, ALLOCATABLE, DIMENSION(:):: DXJ,WPDARJ,CPGFUJ,CURVJ, &
  96. FCPJ,FDIVJ,EMJ,EMTJ,FADJ, &
  97. HDACJ,DDMPUJ,DDMPVJ
  98. REAL, ALLOCATABLE,DIMENSION(:),SAVE:: SG1,SG2,DSG1,DSG2, &
  99. SGML1,SGML2
  100. !-- Carsel and Parrish [1988]
  101. REAL , DIMENSION(100) :: lqmi
  102. integer iicount
  103. REAL:: TPH0D,TLM0D
  104. REAL:: TPH0,WB,SB,TDLM,TDPH
  105. REAL:: WBI,SBI,EBI,ANBI,STPH0,CTPH0
  106. REAL:: TSPH,DTAD,DTCF
  107. REAL:: ACDT,CDDAMP,DXP,FP
  108. REAL:: WBD,SBD
  109. REAL:: RSNOW,SNOFAC
  110. REAL, PARAMETER:: SALP=2.60
  111. REAL, PARAMETER:: SNUP=0.040
  112. REAL:: SMCSUM,STCSUM,SEAICESUM,FISX
  113. REAL:: cur_smc, aposs_smc
  114. REAL:: COAC, CODAMP
  115. INTEGER,PARAMETER:: DOUBLE=SELECTED_REAL_KIND(15,300)
  116. REAL(KIND=DOUBLE):: TERM1,APH,TLM,TPH,DLM,DPH,STPH,CTPH
  117. INTEGER:: KHH,KVH,JAM,JA, IHL, IHH, L
  118. INTEGER:: II,JJ,ISRCH,ISUM,ITER,Ilook,Jlook
  119. REAL, PARAMETER:: DTR=0.01745329
  120. REAL, PARAMETER:: W_NMM=0.08
  121. #if defined(HWRF)
  122. REAL, PARAMETER:: DDFC=1.0
  123. #else
  124. REAL, PARAMETER:: DDFC=8.0
  125. #endif
  126. REAL, PARAMETER:: TWOM=.00014584
  127. REAL, PARAMETER:: CP=1004.6
  128. REAL, PARAMETER:: DFC=1.0
  129. REAL, PARAMETER:: ROI=916.6
  130. REAL, PARAMETER:: R=287.04
  131. REAL, PARAMETER:: CI=2060.0
  132. REAL, PARAMETER:: ROS=1500.
  133. REAL, PARAMETER:: CS=1339.2
  134. REAL, PARAMETER:: DS=0.050
  135. REAL, PARAMETER:: AKS=.0000005
  136. REAL, PARAMETER:: DZG=2.85
  137. REAL, PARAMETER:: DI=.1000
  138. REAL, PARAMETER:: AKI=0.000001075
  139. REAL, PARAMETER:: DZI=2.0
  140. REAL, PARAMETER:: THL=210.
  141. REAL, PARAMETER:: PLQ=70000.
  142. REAL, PARAMETER:: ERAD=6371200.
  143. REAL, PARAMETER:: TG0=258.16
  144. REAL, PARAMETER:: TGA=30.0
  145. integer :: numzero,numexamined
  146. #ifdef HWRF
  147. !============================================================================
  148. ! gopal's doing for ocean coupling
  149. !============================================================================
  150. REAL, DIMENSION(:,:), ALLOCATABLE :: NHLAT,NHLON,NVLAT,NVLON,HRES_SM
  151. REAL :: NDLMD,NDPHD,NWBD,NSBD
  152. INTEGER :: NIDE,NJDE,ILOC,JLOC
  153. INTEGER fid, ierr, nprocs
  154. CHARACTER*255 f65name, SysString
  155. !============================================================================
  156. ! end gopal's doing for ocean coupling
  157. !============================================================================
  158. #endif
  159. if (ALLOCATED(ADUM2D)) DEALLOCATE(ADUM2D)
  160. if (ALLOCATED(TG_ALT)) DEALLOCATE(TG_ALT)
  161. !#define COPY_IN
  162. !#include <scalar_derefs.inc>
  163. #ifdef DM_PARALLEL
  164. # include <data_calls.inc>
  165. #endif
  166. SELECT CASE ( model_data_order )
  167. CASE ( DATA_ORDER_ZXY )
  168. kds = grid%sd31 ; kde = grid%ed31 ;
  169. ids = grid%sd32 ; ide = grid%ed32 ;
  170. jds = grid%sd33 ; jde = grid%ed33 ;
  171. kms = grid%sm31 ; kme = grid%em31 ;
  172. ims = grid%sm32 ; ime = grid%em32 ;
  173. jms = grid%sm33 ; jme = grid%em33 ;
  174. kts = grid%sp31 ; kte = grid%ep31 ; ! tile is entire patch
  175. its = grid%sp32 ; ite = grid%ep32 ; ! tile is entire patch
  176. jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
  177. CASE ( DATA_ORDER_XYZ )
  178. ids = grid%sd31 ; ide = grid%ed31 ;
  179. jds = grid%sd32 ; jde = grid%ed32 ;
  180. kds = grid%sd33 ; kde = grid%ed33 ;
  181. ims = grid%sm31 ; ime = grid%em31 ;
  182. jms = grid%sm32 ; jme = grid%em32 ;
  183. kms = grid%sm33 ; kme = grid%em33 ;
  184. its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
  185. jts = grid%sp32 ; jte = grid%ep32 ; ! tile is entire patch
  186. kts = grid%sp33 ; kte = grid%ep33 ; ! tile is entire patch
  187. CASE ( DATA_ORDER_XZY )
  188. ids = grid%sd31 ; ide = grid%ed31 ;
  189. kds = grid%sd32 ; kde = grid%ed32 ;
  190. jds = grid%sd33 ; jde = grid%ed33 ;
  191. ims = grid%sm31 ; ime = grid%em31 ;
  192. kms = grid%sm32 ; kme = grid%em32 ;
  193. jms = grid%sm33 ; jme = grid%em33 ;
  194. its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
  195. kts = grid%sp32 ; kte = grid%ep32 ; ! tile is entire patch
  196. jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
  197. END SELECT
  198. #ifdef DM_PARALLEL
  199. CALL WRF_GET_MYPROC(MYPE)
  200. CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP)
  201. call wrf_get_nprocx(inpes)
  202. call wrf_get_nprocy(jnpes)
  203. !
  204. allocate(itemp(inpes,jnpes),stat=istat)
  205. npe=0
  206. !
  207. do j=1,jnpes
  208. do i=1,inpes
  209. itemp(i,j)=npe
  210. if(npe==mype)then
  211. myi=i
  212. myj=j
  213. endif
  214. npe=npe+1
  215. enddo
  216. enddo
  217. !
  218. my_n=-1
  219. if(myj+1<=jnpes)my_n=itemp(myi,myj+1)
  220. !
  221. my_e=-1
  222. if(myi+1<=inpes)my_e=itemp(myi+1,myj)
  223. !
  224. my_s=-1
  225. if(myj-1>=1)my_s=itemp(myi,myj-1)
  226. !
  227. my_w=-1
  228. if(myi-1>=1)my_w=itemp(myi-1,myj)
  229. !
  230. my_ne=-1
  231. if((myi+1<=inpes).and.(myj+1<=jnpes)) &
  232. my_ne=itemp(myi+1,myj+1)
  233. !
  234. my_se=-1
  235. if((myi+1<=inpes).and.(myj-1>=1)) &
  236. my_se=itemp(myi+1,myj-1)
  237. !
  238. my_sw=-1
  239. if((myi-1>=1).and.(myj-1>=1)) &
  240. my_sw=itemp(myi-1,myj-1)
  241. !
  242. my_nw=-1
  243. if((myi-1>=1).and.(myj+1<=jnpes)) &
  244. my_nw=itemp(myi-1,myj+1)
  245. !
  246. ! my_neb(1)=my_n
  247. ! my_neb(2)=my_e
  248. ! my_neb(3)=my_s
  249. ! my_neb(4)=my_w
  250. ! my_neb(5)=my_ne
  251. ! my_neb(6)=my_se
  252. ! my_neb(7)=my_sw
  253. ! my_neb(8)=my_nw
  254. !
  255. deallocate(itemp)
  256. #endif
  257. NNXP=min(ITE,IDE-1)
  258. NNYP=min(JTE,JDE-1)
  259. write(message,*) 'IDE, JDE: ', IDE,JDE
  260. write(message,*) 'NNXP, NNYP: ', NNXP,NNYP
  261. CALL wrf_message(message)
  262. JAM=6+2*(JDE-JDS-10)
  263. if (internal_time_loop .eq. 1) then
  264. ALLOCATE(ADUM2D(grid%sm31:grid%em31,jms:jme))
  265. ALLOCATE(KHL2(JTS:NNYP),KVL2(JTS:NNYP),KHH2(JTS:NNYP),KVH2(JTS:NNYP))
  266. ALLOCATE(DXJ(JTS:NNYP),WPDARJ(JTS:NNYP),CPGFUJ(JTS:NNYP),CURVJ(JTS:NNYP))
  267. ALLOCATE(FCPJ(JTS:NNYP),FDIVJ(JTS:NNYP),&
  268. FADJ(JTS:NNYP))
  269. ALLOCATE(HDACJ(JTS:NNYP),DDMPUJ(JTS:NNYP),DDMPVJ(JTS:NNYP))
  270. ALLOCATE(KHLA(JAM),KHHA(JAM))
  271. ALLOCATE(KVLA(JAM),KVHA(JAM))
  272. endif
  273. CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
  274. IF ( CONFIG_FLAGS%FRACTIONAL_SEAICE == 1 ) THEN
  275. CALL WRF_ERROR_FATAL("NMM cannot use FRACTIONAL_SEAICE = 1")
  276. ENDIF
  277. if ( config_flags%bl_pbl_physics == BOULACSCHEME ) then
  278. call wrf_error_fatal("Cannot use BOULAC PBL with NMM")
  279. endif
  280. write(message,*) 'cen_lat: ', config_flags%cen_lat
  281. CALL wrf_debug(100,message)
  282. write(message,*) 'cen_lon: ', config_flags%cen_lon
  283. CALL wrf_debug(100,message)
  284. write(message,*) 'dx: ', config_flags%dx
  285. CALL wrf_debug(100,message)
  286. write(message,*) 'dy: ', config_flags%dy
  287. CALL wrf_debug(100,message)
  288. write(message,*) 'config_flags%start_year: ', config_flags%start_year
  289. CALL wrf_debug(100,message)
  290. write(message,*) 'config_flags%start_month: ', config_flags%start_month
  291. CALL wrf_debug(100,message)
  292. write(message,*) 'config_flags%start_day: ', config_flags%start_day
  293. CALL wrf_debug(100,message)
  294. write(message,*) 'config_flags%start_hour: ', config_flags%start_hour
  295. CALL wrf_debug(100,message)
  296. write(start_date,435) config_flags%start_year, config_flags%start_month, &
  297. config_flags%start_day, config_flags%start_hour
  298. 435 format(I4,'-',I2.2,'-',I2.2,'_',I2.2,':00:00')
  299. grid%dlmd=config_flags%dx
  300. grid%dphd=config_flags%dy
  301. tph0d=config_flags%cen_lat
  302. tlm0d=config_flags%cen_lon
  303. !==========================================================================
  304. !!
  305. ! Check to see if the boundary conditions are set
  306. ! properly in the namelist file.
  307. ! This checks for sufficiency and redundancy.
  308. CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
  309. ! Some sort of "this is the first time" initialization. Who knows.
  310. grid%itimestep=0
  311. ! Pull in the info in the namelist to compare it to the input data.
  312. grid%real_data_init_type = model_config_rec%real_data_init_type
  313. write(message,*) 'what is flag_metgrid: ', flag_metgrid
  314. CALL wrf_message(message)
  315. IF ( flag_metgrid .EQ. 1 ) THEN ! <----- START OF VERTICAL INTERPOLATION PART ---->
  316. num_metgrid_levels = grid%num_metgrid_levels
  317. IF (grid%ght_gc(its,jts,num_metgrid_levels/2) .lt. grid%ght_gc(its,jts,num_metgrid_levels/2+1)) THEN
  318. write(message,*) 'normal ground up file order'
  319. hyb_coor=.false.
  320. CALL wrf_message(message)
  321. ELSE
  322. hyb_coor=.true.
  323. write(message,*) 'reverse the order of coordinate'
  324. CALL wrf_message(message)
  325. CALL reverse_vert_coord(grid%ght_gc, 2, num_metgrid_levels &
  326. &, IDS,IDE,JDS,JDE,KDS,KDE &
  327. &, IMS,IME,JMS,JME,KMS,KME &
  328. &, ITS,ITE,JTS,JTE,KTS,KTE )
  329. #if defined(HWRF)
  330. if(.not. grid%use_prep_hybrid) then
  331. #endif
  332. CALL reverse_vert_coord(grid%p_gc, 2, num_metgrid_levels &
  333. &, IDS,IDE,JDS,JDE,KDS,KDE &
  334. &, IMS,IME,JMS,JME,KMS,KME &
  335. &, ITS,ITE,JTS,JTE,KTS,KTE )
  336. CALL reverse_vert_coord(grid%t_gc, 2, num_metgrid_levels &
  337. &, IDS,IDE,JDS,JDE,KDS,KDE &
  338. &, IMS,IME,JMS,JME,KMS,KME &
  339. &, ITS,ITE,JTS,JTE,KTS,KTE )
  340. CALL reverse_vert_coord(grid%u_gc, 2, num_metgrid_levels &
  341. &, IDS,IDE,JDS,JDE,KDS,KDE &
  342. &, IMS,IME,JMS,JME,KMS,KME &
  343. &, ITS,ITE,JTS,JTE,KTS,KTE )
  344. CALL reverse_vert_coord(grid%v_gc, 2, num_metgrid_levels &
  345. &, IDS,IDE,JDS,JDE,KDS,KDE &
  346. &, IMS,IME,JMS,JME,KMS,KME &
  347. &, ITS,ITE,JTS,JTE,KTS,KTE )
  348. CALL reverse_vert_coord(grid%rh_gc, 2, 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. #if defined(HWRF)
  353. endif
  354. #endif
  355. endif
  356. IF (hyb_coor) THEN
  357. ! limit extreme deviations from source model topography
  358. ! due to potential for nasty extrapolation/interpolation issues
  359. !
  360. write(message,*) 'min, max of grid%ht_gc before adjust: ', minval(grid%ht_gc), maxval(grid%ht_gc)
  361. CALL wrf_debug(100,message)
  362. ICOUNT=0
  363. DO J=JTS,min(JTE,JDE-1)
  364. DO I=ITS,min(ITE,IDE-1)
  365. IF ((grid%ht_gc(I,J) - grid%ght_gc(I,J,2)) .LT. -150.) THEN
  366. grid%ht_gc(I,J)=grid%ght_gc(I,J,2)-150.
  367. IF (ICOUNT .LT. 20) THEN
  368. write(message,*) 'increasing NMM topo toward RUC ', I,J
  369. CALL wrf_debug(100,message)
  370. ICOUNT=ICOUNT+1
  371. ENDIF
  372. ELSEIF ((grid%ht_gc(I,J) - grid%ght_gc(I,J,2)) .GT. 150.) THEN
  373. grid%ht_gc(I,J)=grid%ght_gc(I,J,2)+150.
  374. IF (ICOUNT .LT. 20) THEN
  375. write(message,*) 'decreasing NMM topo toward RUC ', I,J
  376. CALL wrf_debug(100,message)
  377. ICOUNT=ICOUNT+1
  378. ENDIF
  379. ENDIF
  380. END DO
  381. END DO
  382. write(message,*) 'min, max of ht_gc after correct: ', minval(grid%ht_gc), maxval(grid%ht_gc)
  383. CALL wrf_debug(100,message)
  384. ENDIF
  385. CALL boundary_smooth(grid%ht_gc,grid%landmask, grid, 12 , 12 &
  386. &, IDS,IDE,JDS,JDE,KDS,KDE &
  387. &, IMS,IME,JMS,JME,KMS,KME &
  388. &, ITS,ITE,JTS,JTE,KTS,KTE )
  389. DO j = jts, MIN(jte,jde-1)
  390. DO i = its, MIN(ite,ide-1)
  391. if (grid%landmask(I,J) .gt. 0.5) grid%sm(I,J)=0.
  392. if (grid%landmask(I,J) .le. 0.5) grid%sm(I,J)=1.
  393. if (grid%tsk_gc(I,J) .gt. 0.) then
  394. grid%nmm_tsk(I,J)=grid%tsk_gc(I,J)
  395. else
  396. #if defined(HWRF)
  397. if(grid%use_prep_hybrid) then
  398. if(grid%t(I,J,1)<100) then
  399. write(*,*) 'NO VALID SURFACE TEMPERATURE: I,J,TSK_GC(I,J),T(I,J,1) = ', &
  400. I,J,grid%TSK_GC(I,J),grid%T(I,J,1)
  401. else
  402. grid%nmm_tsk(I,J)=grid%t(I,J,1) ! stopgap measure
  403. end if
  404. else
  405. #endif
  406. grid%nmm_tsk(I,J)=grid%t_gc(I,J,1) ! stopgap measure
  407. #if defined(HWRF)
  408. endif
  409. #endif
  410. endif
  411. !
  412. grid%glat(I,J)=grid%hlat_gc(I,J)*DEGRAD
  413. grid%glon(I,J)=grid%hlon_gc(I,J)*DEGRAD
  414. grid%weasd(I,J)=grid%snow(I,J)
  415. grid%xice(I,J)=grid%xice_gc(I,J)
  416. ENDDO
  417. ENDDO
  418. ! First item is to define the target vertical coordinate
  419. num_metgrid_levels = grid%num_metgrid_levels
  420. eta_levels(1:kde) = model_config_rec%eta_levels(1:kde)
  421. ptsgm = model_config_rec%ptsgm
  422. p_top_requested = grid%p_top_requested
  423. grid%pt=p_top_requested
  424. if (internal_time_loop .eq. 1) then
  425. if (eta_levels(1) .ne. 1.0) then
  426. #if defined(HWRF)
  427. if(grid%use_prep_hybrid) then
  428. call wrf_error_fatal('PREP_HYBRID ERROR: eta_levels is not specified, but use_prep_hybrid=.true.')
  429. end if
  430. #endif
  431. write(message,*) '********************************************************************* '
  432. CALL wrf_message(message)
  433. write(message,*) '** eta_levels appears not to be specified in the namelist'
  434. CALL wrf_message(message)
  435. write(message,*) '** We will call compute_nmm_levels to define layer thicknesses.'
  436. CALL wrf_message(message)
  437. write(message,*) '** These levels should be reasonable for running the model, '
  438. CALL wrf_message(message)
  439. write(message,*) '** but may not be ideal for the simulation being made. Consider '
  440. CALL wrf_message(message)
  441. write(message,*) '** defining your own levels by specifying eta_levels in the model '
  442. CALL wrf_message(message)
  443. write(message,*) '** namelist. '
  444. CALL wrf_message(message)
  445. write(message,*) '********************************************************************** '
  446. CALL wrf_message(message)
  447. CALL compute_nmm_levels(KDE,p_top_requested,eta_levels)
  448. DO L=1,KDE
  449. write(message,*) 'L, eta_levels(L) returned :: ', L,eta_levels(L)
  450. CALL wrf_message(message)
  451. ENDDO
  452. endif
  453. write(message,*) 'KDE-1: ', KDE-1
  454. CALL wrf_debug(1,message)
  455. allocate(SG1(1:KDE-1))
  456. allocate(SG2(1:KDE-1))
  457. allocate(DSG1(1:KDE-1))
  458. allocate(DSG2(1:KDE-1))
  459. allocate(SGML1(1:KDE))
  460. allocate(SGML2(1:KDE))
  461. CALL define_nmm_vertical_coord (kde-1, ptsgm, grid%pt,grid%pdtop, eta_levels, &
  462. grid%eta1,grid%deta1,grid%aeta1, &
  463. grid%eta2,grid%deta2,grid%aeta2, grid%dfl, grid%dfrlg )
  464. DO L=KDS,KDE-1
  465. grid%deta(L)=eta_levels(L)-eta_levels(L+1)
  466. ENDDO
  467. endif
  468. write(message,*) 'num_metgrid_levels: ', num_metgrid_levels
  469. CALL wrf_message(message)
  470. DO j = jts, MIN(jte,jde-1)
  471. DO i = its, MIN(ite,ide-1)
  472. grid%fis(I,J)=grid%ht_gc(I,J)*g
  473. !
  474. ! 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
  475. 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
  476. IF (mod(I,10) .eq. 0 .and. mod(J,10) .eq. 0) THEN
  477. write(message,*) 'grid%ht_gc and grid%toposoil to swap, flag_soilhgt ::: ', &
  478. I,J, grid%ht_gc(I,J),grid%toposoil(I,J),flag_soilhgt
  479. CALL wrf_debug(10,message)
  480. ENDIF
  481. IF ( ( flag_soilhgt.EQ. 1 ) ) THEN
  482. grid%ght_gc(I,J,1)=grid%toposoil(I,J)
  483. ENDIF
  484. ENDIF
  485. ENDDO
  486. ENDDO
  487. numzero=0
  488. numexamined=0
  489. DO j = jts, MIN(jte,jde-1)
  490. DO i = its, MIN(ite,ide-1)
  491. numexamined=numexamined+1
  492. if(grid%fis(i,j)<1e-5 .and. grid%fis(i,j)>-1e5 ) then
  493. numzero=numzero+1
  494. end if
  495. enddo
  496. enddo
  497. write(message,*) 'TOTAL NEAR-ZERO FIS POINTS: ',numzero,' OF ',numexamined
  498. call wrf_debug(10,message)
  499. #if defined(HWRF)
  500. interp_notph: if(.not. grid%use_prep_hybrid) then
  501. #endif
  502. if (.NOT. allocated(PDVP)) allocate(PDVP(IMS:IME,JMS:JME))
  503. if (.NOT. allocated(P3D_OUT)) allocate(P3D_OUT(IMS:IME,JMS:JME,KDS:KDE-1))
  504. if (.NOT. allocated(PSFC_OUTV)) allocate(PSFC_OUTV(IMS:IME,JMS:JME))
  505. if (.NOT. allocated(P3DV_OUT)) allocate(P3DV_OUT(IMS:IME,JMS:JME,KDS:KDE-1))
  506. if (.NOT. allocated(P3DV_IN)) allocate(P3DV_IN(IMS:IME,JMS:JME,num_metgrid_levels))
  507. CALL compute_nmm_surfacep (grid%ht_gc, grid%ght_gc, grid%p_gc , grid%t_gc &
  508. &, grid%psfc_out, num_metgrid_levels &
  509. &, IDS,IDE,JDS,JDE,KDS,KDE &
  510. &, IMS,IME,JMS,JME,KMS,KME &
  511. &, ITS,ITE,JTS,JTE,KTS,KTE ) ! H points
  512. if (internal_time_loop .eq. 1) then
  513. write(message,*) 'psfc points (final combined)'
  514. loopinc=max( (JTE-JTS)/20,1)
  515. iloopinc=max( (ITE-ITS)/10,1)
  516. CALL wrf_message(message)
  517. DO J=min(JTE,JDE-1),JTS,-loopinc
  518. write(message,633) (grid%psfc_out(I,J)/100.,I=its,min(ite,IDE-1),iloopinc)
  519. CALL wrf_message(message)
  520. ENDDO
  521. endif
  522. 633 format(35(f5.0,1x))
  523. CALL compute_3d_pressure (grid%psfc_out,grid%aeta1,grid%aeta2 &
  524. &, grid%pdtop,grid%pt,grid%pd,p3d_out &
  525. &, IDS,IDE,JDS,JDE,KDS,KDE &
  526. &, IMS,IME,JMS,JME,KMS,KME &
  527. &, ITS,ITE,JTS,JTE,KTS,KTE )
  528. #ifdef DM_PARALLEL
  529. ips=its ; ipe=ite ; jps=jts ; jpe=jte ; kps=kts ; kpe=kte
  530. # include "HALO_NMM_MG2.inc"
  531. #endif
  532. #ifdef DM_PARALLEL
  533. # include "HALO_NMM_MG3.inc"
  534. #endif
  535. do K=1,num_metgrid_levels
  536. do J=JTS,min(JTE,JDE-1)
  537. do I=ITS,min(ITE,IDE-1)
  538. IF (K .eq. KTS) THEN
  539. IF (J .eq. JDS .and. I .lt. IDE-1) THEN ! S boundary
  540. PDVP(I,J)=0.5*(grid%pd(I,J)+grid%pd(I+1,J))
  541. PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J)+grid%psfc_out(I+1,J))
  542. ELSEIF (J .eq. JDE-1 .and. I .lt. IDE-1) THEN ! N boundary
  543. PDVP(I,J)=0.5*(grid%pd(I,J)+grid%pd(I+1,J))
  544. PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J)+grid%psfc_out(I+1,J))
  545. ELSEIF (I .eq. IDS .and. mod(J,2) .eq. 0) THEN ! W boundary
  546. PDVP(I,J)=0.5*(grid%pd(I,J-1)+grid%pd(I,J+1))
  547. PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J-1)+grid%psfc_out(I,J+1))
  548. ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 0) THEN ! E boundary
  549. PDVP(I,J)=0.5*(grid%pd(I,J-1)+grid%pd(I,J+1))
  550. PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J-1)+grid%psfc_out(I,J+1))
  551. ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 1) THEN ! phantom E boundary
  552. PDVP(I,J)=grid%pd(I,J)
  553. PSFC_OUTV(I,J)=grid%psfc_out(I,J)
  554. ELSEIF (mod(J,2) .eq. 0) THEN ! interior even row
  555. PDVP(I,J)=0.25*(grid%pd(I,J)+grid%pd(I-1,J)+grid%pd(I,J+1)+grid%pd(I,J-1))
  556. PSFC_OUTV(I,J)=0.25*(grid%psfc_out(I,J)+grid%psfc_out(I-1,J)+ &
  557. grid%psfc_out(I,J+1)+grid%psfc_out(I,J-1))
  558. ELSE ! interior odd row
  559. PDVP(I,J)=0.25*(grid%pd(I,J)+grid%pd(I+1,J)+grid%pd(I,J+1)+grid%pd(I,J-1))
  560. PSFC_OUTV(I,J)=0.25*(grid%psfc_out(I,J)+grid%psfc_out(I+1,J)+ &
  561. grid%psfc_out(I,J+1)+grid%psfc_out(I,J-1))
  562. ENDIF
  563. ENDIF
  564. IF (J .eq. JDS .and. I .lt. IDE-1) THEN ! S boundary
  565. P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J,K)+grid%p_gc(I+1,J,K))
  566. ELSEIF (J .eq. JDE-1 .and. I .lt. IDE-1) THEN ! N boundary
  567. P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J,K)+grid%p_gc(I+1,J,K))
  568. ELSEIF (I .eq. IDS .and. mod(J,2) .eq. 0) THEN ! W boundary
  569. P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J-1,K)+grid%p_gc(I,J+1,K))
  570. ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 0) THEN ! E boundary
  571. P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J-1,K)+grid%p_gc(I,J+1,K))
  572. ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 1) THEN ! phantom E boundary
  573. P3DV_IN(I,J,K)=grid%p_gc(I,J,K)
  574. ELSEIF (mod(J,2) .eq. 0) THEN ! interior even row
  575. P3DV_IN(I,J,K)=0.25*(grid%p_gc(I,J,K)+grid%p_gc(I-1,J,K) + &
  576. grid%p_gc(I,J+1,K)+grid%p_gc(I,J-1,K))
  577. ELSE ! interior odd row
  578. P3DV_IN(I,J,K)=0.25*(grid%p_gc(I,J,K)+grid%p_gc(I+1,J,K) + &
  579. grid%p_gc(I,J+1,K)+grid%p_gc(I,J-1,K))
  580. ENDIF
  581. enddo
  582. enddo
  583. enddo
  584. CALL compute_3d_pressure (psfc_outv,grid%aeta1,grid%aeta2 &
  585. &, grid%pdtop,grid%pt,pdvp,p3dv_out &
  586. &, IDS,IDE,JDS,JDE,KDS,KDE &
  587. &, IMS,IME,JMS,JME,KMS,KME &
  588. &, ITS,ITE,JTS,JTE,KTS,KTE )
  589. CALL interp_press2press_lin(grid%p_gc, p3d_out &
  590. &, grid%t_gc, grid%t,num_metgrid_levels &
  591. &, .TRUE.,.TRUE.,.TRUE. & ! extrap, ignore_lowest, t_field
  592. &, IDS,IDE,JDS,JDE,KDS,KDE &
  593. &, IMS,IME,JMS,JME,KMS,KME &
  594. &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
  595. CALL interp_press2press_lin(p3dv_in, p3dv_out &
  596. &, grid%u_gc, grid%u,num_metgrid_levels &
  597. &, .FALSE.,.TRUE.,.FALSE. &
  598. &, IDS,IDE,JDS,JDE,KDS,KDE &
  599. &, IMS,IME,JMS,JME,KMS,KME &
  600. &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
  601. CALL interp_press2press_lin(p3dv_in, p3dv_out &
  602. &, grid%v_gc, grid%v,num_metgrid_levels &
  603. &, .FALSE.,.TRUE.,.FALSE. &
  604. &, IDS,IDE,JDS,JDE,KDS,KDE &
  605. &, IMS,IME,JMS,JME,KMS,KME &
  606. &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
  607. IF (hyb_coor) THEN
  608. CALL wind_adjust(p3dv_in,p3dv_out,grid%u_gc,grid%v_gc,grid%u,grid%v &
  609. &, num_metgrid_levels,5000. &
  610. &, IDS,IDE,JDS,JDE,KDS,KDE &
  611. &, IMS,IME,JMS,JME,KMS,KME &
  612. &, ITS,ITE,JTS,JTE,KTS,KTE )
  613. ENDIF
  614. ALLOCATE(qtmp(IMS:IME,JMS:JME,num_metgrid_levels))
  615. ALLOCATE(qtmp2(IMS:IME,JMS:JME,num_metgrid_levels))
  616. CALL rh_to_mxrat (grid%rh_gc, grid%t_gc, grid%p_gc, qtmp , .TRUE. , &
  617. ids , ide , jds , jde , 1 , num_metgrid_levels , &
  618. ims , ime , jms , jme , 1 , num_metgrid_levels , &
  619. its , ite , jts , jte , 1 , num_metgrid_levels )
  620. do K=1,num_metgrid_levels
  621. do J=JTS,min(JTE,JDE-1)
  622. do I=ITS,min(ITE,IDE-1)
  623. QTMP2(I,J,K)=QTMP(I,J,K)/(1.0+QTMP(I,J,K))
  624. end do
  625. end do
  626. end do
  627. CALL interp_press2press_log(grid%p_gc, p3d_out &
  628. &, QTMP2, grid%q,num_metgrid_levels &
  629. &, .FALSE.,.TRUE. &
  630. &, IDS,IDE,JDS,JDE,KDS,KDE &
  631. &, IMS,IME,JMS,JME,KMS,KME &
  632. &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
  633. IF (ALLOCATED(QTMP)) DEALLOCATE(QTMP)
  634. IF (ALLOCATED(QTMP)) DEALLOCATE(QTMP2)
  635. #if defined(HWRF)
  636. else ! we are using prep_hybrid
  637. ! Compute surface pressure:
  638. grid%psfc_out=grid%pdtop+grid%pd
  639. end if interp_notph
  640. #endif
  641. ! Get the monthly values interpolated to the current date
  642. ! for the traditional monthly
  643. ! fields of green-ness fraction and background grid%albedo.
  644. if (internal_time_loop .eq. 1 .or. config_flags%sst_update .eq. 1) then
  645. CALL monthly_interp_to_date ( grid%greenfrac_gc , current_date , grid%vegfra , &
  646. ids , ide , jds , jde , kds , kde , &
  647. ims , ime , jms , jme , kms , kme , &
  648. its , ite , jts , jte , kts , kte )
  649. CALL monthly_interp_to_date ( grid%albedo12m_gc , current_date , grid%albbck , &
  650. ids , ide , jds , jde , kds , kde , &
  651. ims , ime , jms , jme , kms , kme , &
  652. its , ite , jts , jte , kts , kte )
  653. ! Get the min/max of each i,j for the monthly green-ness fraction.
  654. CALL monthly_min_max ( grid%greenfrac_gc , grid%shdmin , grid%shdmax , &
  655. ids , ide , jds , jde , kds , kde , &
  656. ims , ime , jms , jme , kms , kme , &
  657. its , ite , jts , jte , kts , kte )
  658. ! The model expects the green-ness values in percent, not fraction.
  659. DO j = jts, MIN(jte,jde-1)
  660. DO i = its, MIN(ite,ide-1)
  661. !! grid%vegfra(i,j) = grid%vegfra(i,j) * 100.
  662. grid%shdmax(i,j) = grid%shdmax(i,j) * 100.
  663. grid%shdmin(i,j) = grid%shdmin(i,j) * 100.
  664. grid%vegfrc(I,J)=grid%vegfra(I,J)
  665. END DO
  666. END DO
  667. ! The model expects the albedo fields as
  668. ! a fraction, not a percent. Set the water values to 8%.
  669. DO j = jts, MIN(jte,jde-1)
  670. DO i = its, MIN(ite,ide-1)
  671. if (grid%albbck(i,j) .lt. 5.) then
  672. write(message,*) 'reset albedo to 8%... I,J,albbck:: ', I,J,grid%albbck(I,J)
  673. CALL wrf_debug(10,message)
  674. grid%albbck(I,J)=8.
  675. endif
  676. grid%albbck(i,j) = grid%albbck(i,j) / 100.
  677. grid%snoalb(i,j) = grid%snoalb(i,j) / 100.
  678. IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
  679. grid%albbck(i,j) = 0.08
  680. grid%snoalb(i,j) = 0.08
  681. END IF
  682. grid%albase(i,j)=grid%albbck(i,j)
  683. grid%mxsnal(i,j)=grid%snoalb(i,j)
  684. END DO
  685. END DO
  686. endif
  687. #if defined(HWRF)
  688. if(.not.grid%use_prep_hybrid) then
  689. #endif
  690. ! new deallocs
  691. DEALLOCATE(p3d_out,p3dv_out,p3dv_in)
  692. #if defined(HWRF)
  693. end if
  694. #endif
  695. END IF ! <----- END OF VERTICAL INTERPOLATION PART ---->
  696. !! compute SST at each time if updating SST
  697. if (config_flags%sst_update == 1) then
  698. DO j = jts, MIN(jde-1,jte)
  699. DO i = its, MIN(ide-1,ite)
  700. if (grid%SM(I,J) .lt. 0.5) then
  701. grid%SST(I,J)=0.
  702. endif
  703. if (grid%SM(I,J) .gt. 0.5) then
  704. grid%SST(I,J)=grid%NMM_TSK(I,J)
  705. grid%NMM_TSK(I,J)=0.
  706. endif
  707. IF ( (grid%NMM_TSK(I,J)+grid%SST(I,J)) .lt. 200. .or. &
  708. (grid%NMM_TSK(I,J)+grid%SST(I,J)) .gt. 350. ) THEN
  709. write(message,*) 'TSK, SST trouble at : ', I,J
  710. CALL wrf_message(message)
  711. write(message,*) 'SM, NMM_TSK,SST ', grid%SM(I,J),grid%NMM_TSK(I,J),grid%SST(I,J)
  712. CALL wrf_message(message)
  713. ENDIF
  714. ENDDO
  715. ENDDO
  716. endif ! sst_update test
  717. if (internal_time_loop .eq. 1) then
  718. !!! weasd has "snow water equivalent" in mm
  719. DO j = jts, MIN(jte,jde-1)
  720. DO i = its, MIN(ite,ide-1)
  721. IF(grid%sm(I,J).GT.0.9) THEN
  722. IF (grid%xice(I,J) .gt. 0) then
  723. grid%si(I,J)=1.0
  724. ENDIF
  725. ! SEA
  726. grid%epsr(I,J)=.97
  727. grid%embck(I,J)=.97
  728. grid%gffc(I,J)=0.
  729. grid%albedo(I,J)=.06
  730. grid%albase(I,J)=.06
  731. IF(grid%si (I,J).GT.0. ) THEN
  732. ! SEA-ICE
  733. grid%sm(I,J)=0.
  734. grid%si(I,J)=0.
  735. grid%sice(I,J)=1.
  736. grid%gffc(I,J)=0. ! just leave zero as irrelevant
  737. grid%albedo(I,J)=.60
  738. grid%albase(I,J)=.60
  739. ENDIF
  740. ELSE
  741. grid%si(I,J)=5.0*grid%weasd(I,J)/1000.
  742. ! LAND
  743. grid%epsr(I,J)=1.0
  744. grid%embck(I,J)=1.0
  745. grid%gffc(I,J)=0.0 ! just leave zero as irrelevant
  746. grid%sice(I,J)=0.
  747. grid%sno(I,J)=grid%si(I,J)*.20
  748. ENDIF
  749. ENDDO
  750. ENDDO
  751. ! DETERMINE grid%albedo OVER LAND
  752. DO j = jts, MIN(jte,jde-1)
  753. DO i = its, MIN(ite,ide-1)
  754. IF(grid%sm(I,J).LT.0.9.AND.grid%sice(I,J).LT.0.9) THEN
  755. ! SNOWFREE albedo
  756. IF ( (grid%sno(I,J) .EQ. 0.0) .OR. &
  757. (grid%albase(I,J) .GE. grid%mxsnal(I,J) ) ) THEN
  758. grid%albedo(I,J) = grid%albase(I,J)
  759. ELSE
  760. ! MODIFY albedo IF SNOWCOVER:
  761. ! BELOW SNOWDEPTH THRESHOLD...
  762. IF (grid%sno(I,J) .LT. SNUP) THEN
  763. RSNOW = grid%sno(I,J)/SNUP
  764. SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP))
  765. ! ABOVE SNOWDEPTH THRESHOLD...
  766. ELSE
  767. SNOFAC = 1.0
  768. ENDIF
  769. ! CALCULATE grid%albedo ACCOUNTING FOR SNOWDEPTH AND VGFRCK
  770. grid%albedo(I,J) = grid%albase(I,J) &
  771. + (1.0-grid%vegfra(I,J))*SNOFAC*(grid%mxsnal(I,J)-grid%albase(I,J))
  772. ENDIF
  773. END IF
  774. grid%si(I,J)=5.0*grid%weasd(I,J)
  775. grid%sno(I,J)=grid%weasd(I,J)
  776. !! convert vegfra
  777. grid%vegfra(I,J)=grid%vegfra(I,J)*100.
  778. !
  779. ENDDO
  780. ENDDO
  781. #ifdef DM_PARALLEL
  782. ALLOCATE(SM_G(IDS:IDE,JDS:JDE),SICE_G(IDS:IDE,JDS:JDE))
  783. CALL WRF_PATCH_TO_GLOBAL_REAL( grid%sice(IMS,JMS) &
  784. &, SICE_G,grid%DOMDESC &
  785. &, 'z','xy' &
  786. &, IDS,IDE-1,JDS,JDE-1,1,1 &
  787. &, IMS,IME,JMS,JME,1,1 &
  788. &, ITS,ITE,JTS,JTE,1,1 )
  789. CALL WRF_PATCH_TO_GLOBAL_REAL( grid%sm(IMS,JMS) &
  790. &, SM_G,grid%DOMDESC &
  791. &, 'z','xy' &
  792. &, IDS,IDE-1,JDS,JDE-1,1,1 &
  793. &, IMS,IME,JMS,JME,1,1 &
  794. &, ITS,ITE,JTS,JTE,1,1 )
  795. IF (WRF_DM_ON_MONITOR()) THEN
  796. 637 format(40(f3.0,1x))
  797. allocate(IHE_G(JDS:JDE-1),IHW_G(JDS:JDE-1))
  798. DO j = JDS, JDE-1
  799. IHE_G(J)=MOD(J+1,2)
  800. IHW_G(J)=IHE_G(J)-1
  801. ENDDO
  802. DO ITER=1,10
  803. DO j = jds+1, (jde-1)-1
  804. DO i = ids+1, (ide-1)-1
  805. ! any sea ice around point in question?
  806. IF (SM_G(I,J) .ge. 0.9) THEN
  807. SEAICESUM=SICE_G(I+IHE_G(J),J+1)+SICE_G(I+IHW_G(J),J+1)+ &
  808. SICE_G(I+IHE_G(J),J-1)+SICE_G(I+IHW_G(J),J-1)
  809. IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN
  810. 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. &
  811. (SICE_G(I+IHW_G(J),J+1).lt.0.1 .and. SM_G(I+IHW_G(J),J+1).lt.0.1) .OR. &
  812. (SICE_G(I+IHE_G(J),J-1).lt.0.1 .and. SM_G(I+IHE_G(J),J-1).lt.0.1) .OR. &
  813. (SICE_G(I+IHW_G(J),J-1).lt.0.1 .and. SM_G(I+IHW_G(J),J-1).lt.0.1)) THEN
  814. ! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE
  815. write(message,*) 'making seaice (1): ', I,J
  816. CALL wrf_debug(100,message)
  817. SICE_G(I,J)=1.0
  818. SM_G(I,J)=0.
  819. ENDIF
  820. ELSEIF (SEAICESUM .ge. 3) THEN
  821. ! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE
  822. write(message,*) 'making seaice (2): ', I,J
  823. CALL wrf_debug(100,message)
  824. SICE_G(I,J)=1.0
  825. SM_G(I,J)=0.
  826. ENDIF
  827. ENDIF
  828. ENDDO
  829. ENDDO
  830. ENDDO
  831. ENDIF
  832. CALL WRF_GLOBAL_TO_PATCH_REAL( SICE_G, grid%sice &
  833. &, grid%DOMDESC &
  834. &, 'z','xy' &
  835. &, IDS,IDE-1,JDS,JDE-1,1,1 &
  836. &, IMS,IME,JMS,JME,1,1 &
  837. &, ITS,ITE,JTS,JTE,1,1 )
  838. CALL WRF_GLOBAL_TO_PATCH_REAL( SM_G,grid%sm &
  839. &, grid%DOMDESC &
  840. &, 'z','xy' &
  841. &, IDS,IDE-1,JDS,JDE-1,1,1 &
  842. &, IMS,IME,JMS,JME,1,1 &
  843. &, ITS,ITE,JTS,JTE,1,1 )
  844. IF (WRF_DM_ON_MONITOR()) THEN
  845. #if defined(HWRF)
  846. ! SM_G is still needed for the high-res grid
  847. #else
  848. DEALLOCATE(SM_G)
  849. #endif
  850. deallocate(SICE_G)
  851. DEALLOCATE(IHE_G,IHW_G)
  852. ENDIF
  853. ! write(message,*) 'revised sea ice on patch'
  854. ! CALL wrf_debug(100,message)
  855. ! DO J=JTE,JTS,-(((JTE-JTS)/25)+1)
  856. ! write(message,637) (grid%sice(I,J),I=ITS,ITE,ITE/20)
  857. ! CALL wrf_debug(100,message)
  858. ! END DO
  859. #else
  860. ! serial sea ice reprocessing
  861. allocate(IHE(JDS:JDE-1),IHW(JDS:JDE-1))
  862. DO j = jts, MIN(jte,jde-1)
  863. IHE(J)=MOD(J+1,2)
  864. IHW(J)=IHE(J)-1
  865. ENDDO
  866. DO ITER=1,10
  867. DO j = jts+1, MIN(jte,jde-1)-1
  868. DO i = its+1, MIN(ite,ide-1)-1
  869. ! any sea ice around point in question?
  870. IF (grid%sm(I,J) .gt. 0.9) THEN
  871. SEAICESUM=grid%sice(I+IHE(J),J+1)+grid%sice(I+IHW(J),J+1)+ &
  872. grid%sice(I+IHE(J),J-1)+grid%sice(I+IHW(J),J-1)
  873. IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN
  874. IF ((grid%sice(I+IHE(J),J+1).lt.0.1 .and. grid%sm(I+IHE(J),J+1).lt.0.1) .OR. &
  875. (grid%sice(I+IHW(J),J+1).lt.0.1 .and. grid%sm(I+IHW(J),J+1).lt.0.1) .OR. &
  876. (grid%sice(I+IHE(J),J-1).lt.0.1 .and. grid%sm(I+IHE(J),J-1).lt.0.1) .OR. &
  877. (grid%sice(I+IHW(J),J-1).lt.0.1 .and. grid%sm(I+IHW(J),J-1).lt.0.1)) THEN
  878. ! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE
  879. grid%sice(I,J)=1.0
  880. grid%sm(I,J)=0.
  881. ENDIF
  882. ELSEIF (SEAICESUM .ge. 3) THEN
  883. ! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE
  884. grid%sice(I,J)=1.0
  885. grid%sm(I,J)=0.
  886. ENDIF
  887. ENDIF
  888. ENDDO
  889. ENDDO
  890. ENDDO
  891. DEALLOCATE(IHE,IHW)
  892. #endif
  893. ! this block meant to guarantee land/sea agreement between sm and landmask
  894. DO j = jts, MIN(jte,jde-1)
  895. DO i = its, MIN(ite,ide-1)
  896. IF (grid%sm(I,J) .gt. 0.5) THEN
  897. grid%landmask(I,J)=0.0
  898. ELSEIF (grid%sm(I,J) .lt. 0.5 .and. grid%sice(I,J) .gt. 0.9) then
  899. grid%landmask(I,J)=0.0
  900. ELSEIF (grid%sm(I,J) .lt. 0.5 .and. grid%sice(I,J) .lt. 0.1) then
  901. grid%landmask(I,J)=1.0
  902. ELSE
  903. write(message,*) 'missed point in grid%landmask definition ' , I,J
  904. CALL wrf_message(message)
  905. grid%landmask(I,J)=0.0
  906. ENDIF
  907. !
  908. 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
  909. write(message,*) 'set grid%nmm_tsk to: ', grid%sst(I,J)
  910. CALL wrf_message(message)
  911. grid%nmm_tsk(I,J)=grid%sst(I,J)
  912. grid%sst(I,J)=0.
  913. endif
  914. ENDDO
  915. ENDDO
  916. ! For sf_surface_physics = 1, we want to use close to a 10 cm value
  917. ! for the bottom level of the soil temps.
  918. IF ( ( model_config_rec%sf_surface_physics(grid%id) .EQ. 1 ) .AND. &
  919. ( flag_st000010 .EQ. 1 ) ) THEN
  920. DO j = jts , MIN(jde-1,jte)
  921. DO i = its , MIN(ide-1,ite)
  922. grid%soiltb(i,j) = grid%st000010(i,j)
  923. END DO
  924. END DO
  925. END IF
  926. ! Adjust the various soil temperature values depending on the difference in
  927. ! in elevation between the current model's elevation and the incoming data's
  928. ! orography.
  929. IF ( ( flag_toposoil .EQ. 1 ) ) THEN
  930. ALLOCATE(HT(ims:ime,jms:jme))
  931. DO J=jms,jme
  932. DO I=ims,ime
  933. HT(I,J)=grid%fis(I,J)/9.81
  934. END DO
  935. END DO
  936. ! if (maxval(grid%toposoil) .gt. 100.) then
  937. !
  938. ! Being avoided. Something to revisit eventually.
  939. !
  940. !1219 might be simply a matter of including toposoil
  941. !
  942. ! CODE NOT TESTED AT NCEP USING THIS FUNCTIONALITY,
  943. ! SO TO BE SAFE WILL AVOID FOR RETRO RUNS.
  944. !
  945. ! CALL adjust_soil_temp_new ( grid%soiltb , 2 , &
  946. ! grid%nmm_tsk , ht , grid%toposoil , grid%landmask, flag_toposoil , &
  947. ! grid%st000010 , st010040 , st040100 , st100200 , st010200 , &
  948. ! flag_st000010 , flag_st010040 , flag_st040100 , &
  949. ! flag_st100200 , flag_st010200 , &
  950. ! soilt000 , soilt005 , soilt020 , soilt040 , &
  951. ! soilt160 , soilt300 , &
  952. ! flag_soilt000 , flag_soilt005 , flag_soilt020 , &
  953. ! flag_soilt040 , flag_soilt160 , flag_soilt300 , &
  954. ! ids , ide , jds , jde , kds , kde , &
  955. ! ims , ime , jms , jme , kms , kme , &
  956. ! its , ite , jts , jte , kts , kte )
  957. ! endif
  958. END IF
  959. ! Process the LSM data.
  960. ! surface_input_source=1 => use data from static file
  961. ! (fractional category as input)
  962. ! surface_input_source=2 => use data from grib file
  963. ! (dominant category as input)
  964. IF ( config_flags%surface_input_source .EQ. 1 ) THEN
  965. grid%vegcat (its,jts) = 0
  966. grid%soilcat(its,jts) = 0
  967. END IF
  968. ! Generate the vegetation and soil category information
  969. ! from the fractional input
  970. ! data, or use the existing dominant category fields if they exist.
  971. IF ((grid%soilcat(its,jts) .LT. 0.5) .AND. (grid%vegcat(its,jts) .LT. 0.5)) THEN
  972. num_veg_cat = SIZE ( grid%landusef_gc , DIM=3 )
  973. num_soil_top_cat = SIZE ( grid%soilctop_gc , DIM=3 )
  974. num_soil_bot_cat = SIZE ( grid%soilcbot_gc , DIM=3 )
  975. do J=JMS,JME
  976. do K=1,num_veg_cat
  977. do I=IMS,IME
  978. grid%landusef(I,K,J)=grid%landusef_gc(I,J,K)
  979. enddo
  980. enddo
  981. enddo
  982. do J=JMS,JME
  983. do K=1,num_soil_top_cat
  984. do I=IMS,IME
  985. grid%soilctop(I,K,J)=grid%soilctop_gc(I,J,K)
  986. enddo
  987. enddo
  988. enddo
  989. do J=JMS,JME
  990. do K=1,num_soil_bot_cat
  991. do I=IMS,IME
  992. grid%soilcbot(I,K,J)=grid%soilcbot_gc(I,J,K)
  993. enddo
  994. enddo
  995. enddo
  996. ! grid%sm (1=water, 0=land)
  997. ! grid%landmask(0=water, 1=land)
  998. write(message,*) 'landmask into process_percent_cat_new'
  999. CALL wrf_debug(1,message)
  1000. do J=JTE,JTS,-(((JTE-JTS)/20)+1)
  1001. write(message,641) (grid%landmask(I,J),I=ITS,min(ITE,IDE-1),((ITE-ITS)/15)+1)
  1002. CALL wrf_debug(1,message)
  1003. enddo
  1004. 641 format(25(f3.0,1x))
  1005. CALL process_percent_cat_new ( grid%landmask , &
  1006. grid%landusef , grid%soilctop , grid%soilcbot , &
  1007. grid%isltyp , grid%ivgtyp , &
  1008. num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
  1009. ids , ide , jds , jde , kds , kde , &
  1010. ims , ime , jms , jme , kms , kme , &
  1011. its , ite , jts , jte , kts , kte , &
  1012. model_config_rec%iswater(grid%id) )
  1013. DO j = jts , MIN(jde-1,jte)
  1014. DO i = its , MIN(ide-1,ite)
  1015. grid%vegcat(i,j) = grid%ivgtyp(i,j)
  1016. grid%soilcat(i,j) = grid%isltyp(i,j)
  1017. END DO
  1018. END DO
  1019. ELSE
  1020. ! Do we have dominant soil and veg data from the input already?
  1021. IF ( grid%soilcat(its,jts) .GT. 0.5 ) THEN
  1022. DO j = jts, MIN(jde-1,jte)
  1023. DO i = its, MIN(ide-1,ite)
  1024. grid%isltyp(i,j) = NINT( grid%soilcat(i,j) )
  1025. END DO
  1026. END DO
  1027. END IF
  1028. IF ( grid%vegcat(its,jts) .GT. 0.5 ) THEN
  1029. DO j = jts, MIN(jde-1,jte)
  1030. DO i = its, MIN(ide-1,ite)
  1031. grid%ivgtyp(i,j) = NINT( grid%vegcat(i,j) )
  1032. END DO
  1033. END DO
  1034. END IF
  1035. ENDIF
  1036. DO j = jts, MIN(jde-1,jte)
  1037. DO i = its, MIN(ide-1,ite)
  1038. IF (grid%sice(I,J) .lt. 0.1) THEN
  1039. IF (grid%landmask(I,J) .gt. 0.5 .and. grid%sm(I,J) .gt. 0.5) THEN
  1040. write(message,*) 'land mask and grid%sm both > 0.5: ', &
  1041. I,J,grid%landmask(I,J),grid%sm(I,J)
  1042. CALL wrf_message(message)
  1043. grid%sm(I,J)=0.
  1044. ELSEIF (grid%landmask(I,J) .lt. 0.5 .and. grid%sm(I,J) .lt. 0.5) THEN
  1045. write(message,*) 'land mask and grid%sm both < 0.5: ', &
  1046. I,J, grid%landmask(I,J),grid%sm(I,J)
  1047. CALL wrf_message(message)
  1048. grid%sm(I,J)=1.
  1049. ENDIF
  1050. ELSE
  1051. IF (grid%landmask(I,J) .gt. 0.5 .and. grid%sm(I,J)+grid%sice(I,J) .gt. 0.9) then
  1052. write(message,*) 'landmask says LAND, sm/sice say SEAICE: ', I,J
  1053. ENDIF
  1054. ENDIF
  1055. ENDDO
  1056. ENDDO
  1057. DO j = jts, MIN(jde-1,jte)
  1058. DO i = its, MIN(ide-1,ite)
  1059. if (grid%sice(I,J) .gt. 0.9) then
  1060. grid%isltyp(I,J)=16
  1061. grid%ivgtyp(I,J)=24
  1062. endif
  1063. ENDDO
  1064. ENDDO
  1065. DO j = jts, MIN(jde-1,jte)
  1066. DO i = its, MIN(ide-1,ite)
  1067. if (grid%sm(I,J) .lt. 0.5) then
  1068. grid%sst(I,J)=0.
  1069. endif
  1070. if (grid%sm(I,J) .gt. 0.5) then
  1071. if (grid%sst(I,J) .lt. 0.1) then
  1072. grid%sst(I,J)=grid%nmm_tsk(I,J)
  1073. endif
  1074. grid%nmm_tsk(I,J)=0.
  1075. endif
  1076. IF ( (grid%nmm_tsk(I,J)+grid%sst(I,J)) .lt. 200. .or. &
  1077. (grid%nmm_tsk(I,J)+grid%sst(I,J)) .gt. 350. ) THEN
  1078. write(message,*) 'TSK, sst trouble at : ', I,J
  1079. CALL wrf_message(message)
  1080. write(message,*) 'sm, nmm_tsk,sst ', grid%sm(I,J),grid%nmm_tsk(I,J),grid%sst(I,J)
  1081. CALL wrf_message(message)
  1082. ENDIF
  1083. ENDDO
  1084. ENDDO
  1085. write(message,*) 'grid%sm'
  1086. CALL wrf_message(message)
  1087. DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
  1088. write(message,635) (grid%sm(i,J),I=its,ite,((ite-its)/10)+1)
  1089. CALL wrf_message(message)
  1090. END DO
  1091. write(message,*) 'sst/nmm_tsk'
  1092. CALL wrf_debug(10,message)
  1093. DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
  1094. write(message,635) (grid%sst(I,J)+grid%nmm_tsk(I,J),I=ITS,min(ide-1,ite),((ite-its)/10)+1)
  1095. CALL wrf_debug(10,message)
  1096. END DO
  1097. 635 format(20(f5.1,1x))
  1098. DO j = jts, MIN(jde-1,jte)
  1099. DO i = its, MIN(ide-1,ite)
  1100. IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
  1101. grid%soiltb(i,j) = grid%sst(i,j)
  1102. ELSE IF ( grid%landmask(i,j) .GT. 0.5 ) THEN
  1103. grid%soiltb(i,j) = grid%nmm_tsk(i,j)
  1104. END IF
  1105. END DO
  1106. END DO
  1107. ! END IF
  1108. ! Land use categories, dominant soil and vegetation types (if available).
  1109. ! allocate(grid%lu_index(ims:ime,jms:jme))
  1110. DO j = jts, MIN(jde-1,jte)
  1111. DO i = its, MIN(ide-1,ite)
  1112. grid%lu_index(i,j) = grid%ivgtyp(i,j)
  1113. END DO
  1114. END DO
  1115. if (flag_sst .eq. 1) log_flag_sst=.true.
  1116. if (flag_sst .eq. 0) log_flag_sst=.false.
  1117. write(message,*) 'st_input dimensions: ', size(st_input,dim=1), &
  1118. size(st_input,dim=2),size(st_input,dim=3)
  1119. CALL wrf_debug(100,message)
  1120. ! write(message,*) 'maxval st_input(1): ', maxval(st_input(:,1,:))
  1121. ! CALL wrf_message(message)
  1122. ! write(message,*) 'maxval st_input(2): ', maxval(st_input(:,2,:))
  1123. ! CALL wrf_message(message)
  1124. ! write(message,*) 'maxval st_input(3): ', maxval(st_input(:,3,:))
  1125. ! CALL wrf_message(message)
  1126. ! write(message,*) 'maxval st_input(4): ', maxval(st_input(:,4,:))
  1127. ! CALL wrf_message(message)
  1128. ! =============================================================
  1129. IF (.NOT. ALLOCATED(TG_ALT))ALLOCATE(TG_ALT(grid%sm31:grid%em31,jms:jme))
  1130. TPH0=TPH0D*DTR
  1131. WBD=-(((ide-1)-1)*grid%dlmd)
  1132. WB= WBD*DTR
  1133. SBD=-(((jde-1)/2)*grid%dphd)
  1134. SB= SBD*DTR
  1135. DLM=grid%dlmd*DTR
  1136. DPH=grid%dphd*DTR
  1137. TDLM=DLM+DLM
  1138. TDPH=DPH+DPH
  1139. WBI=WB+TDLM
  1140. SBI=SB+TDPH
  1141. EBI=WB+(ide-2)*TDLM
  1142. ANBI=SB+(jde-2)*DPH
  1143. STPH0=SIN(TPH0)
  1144. CTPH0=COS(TPH0)
  1145. TSPH=3600./GRID%DT
  1146. DO J=JTS,min(JTE,JDE-1)
  1147. TLM=WB-TDLM+MOD(J,2)*DLM !For velocity points on the E grid
  1148. TPH=SB+float(J-1)*DPH
  1149. STPH=SIN(TPH)
  1150. CTPH=COS(TPH)
  1151. DO I=ITS,MIN(ITE,IDE-1)
  1152. if (I .eq. ITS) THEN
  1153. TLM=TLM+TDLM*ITS
  1154. else
  1155. TLM=TLM+TDLM
  1156. endif
  1157. TERM1=(STPH0*CTPH*COS(TLM)+CTPH0*STPH)
  1158. FP=TWOM*(TERM1)
  1159. grid%f(I,J)=0.5*GRID%DT*FP
  1160. ENDDO
  1161. ENDDO
  1162. DO J=JTS,min(JTE,JDE-1)
  1163. TLM=WB-TDLM+MOD(J+1,2)*DLM !For mass points on the E grid
  1164. TPH=SB+float(J-1)*DPH
  1165. STPH=SIN(TPH)
  1166. CTPH=COS(TPH)
  1167. DO I=ITS,MIN(ITE,IDE-1)
  1168. if (I .eq. ITS) THEN
  1169. TLM=TLM+TDLM*ITS
  1170. else
  1171. TLM=TLM+TDLM
  1172. endif
  1173. TERM1=(STPH0*CTPH*COS(TLM)+CTPH0*STPH)
  1174. TERM1=MIN(TERM1,1.0D0)
  1175. TERM1=MAX(TERM1,-1.0D0)
  1176. APH=ASIN(TERM1)
  1177. TG_ALT(I,J)=TG0+TGA*COS(APH)-grid%fis(I,J)/3333.
  1178. ENDDO
  1179. ENDDO
  1180. DO j = jts, MIN(jde-1,jte)
  1181. DO i = its, MIN(ide-1,ite)
  1182. ! IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
  1183. ! grid%sice(I,J) .eq. 0. ) THEN
  1184. ! grid%tg(i,j) = grid%sst(i,j)
  1185. ! ELSEIF (grid%sice(I,J) .eq. 1) THEN
  1186. ! grid%tg(i,j) = 271.16
  1187. ! END IF
  1188. if (grid%tg(I,J) .lt. 200.) then ! only use default TG_ALT definition if
  1189. ! not getting TGROUND from grid%si
  1190. grid%tg(I,J)=TG_ALT(I,J)
  1191. endif
  1192. if (grid%tg(I,J) .lt. 200. .or. grid%tg(I,J) .gt. 320.) then
  1193. write(message,*) 'problematic grid%tg point at : ', I,J
  1194. CALL wrf_message( message )
  1195. endif
  1196. adum2d(i,j)=grid%nmm_tsk(I,J)+grid%sst(I,J)
  1197. END DO
  1198. END DO
  1199. DEALLOCATE(TG_ALT)
  1200. write(message,*) 'call process_soil_real with num_st_levels_input: ', num_st_levels_input
  1201. CALL wrf_message( message )
  1202. ! =============================================================
  1203. CALL process_soil_real ( adum2d, grid%tg , &
  1204. grid%landmask, grid%sst, &
  1205. st_input, sm_input, sw_input, &
  1206. st_levels_input , sm_levels_input , &
  1207. sw_levels_input , &
  1208. grid%sldpth , grid%dzsoil , grid%stc , grid%smc , grid%sh2o, &
  1209. flag_sst , flag_soilt000, flag_soilm000, &
  1210. ids , ide , jds , jde , kds , kde , &
  1211. ims , ime , jms , jme , kms , kme , &
  1212. its , ite , jts , jte , kts , kte , &
  1213. model_config_rec%sf_surface_physics(grid%id) , &
  1214. model_config_rec%num_soil_layers , &
  1215. model_config_rec%real_data_init_type , &
  1216. num_st_levels_input , num_sm_levels_input , &
  1217. num_sw_levels_input , &
  1218. num_st_levels_alloc , num_sm_levels_alloc , &
  1219. num_sw_levels_alloc )
  1220. ! =============================================================
  1221. ! Minimum soil values, residual, from RUC LSM scheme.
  1222. ! For input from Noah and using
  1223. ! RUC LSM scheme, this must be subtracted from the input
  1224. ! total soil moisture. For input RUC data and using the Noah LSM scheme,
  1225. ! this value must be added to the soil moisture_input.
  1226. lqmi(1:num_soil_top_cat) = &
  1227. (/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
  1228. 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
  1229. 0.004, 0.065 /) !dusan , 0.020, 0.004, 0.008 /)
  1230. ! At the initial time we care about values of soil moisture and temperature,
  1231. ! other times are ignored by the model, so we ignore them, too.
  1232. account_for_zero_soil_moisture : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
  1233. CASE ( LSMSCHEME )
  1234. iicount = 0
  1235. IF ( FLAG_SM000010 .EQ. 1 ) THEN
  1236. DO j = jts, MIN(jde-1,jte)
  1237. DO i = its, MIN(ide-1,ite)
  1238. IF ((grid%landmask(i,j).gt.0.5) .and. (grid%stc(i,1,j) .gt. 200) .and. &
  1239. (grid%stc(i,1,j) .lt. 400) .and. (grid%smc(i,1,j) .lt. 0.005)) then
  1240. write(message,*) 'Noah > Noah: bad soil moisture at i,j = ',i,j,grid%smc(i,:,j)
  1241. CALL wrf_message(message)
  1242. iicount = iicount + 1
  1243. grid%smc(i,:,j) = 0.005
  1244. END IF
  1245. END DO
  1246. END DO
  1247. IF ( iicount .GT. 0 ) THEN
  1248. write(message,*) 'Noah -> Noah: total number of small soil moisture locations= ',&
  1249. iicount
  1250. CALL wrf_message(message)
  1251. END IF
  1252. ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
  1253. DO j = jts, MIN(jde-1,jte)
  1254. DO i = its, MIN(ide-1,ite)
  1255. grid%smc(i,:,j) = grid%smc(i,:,j) + lqmi(grid%isltyp(i,j))
  1256. END DO
  1257. END DO
  1258. DO j = jts, MIN(jde-1,jte)
  1259. DO i = its, MIN(ide-1,ite)
  1260. IF ((grid%landmask(i,j).gt.0.5) .and. (grid%stc(i,1,j) .gt. 200) .and. &
  1261. (grid%stc(i,1,j) .lt. 400) .and. (grid%smc(i,1,j) .lt. 0.004)) then
  1262. write(message,*) 'RUC -> Noah: bad soil moisture at i,j = ' &
  1263. ,i,j,grid%smc(i,:,j)
  1264. CALL wrf_message(message)
  1265. iicount = iicount + 1
  1266. grid%smc(i,:,j) = 0.004
  1267. END IF
  1268. END DO
  1269. END DO
  1270. IF ( iicount .GT. 0 ) THEN
  1271. write(message,*) 'RUC -> Noah: total number of small soil moisture locations = ',&
  1272. iicount
  1273. CALL wrf_message(message)
  1274. END IF
  1275. END IF
  1276. CASE ( RUCLSMSCHEME )
  1277. iicount = 0
  1278. IF ( FLAG_SM000010 .EQ. 1 ) THEN
  1279. DO j = jts, MIN(jde-1,jte)
  1280. DO i = its, MIN(ide-1,ite)
  1281. grid%smc(i,:,j) = MAX ( grid%smc(i,:,j) - lqmi(grid%isltyp(i,j)) , 0. )
  1282. END DO
  1283. END DO
  1284. ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
  1285. ! no op
  1286. END IF
  1287. END SELECT account_for_zero_soil_moisture
  1288. !!! zero out grid%nmm_tsk at water points again
  1289. DO j = jts, MIN(jde-1,jte)
  1290. DO i = its, MIN(ide-1,ite)
  1291. if (grid%sm(I,J) .gt. 0.5) then
  1292. grid%nmm_tsk(I,J)=0.
  1293. endif
  1294. END DO
  1295. END DO
  1296. !! check on grid%stc
  1297. DO j = jts, MIN(jde-1,jte)
  1298. DO i = its, MIN(ide-1,ite)
  1299. IF (grid%sice(I,J) .gt. 0.9) then
  1300. DO L = 1, grid%num_soil_layers
  1301. grid%stc(I,L,J)=271.16 ! grid%tg value used by Eta/NMM
  1302. END DO
  1303. END IF
  1304. IF (grid%sm(I,J) .gt. 0.9) then
  1305. DO L = 1, grid%num_soil_layers
  1306. grid%stc(I,L,J)=273.16 ! grid%tg value used by Eta/NMM
  1307. END DO
  1308. END IF
  1309. END DO
  1310. END DO
  1311. DO j = jts, MIN(jde-1,jte)
  1312. DO i = its, MIN(ide-1,ite)
  1313. if (grid%sm(I,J) .lt. 0.1 .and. grid%stc(I,1,J) .lt. 0.1) THEN
  1314. write(message,*) 'troublesome grid%sm,grid%stc,grid%smc value: ', I,J,grid%sm(I,J), grid%stc(I,1,J),grid%smc(I,1,J)
  1315. CALL wrf_message(message)
  1316. do JJ=J-1,J+1
  1317. do L=1, grid%num_soil_layers
  1318. do II=I-1,I+1
  1319. if (II .ge. its .and. II .le. MIN(ide-1,ite) .and. &
  1320. JJ .ge. jts .and. JJ .le. MIN(jde-1,jte)) then
  1321. grid%stc(I,L,J)=amax1(grid%stc(I,L,J),grid%stc(II,L,JJ))
  1322. cur_smc=grid%smc(I,L,J)
  1323. if ( grid%smc(II,L,JJ) .gt. 0.005 .and. grid%smc(II,L,JJ) .lt. 1.0) then
  1324. aposs_smc=grid%smc(II,L,JJ)
  1325. if ( cur_smc .eq. 0 ) then
  1326. cur_smc=aposs_smc
  1327. grid%smc(I,L,J)=cur_smc
  1328. else
  1329. cur_smc=amin1(cur_smc,aposs_smc)
  1330. cur_smc=amin1(cur_smc,aposs_smc)
  1331. grid%smc(I,L,J)=cur_smc
  1332. endif
  1333. endif
  1334. endif ! bounds check
  1335. enddo
  1336. enddo
  1337. enddo
  1338. write(message,*) 'grid%stc, grid%smc(1) now: ', grid%stc(I,1,J),grid%smc(I,1,J)
  1339. CALL wrf_message(message)
  1340. endif
  1341. if (grid%stc(I,1,J) .lt. 0.1) then
  1342. write(message,*) 'QUITTING DUE TO STILL troublesome grid%stc value: ', I,J, grid%stc(I,1,J),grid%smc(I,1,J)
  1343. call wrf_error_fatal(message)
  1344. endif
  1345. ENDDO
  1346. ENDDO
  1347. !hardwire soil stuff for time being
  1348. ! RTDPTH=0.
  1349. ! RTDPTH(1)=0.1
  1350. ! RTDPTH(2)=0.3
  1351. ! RTDPTH(3)=0.6
  1352. ! grid%sldpth=0.
  1353. ! grid%sldpth(1)=0.1
  1354. ! grid%sldpth(2)=0.3
  1355. ! grid%sldpth(3)=0.6
  1356. ! grid%sldpth(4)=1.0
  1357. !!! main body of nmm_specific starts here
  1358. !
  1359. do J=jts,min(jte,jde-1)
  1360. do I=its,min(ite,ide-1)
  1361. grid%res(I,J)=1.
  1362. enddo
  1363. enddo
  1364. !! grid%hbm2
  1365. grid%hbm2=0.
  1366. do J=jts,min(jte,jde-1)
  1367. do I=its,min(ite,ide-1)
  1368. IF ( (J .ge. 3 .and. J .le. (jde-1)-2) .AND. &
  1369. (I .ge. 2 .and. I .le. (ide-1)-2+mod(J,2)) ) THEN
  1370. grid%hbm2(I,J)=1.
  1371. ENDIF
  1372. enddo
  1373. enddo
  1374. !! grid%hbm3
  1375. grid%hbm3=0.
  1376. !! LOOP OVER LOCAL DIMENSIONS
  1377. do J=jts,min(jte,jde-1)
  1378. grid%ihwg(J)=mod(J+1,2)-1
  1379. IF (J .ge. 4 .and. J .le. (jde-1)-3) THEN
  1380. IHL=(ids+1)-grid%ihwg(J)
  1381. IHH=(ide-1)-2
  1382. do I=its,min(ite,ide-1)
  1383. IF (I .ge. IHL .and. I .le. IHH) grid%hbm3(I,J)=1.
  1384. enddo
  1385. ENDIF
  1386. enddo
  1387. !! grid%vbm2
  1388. grid%vbm2=0.
  1389. do J=jts,min(jte,jde-1)
  1390. do I=its,min(ite,ide-1)
  1391. IF ( (J .ge. 3 .and. J .le. (jde-1)-2) .AND. &
  1392. (I .ge. 2 .and. I .le. (ide-1)-1-mod(J,2)) ) THEN
  1393. grid%vbm2(I,J)=1.
  1394. ENDIF
  1395. enddo
  1396. enddo
  1397. !! grid%vbm3
  1398. grid%vbm3=0.
  1399. do J=jts,min(jte,jde-1)
  1400. do I=its,min(ite,ide-1)
  1401. IF ( (J .ge. 4 .and. J .le. (jde-1)-3) .AND. &
  1402. (I .ge. 3-mod(J,2) .and. I .le. (ide-1)-2) ) THEN
  1403. grid%vbm3(I,J)=1.
  1404. ENDIF
  1405. enddo
  1406. enddo
  1407. COAC=model_config_rec%coac(grid%id)
  1408. CODAMP=model_config_rec%codamp(grid%id)
  1409. DTAD=1.0
  1410. ! IDTCF=DTCF, IDTCF=4
  1411. DTCF=4.0 ! used?
  1412. grid%dy_nmm=ERAD*DPH
  1413. grid%cpgfv=-GRID%DT/(48.*grid%dy_nmm)
  1414. grid%en= GRID%DT/( 4.*grid%dy_nmm)*DTAD
  1415. grid%ent=GRID%DT/(16.*grid%dy_nmm)*DTAD
  1416. DO J=jts,nnyp
  1417. KHL2(J)=(IDE-1)*(J-1)-(J-1)/2+2
  1418. KVL2(J)=(IDE-1)*(J-1)-J/2+2
  1419. KHH2(J)=(IDE-1)*J-J/2-1
  1420. KVH2(J)=(IDE-1)*J-(J+1)/2-1
  1421. ENDDO
  1422. TPH=SB-DPH
  1423. DO J=jts,min(jte,jde-1)
  1424. TPH=SB+float(J-1)*DPH
  1425. DXP=ERAD*DLM*COS(TPH)
  1426. DXJ(J)=DXP
  1427. WPDARJ(J)=-W_NMM * &
  1428. ((ERAD*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+grid%dy_nmm**2)/ &
  1429. (GRID%DT*32.*DXP*grid%dy_nmm)
  1430. CPGFUJ(J)=-GRID%DT/(48.*DXP)
  1431. CURVJ(J)=.5*GRID%DT*TAN(TPH)/ERAD
  1432. FCPJ(J)=GRID%DT/(CP*192.*DXP*grid%dy_nmm)
  1433. FDIVJ(J)=1./(12.*DXP*grid%dy_nmm)
  1434. ! EMJ(J)= GRID%DT/( 4.*DXP)*DTAD
  1435. ! EMTJ(J)=GRID%DT/(16.*DXP)*DTAD
  1436. FADJ(J)=-GRID%DT/(48.*DXP*grid%dy_nmm)*DTAD
  1437. ACDT=GRID%DT*SQRT((ERAD*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+grid%dy_nmm**2)
  1438. CDDAMP=CODAMP*ACDT
  1439. HDACJ(J)=COAC*ACDT/(4.*DXP*grid%dy_nmm)
  1440. DDMPUJ(J)=CDDAMP/DXP
  1441. DDMPVJ(J)=CDDAMP/grid%dy_nmm
  1442. ENDDO
  1443. DO J=JTS,min(JTE,JDE-1)
  1444. TLM=WB-TDLM+MOD(J,2)*DLM
  1445. TPH=SB+float(J-1)*DPH
  1446. STPH=SIN(TPH)
  1447. CTPH=COS(TPH)
  1448. DO I=ITS,MIN(ITE,IDE-1)
  1449. if (I .eq. ITS) THEN
  1450. TLM=TLM+TDLM*ITS
  1451. else
  1452. TLM=TLM+TDLM
  1453. endif
  1454. FP=TWOM*(CTPH0*STPH+STPH0*CTPH*COS(TLM))
  1455. grid%f(I,J)=0.5*GRID%DT*FP
  1456. ENDDO
  1457. ENDDO
  1458. ! --------------DERIVED VERTICAL GRID CONSTANTS--------------------------
  1459. grid%ef4t=.5*GRID%DT/CP
  1460. grid%f4q = -GRID%DT*DTAD
  1461. grid%f4d =-.5*GRID%DT*DTAD
  1462. DO L=KDS,KDE-1
  1463. grid%rdeta(L)=1./grid%deta(L)
  1464. grid%f4q2(L)=-.25*GRID%DT*DTAD/grid%deta(L)
  1465. ENDDO
  1466. DO J=JTS,min(JTE,JDE-1)
  1467. DO I=ITS,min(ITE,IDE-1)
  1468. grid%dx_nmm(I,J)=DXJ(J)
  1469. grid%wpdar(I,J)=WPDARJ(J)*grid%hbm2(I,J)
  1470. grid%cpgfu(I,J)=CPGFUJ(J)*grid%vbm2(I,J)
  1471. grid%curv(I,J)=CURVJ(J)*grid%vbm2(I,J)
  1472. grid%fcp(I,J)=FCPJ(J)*grid%hbm2(I,J)
  1473. grid%fdiv(I,J)=FDIVJ(J)*grid%hbm2(I,J)
  1474. grid%fad(I,J)=FADJ(J)
  1475. grid%hdacv(I,J)=HDACJ(J)*grid%vbm2(I,J)
  1476. grid%hdac(I,J)=HDACJ(J)*1.25*grid%hbm2(I,J)
  1477. ENDDO
  1478. ENDDO
  1479. DO J=JTS, MIN(JDE-1,JTE)
  1480. IF (J.LE.5.OR.J.GE.(JDE-1)-4) THEN
  1481. KHH=(IDE-1)-2+MOD(J,2) ! KHH is global...loop over I that have
  1482. DO I=ITS,MIN(IDE-1,ITE)
  1483. IF (I .ge. 2 .and. I .le. KHH) THEN
  1484. grid%hdac(I,J)=grid%hdac(I,J)* DFC
  1485. ENDIF
  1486. ENDDO
  1487. ELSE
  1488. KHH=2+MOD(J,2)
  1489. DO I=ITS,MIN(IDE-1,ITE)
  1490. IF (I .ge. 2 .and. I .le. KHH) THEN
  1491. grid%hdac(I,J)=grid%hdac(I,J)* DFC
  1492. ENDIF
  1493. ENDDO
  1494. KHH=(IDE-1)-2+MOD(J,2)
  1495. DO I=ITS,MIN(IDE-1,ITE)
  1496. IF (I .ge. (IDE-1)-2 .and. I .le. KHH) THEN
  1497. grid%hdac(I,J)=grid%hdac(I,J)* DFC
  1498. ENDIF
  1499. ENDDO
  1500. ENDIF
  1501. ENDDO
  1502. DO J=JTS,min(JTE,JDE-1)
  1503. DO I=ITS,min(ITE,IDE-1)
  1504. grid%ddmpu(I,J)=DDMPUJ(J)*grid%vbm2(I,J)
  1505. grid%ddmpv(I,J)=DDMPVJ(J)*grid%vbm2(I,J)
  1506. grid%hdacv(I,J)=grid%hdacv(I,J)*grid%vbm2(I,J)
  1507. ENDDO
  1508. ENDDO
  1509. ! --------------INCREASING DIFFUSION ALONG THE BOUNDARIES----------------
  1510. DO J=JTS,MIN(JDE-1,JTE)
  1511. IF (J.LE.5.OR.J.GE.JDE-1-4) THEN
  1512. KVH=(IDE-1)-1-MOD(J,2)
  1513. DO I=ITS,min(IDE-1,ITE)
  1514. IF (I .ge. 2 .and. I .le. KVH) THEN
  1515. grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC
  1516. grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC
  1517. grid%hdacv(I,J)=grid%hdacv(I,J)* DFC
  1518. ENDIF
  1519. ENDDO
  1520. ELSE
  1521. KVH=3-MOD(J,2)
  1522. DO I=ITS,min(IDE-1,ITE)
  1523. IF (I .ge. 2 .and. I .le. KVH) THEN
  1524. grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC
  1525. grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC
  1526. grid%hdacv(I,J)=grid%hdacv(I,J)* DFC
  1527. ENDIF
  1528. ENDDO
  1529. KVH=(IDE-1)-1-MOD(J,2)
  1530. DO I=ITS,min(IDE-1,ITE)
  1531. IF (I .ge. IDE-1-2 .and. I .le. KVH) THEN
  1532. grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC
  1533. grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC
  1534. grid%hdacv(I,J)=grid%hdacv(I,J)* DFC
  1535. ENDIF
  1536. ENDDO
  1537. ENDIF
  1538. ENDDO
  1539. write(message,*) 'grid%stc(1)'
  1540. CALL wrf_message(message)
  1541. DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
  1542. write(message,635) (grid%stc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12+1)
  1543. CALL wrf_message(message)
  1544. ENDDO
  1545. write(message,*) 'grid%smc(1)'
  1546. CALL wrf_message(message)
  1547. DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
  1548. write(message,635) (grid%smc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12+1)
  1549. CALL wrf_message(message)
  1550. ENDDO
  1551. DO j = jts, MIN(jde-1,jte)
  1552. DO i= ITS, MIN(IDE-1,ITE)
  1553. if (grid%sm(I,J) .lt. 0.1 .and. grid%smc(I,1,J) .gt. 0.5 .and. grid%sice(I,J) .lt. 0.1) then
  1554. write(message,*) 'very moist on land point: ', I,J,grid%smc(I,1,J)
  1555. CALL wrf_debug(10,message)
  1556. endif
  1557. enddo
  1558. enddo
  1559. !!! compute grid%emt, grid%em on global domain, and only on task 0.
  1560. #ifdef DM_PARALLEL
  1561. IF (wrf_dm_on_monitor()) THEN !!!! NECESSARY TO LIMIT THIS TO TASK ZERO?
  1562. #else
  1563. IF (JDS .eq. JTS) THEN !! set unfailable condition for serial job
  1564. #endif
  1565. ALLOCATE(EMJ(JDS:JDE-1),EMTJ(JDS:JDE-1))
  1566. DO J=JDS,JDE-1
  1567. TPH=SB+float(J-1)*DPH
  1568. DXP=ERAD*DLM*COS(TPH)
  1569. EMJ(J)= GRID%DT/( 4.*DXP)*DTAD
  1570. EMTJ(J)=GRID%DT/(16.*DXP)*DTAD
  1571. ENDDO
  1572. JA=0
  1573. DO 161 J=3,5
  1574. JA=JA+1
  1575. KHLA(JA)=2
  1576. KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
  1577. 161 grid%emt(JA)=EMTJ(J)
  1578. DO 162 J=(JDE-1)-4,(JDE-1)-2
  1579. JA=JA+1
  1580. KHLA(JA)=2
  1581. KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
  1582. 162 grid%emt(JA)=EMTJ(J)
  1583. DO 163 J=6,(JDE-1)-5
  1584. JA=JA+1
  1585. KHLA(JA)=2
  1586. KHHA(JA)=2+MOD(J,2)
  1587. 163 grid%emt(JA)=EMTJ(J)
  1588. DO 164 J=6,(JDE-1)-5
  1589. JA=JA+1
  1590. KHLA(JA)=(IDE-1)-2
  1591. KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
  1592. 164 grid%emt(JA)=EMTJ(J)
  1593. ! --------------SPREADING OF UPSTREAM VELOCITY-POINT ADVECTION FACTOR----
  1594. JA=0
  1595. DO 171 J=3,5
  1596. JA=JA+1
  1597. KVLA(JA)=2
  1598. KVHA(JA)=(IDE-1)-1-MOD(J,2)
  1599. 171 grid%em(JA)=EMJ(J)
  1600. DO 172 J=(JDE-1)-4,(JDE-1)-2
  1601. JA=JA+1
  1602. KVLA(JA)=2
  1603. KVHA(JA)=(IDE-1)-1-MOD(J,2)
  1604. 172 grid%em(JA)=EMJ(J)
  1605. DO 173 J=6,(JDE-1)-5
  1606. JA=JA+1
  1607. KVLA(JA)=2
  1608. KVHA(JA)=2+MOD(J+1,2)
  1609. 173 grid%em(JA)=EMJ(J)
  1610. DO 174 J=6,(JDE-1)-5
  1611. JA=JA+1
  1612. KVLA(JA)=(IDE-1)-2
  1613. KVHA(JA)=(IDE-1)-1-MOD(J,2)
  1614. 174 grid%em(JA)=EMJ(J)
  1615. 696 continue
  1616. ENDIF ! wrf_dm_on_monitor/serial job
  1617. call NMM_SH2O(IMS,IME,JMS,JME,ITS,NNXP,JTS,NNYP,grid%num_soil_layers,grid%isltyp, &
  1618. grid%sm,grid%sice,grid%stc,grid%smc,grid%sh2o)
  1619. !! must be a better place to put this, but will eliminate "phantom"
  1620. !! wind points here (no wind point on eastern boundary of odd numbered rows)
  1621. IF ( abs(IDE-1-ITE) .eq. 1 ) THEN ! along eastern boundary
  1622. write(message,*) 'zero phantom winds'
  1623. CALL wrf_message(message)
  1624. DO K=1,KDE-1
  1625. DO J=JDS,JDE-1,2
  1626. IF (J .ge. JTS .and. J .le. JTE) THEN
  1627. grid%u(IDE-1,J,K)=0.
  1628. grid%v(IDE-1,J,K)=0.
  1629. ENDIF
  1630. ENDDO
  1631. ENDDO
  1632. ENDIF
  1633. 969 continue
  1634. DO j = jms, jme
  1635. DO i = ims, ime
  1636. fisx=max(grid%fis(i,j),0.)
  1637. grid%z0(I,J) =grid%sm(I,J)*Z0SEA+(1.-grid%sm(I,J))* &
  1638. & (0.*Z0MAX+FISx *FCM+Z0LAND)
  1639. ENDDO
  1640. ENDDO
  1641. write(message,*) 'grid%z0 over memory, leaving module_initialize_real'
  1642. CALL wrf_message(message)
  1643. DO J=JME,JMS,-((JME-JMS)/20+1)
  1644. write(message,635) (grid%z0(I,J),I=IMS,IME,(IME-IMS)/14+1)
  1645. CALL wrf_message(message)
  1646. ENDDO
  1647. endif ! on first_time check
  1648. write(message,*) 'leaving init_domain_nmm'
  1649. CALL wrf_message( TRIM(message) )
  1650. !
  1651. write(message,*)'STUFF MOVED TO REGISTRY:',grid%IDTAD, &
  1652. & grid%NSOIL,grid%NRADL,grid%NRADS,grid%NPHS,grid%NCNVC,grid%sigma
  1653. CALL wrf_message( TRIM(message) )
  1654. #ifdef HWRF
  1655. !=========================================================================================
  1656. ! gopal's doing for ocean coupling. Produce a high resolution grid for the entire domain
  1657. !=========================================================================================
  1658. if(internal_time_loop.eq.1) then !Kwon's doing
  1659. NDLMD=grid%dlmd/3.
  1660. NDPHD=grid%dphd/3.
  1661. NIDE=3*(IDE-1)-2
  1662. NJDE=3*(JDE-1)-2
  1663. ILOC=1
  1664. JLOC=1
  1665. NWBD= WBD ! + (ILOC -1)*2.*grid%dlmd + MOD(JLOC+1,2)*grid%dlmd
  1666. NSBD= SBD ! + (JLOC -1)*grid%dphd
  1667. ALLOCATE (NHLAT(NIDE,NJDE))
  1668. ALLOCATE (NHLON(NIDE,NJDE))
  1669. ALLOCATE (NVLAT(NIDE,NJDE))
  1670. ALLOCATE (NVLON(NIDE,NJDE))
  1671. ALLOCATE (HRES_SM(NIDE,NJDE))
  1672. #if defined(DM_PARALLEL)
  1673. if(wrf_dm_on_monitor()) then
  1674. ! Only the monitor process does the actual work (kinda
  1675. ! stupid; should be parallelized, but it's better than
  1676. ! writing garbage like it did before with >1 process)
  1677. ! Get high-res lat & lon:
  1678. CALL EARTH_LATLON_hwrf ( NHLAT,NHLON,NVLAT,NVLON, & ! rotated lat,lon at H and V points
  1679. NDLMD,NDPHD,NWBD,NSBD, &
  1680. tph0d,tlm0d, &
  1681. 1,NIDE,1,NJDE,1,1, &
  1682. 1,NIDE,1,NJDE,1,1, &
  1683. 1,NIDE,1,NJDE,1,1 )
  1684. ! Interpolate landmask to high-res grid:
  1685. CALL G2T2H_hwrf ( SM_G,HRES_SM, & ! output grid index and weights
  1686. NHLAT,NHLON, & ! target (hres) input lat lon in degrees
  1687. grid%DLMD,grid%DPHD,WBD,SBD, & ! parent res, west and south boundaries
  1688. tph0d,tlm0d, & ! parent central lat,lon, all in degrees
  1689. IDE,JDE,IDS,IDE,JDS,JDE, & ! parent imax and jmax, ime,jme
  1690. 1,NIDE,1,NJDE,1,1, &
  1691. 1,NIDE,1,NJDE,1,1, &
  1692. 1,NIDE,1,NJDE,1,1 )
  1693. ! We're done with the low-res sm grid now:
  1694. deallocate(SM_G)
  1695. ! Write out high-res grid for coupler:
  1696. WRITE(65)NHLAT(1:NIDE,1:NJDE)
  1697. WRITE(65)NHLON(1:NIDE,1:NJDE)
  1698. WRITE(65)NVLAT(1:NIDE,1:NJDE)
  1699. WRITE(65)NVLON(1:NIDE,1:NJDE)
  1700. WRITE(65)HRES_SM(1:NIDE,1:NJDE)
  1701. endif
  1702. #else
  1703. ! This code is the same as above, but for the non-mpi version:
  1704. CALL EARTH_LATLON_hwrf ( NHLAT,NHLON,NVLAT,NVLON, & ! rotated lat,lon at H and V points
  1705. NDLMD,NDPHD,NWBD,NSBD, &
  1706. tph0d,tlm0d, &
  1707. 1,NIDE,1,NJDE,1,1, &
  1708. 1,NIDE,1,NJDE,1,1, &
  1709. 1,NIDE,1,NJDE,1,1 )
  1710. CALL G2T2H_hwrf ( grid%SM,HRES_SM, & ! output grid index and weights
  1711. NHLAT,NHLON, & ! target (hres) input lat lon in degrees
  1712. grid%DLMD,grid%DPHD,WBD,SBD, & ! parent res, west and south boundaries
  1713. tph0d,tlm0d, & ! parent central lat,lon, all in degrees
  1714. IDE,JDE,IMS,IME,JMS,JME, & ! parent imax and jmax, ime,jme
  1715. 1,NIDE,1,NJDE,1,1, &
  1716. 1,NIDE,1,NJDE,1,1, &
  1717. 1,NIDE,1,NJDE,1,1 )
  1718. WRITE(65)NHLAT(1:NIDE,1:NJDE)
  1719. WRITE(65)NHLON(1:NIDE,1:NJDE)
  1720. WRITE(65)NVLAT(1:NIDE,1:NJDE)
  1721. WRITE(65)NVLON(1:NIDE,1:NJDE)
  1722. WRITE(65)HRES_SM(1:NIDE,1:NJDE)
  1723. #endif
  1724. DEALLOCATE (NHLAT)
  1725. DEALLOCATE (NHLON)
  1726. DEALLOCATE (NVLAT)
  1727. DEALLOCATE (NVLON)
  1728. DEALLOCATE (HRES_SM)
  1729. endif !Kwon's doing
  1730. !==================================================================================
  1731. ! end gopal's doing for ocean coupling.
  1732. !==================================================================================
  1733. #endif
  1734. !#define COPY_OUT
  1735. !#include <scalar_derefs.inc>
  1736. RETURN
  1737. END SUBROUTINE init_domain_nmm
  1738. !------------------------------------------------------
  1739. SUBROUTINE define_nmm_vertical_coord ( LM, PTSGM, pt, pdtop,HYBLEVS, &
  1740. SG1,DSG1,SGML1, &
  1741. SG2,DSG2,SGML2,dfl, dfrlg )
  1742. IMPLICIT NONE
  1743. ! USE module_model_constants
  1744. !!! certain physical parameters here probably don't need to be defined, as defined
  1745. !!! elsewhere within WRF. Done for initial testing purposes.
  1746. INTEGER :: LM, LPT2, L
  1747. REAL :: PTSGM, pt, PL, PT2, pdtop
  1748. REAL :: RGOG, PSIG,PHYB,PHYBM
  1749. REAL, PARAMETER :: Rd = 287.04 ! J deg{-1} kg{-1}
  1750. REAL, PARAMETER :: CP=1004.6,GAMMA=.0065,PRF0=101325.,T0=288.
  1751. REAL, PARAMETER :: g=9.81
  1752. REAL, DIMENSION(LM) :: DSG,DSG1,DSG2
  1753. REAL, DIMENSION(LM) :: SGML1,SGML2
  1754. REAL, DIMENSION(LM+1) :: SG1,SG2,HYBLEVS,dfl,dfrlg
  1755. CHARACTER(LEN=255) :: message
  1756. LPT2=LM+1
  1757. write(message,*) 'pt= ', pt
  1758. CALL wrf_message(message)
  1759. DO L=LM+1,1,-1
  1760. pl=HYBLEVS(L)*(101325.-pt)+pt
  1761. if(pl.lt.ptSGm) LPT2=l
  1762. ENDDO
  1763. IF(LPT2.lt.LM+1) THEN
  1764. pt2=HYBLEVS(LPT2)*(101325.-pt)+pt
  1765. ELSE
  1766. pt2=pt
  1767. ENDIF
  1768. write(message,*) '*** Sigma system starts at ',pt2,' Pa, from level ',LPT2
  1769. CALL wrf_message(message)
  1770. pdtop=pt2-pt
  1771. write(message,*) 'allocating DSG,DSG1,DSG2 as ', LM
  1772. CALL wrf_debug(10,message)
  1773. DSG=-99.
  1774. DO L=1,LM
  1775. DSG(L)=HYBLEVS(L)- HYBLEVS(L+1)
  1776. ENDDO
  1777. DSG1=0.
  1778. DSG2=0.
  1779. DO L=LM,1,-1
  1780. IF(L.ge.LPT2) then
  1781. DSG1(L)=DSG(L)
  1782. ELSE
  1783. DSG2(L)=DSG(L)
  1784. ENDIF
  1785. ENDDO
  1786. SGML1=-99.
  1787. SGML2=-99.
  1788. IF(LPT2.le.LM+1) THEN
  1789. DO L=LM+1,LPT2,-1
  1790. SG2(L)=0.
  1791. ENDDO
  1792. DO L=LPT2,2,-1
  1793. SG2(L-1)=SG2(L)+DSG2(L-1)
  1794. ENDDO
  1795. DO L=LPT2-1,1,-1
  1796. SG2(L)=SG2(L)/SG2(1)
  1797. ENDDO
  1798. SG2(1)=1.
  1799. DO L=LPT2-1,1,-1
  1800. DSG2(L)=SG2(L)-SG2(L+1)
  1801. SGML2(l)=(SG2(l)+SG2(l+1))*0.5
  1802. ENDDO
  1803. ENDIF
  1804. DO L=LM,LPT2,-1
  1805. DSG2(L)=0.
  1806. SGML2(L)=0.
  1807. ENDDO
  1808. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1809. SG1(LM+1)=0.
  1810. DO L=LM+1,LPT2,-1
  1811. SG1(L-1)=SG1(L)+DSG1(L-1)
  1812. ENDDO
  1813. DO L=LM,LPT2,-1
  1814. SG1(L)=SG1(L)/SG1(LPT2-1)
  1815. ENDDO
  1816. SG1(LPT2-1)=1.
  1817. do l=LPT2-2,1,-1
  1818. SG1(l)=1.
  1819. enddo
  1820. DO L=LM,LPT2,-1
  1821. DSG1(L)=SG1(L)-SG1(L+1)
  1822. SGML1(L)=(SG1(L)+SG1(L+1))*0.5
  1823. ENDDO
  1824. DO L=LPT2-1,1,-1
  1825. DSG1(L)=0.
  1826. SGML1(L)=1.
  1827. ENDDO
  1828. 1000 format('l,hyblevs,psig,SG1,SG2,phyb,phybm')
  1829. 1100 format(' ',i4,f7.4,f10.2,2f7.4,2f10.2)
  1830. write(message,1000)
  1831. CALL wrf_debug(100,message)
  1832. do l=1,LM+1
  1833. psig=HYBLEVS(L)*(101325.-pt)+pt
  1834. phyb=SG1(l)*pdtop+SG2(l)*(101325.-pdtop-pt)+pt
  1835. if(l.lt.LM+1) then
  1836. phybm=SGML1(l)*pdtop+SGML2(l)*(101325.-pdtop-pt)+pt
  1837. else
  1838. phybm=-99.
  1839. endif
  1840. write(message,1100) l,HYBLEVS(L),psig &
  1841. ,SG1(l),SG2(l),phyb,phybm
  1842. CALL wrf_debug(100,message)
  1843. enddo
  1844. 632 format(f9.6)
  1845. write(message,*) 'SG1'
  1846. CALL wrf_debug(100,message)
  1847. do L=LM+1,1,-1
  1848. write(message,632) SG1(L)
  1849. CALL wrf_debug(100,message)
  1850. enddo
  1851. write(message,*) 'SG2'
  1852. CALL wrf_debug(100,message)
  1853. do L=LM+1,1,-1
  1854. write(message,632) SG2(L)
  1855. CALL wrf_debug(100,message)
  1856. enddo
  1857. write(message,*) 'DSG1'
  1858. CALL wrf_debug(100,message)
  1859. do L=LM,1,-1
  1860. write(message,632) DSG1(L)
  1861. CALL wrf_debug(100,message)
  1862. enddo
  1863. write(message,*) 'DSG2'
  1864. CALL wrf_debug(100,message)
  1865. do L=LM,1,-1
  1866. write(message,632) DSG2(L)
  1867. CALL wrf_debug(100,message)
  1868. enddo
  1869. write(message,*) 'SGML1'
  1870. CALL wrf_debug(100,message)
  1871. do L=LM,1,-1
  1872. write(message,632) SGML1(L)
  1873. CALL wrf_debug(100,message)
  1874. enddo
  1875. write(message,*) 'SGML2'
  1876. CALL wrf_debug(100,message)
  1877. do L=LM,1,-1
  1878. write(message,632) SGML2(L)
  1879. CALL wrf_debug(100,message)
  1880. enddo
  1881. rgog=(rd*gamma)/g
  1882. DO L=1,LM+1
  1883. dfl(L)=g*T0*(1.-((pt+SG1(L)*pdtop+SG2(L)*(101325.-pt2)) &
  1884. /101325.)**rgog)/gamma
  1885. dfrlg(L)=dfl(L)/g
  1886. write(message,*) 'L, dfl(L): ', L, dfl(L)
  1887. CALL wrf_debug(100,message)
  1888. ENDDO
  1889. END SUBROUTINE define_nmm_vertical_coord
  1890. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1891. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1892. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1893. SUBROUTINE compute_nmm_surfacep ( TERRAIN_HGT_T, Z3D_IN, PRESS3D_IN, T3D_IN &
  1894. &, psfc_out,generic &
  1895. &, IDS,IDE,JDS,JDE,KDS,KDE &
  1896. &, IMS,IME,JMS,JME,KMS,KME &
  1897. &, ITS,ITE,JTS,JTE,KTS,KTE )
  1898. IMPLICIT NONE
  1899. real, allocatable:: dum2d(:,:),DUM2DB(:,:)
  1900. integer :: IDS,IDE,JDS,JDE,KDS,KDE
  1901. integer :: IMS,IME,JMS,JME,KMS,KME
  1902. integer :: ITS,ITE,JTS,JTE,KTS,KTE,Ilook,Jlook
  1903. integer :: I,J,II,generic,L,KINSERT,K,bot_lev,LL
  1904. integer :: IHE(JMS:JME),IHW(JMS:JME)
  1905. real :: TERRAIN_HGT_T(IMS:IME,JMS:JME)
  1906. real :: Z3D_IN(IMS:IME,JMS:JME,generic)
  1907. real :: T3D_IN(IMS:IME,JMS:JME,generic)
  1908. real :: PRESS3D_IN(IMS:IME,JMS:JME,generic)
  1909. real :: PSFC_IN(IMS:IME,JMS:JME),TOPO_IN(IMS:IME,JMS:JME)
  1910. real :: psfc_out(IMS:IME,JMS:JME),rincr(IMS:IME,JMS:JME)
  1911. real :: dif1,dif2,dif3,dif4,dlnpdz,BOT_INPUT_HGT,BOT_INPUT_PRESS,dpdz,rhs
  1912. real :: zin(generic),pin(generic)
  1913. character (len=255) :: message
  1914. logical :: DEFINED_PSFC(IMS:IME,JMS:JME), DEFINED_PSFCB(IMS:IME,JMS:JME)
  1915. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  1916. Ilook=25
  1917. Jlook=25
  1918. DO j = JMS, JME
  1919. IHE(J)=MOD(J+1,2)
  1920. IHW(J)=IHE(J)-1
  1921. ENDDO
  1922. DO J=JMS,JME
  1923. DO I=IMS,IME
  1924. DEFINED_PSFC(I,J)=.FALSE.
  1925. DEFINED_PSFCB(I,J)=.FALSE.
  1926. IF (PRESS3D_IN(I,J,1) .ne. 200100.) THEN
  1927. PSFC_IN(I,J)=PRESS3D_IN(I,J,1)
  1928. TOPO_IN(I,J)=Z3D_IN(I,J,1)
  1929. ELSE
  1930. PSFC_IN(I,J)=PRESS3D_IN(I,J,2)
  1931. TOPO_IN(I,J)=Z3D_IN(I,J,2)
  1932. ENDIF
  1933. ENDDO
  1934. ENDDO
  1935. ! input surface pressure smoothing over the ocean - still needed for NAM?
  1936. II_loop: do II=1,8
  1937. CYCLE II_loop
  1938. do J=JTS+1,min(JTE,JDE-1)-1
  1939. do I=ITS+1,min(ITE,IDE-1)-1
  1940. rincr(I,J)=0.
  1941. if (PSFC_IN(I,J) .gt. 100000. .and. &
  1942. PSFC_IN(I+IHE(J),J+1) .gt. 100000. .and. &
  1943. PSFC_IN(I+IHE(J),J-1) .gt. 100000. .and. &
  1944. PSFC_IN(I+IHW(J),J+1) .gt. 100000. .and. &
  1945. PSFC_IN(I+IHW(J),J-1) .gt. 100000. ) then
  1946. dif1=abs(PSFC_IN(I,J)-PSFC_IN(I+IHE(J),J+1))
  1947. dif2=abs(PSFC_IN(I,J)-PSFC_IN(I+IHE(J),J-1))
  1948. dif3=abs(PSFC_IN(I,J)-PSFC_IN(I+IHW(J),J+1))
  1949. dif4=abs(PSFC_IN(I,J)-PSFC_IN(I+IHW(J),J-1))
  1950. if (max(dif1,dif2,dif3,dif4) .lt. 200. .and. TOPO_IN(I,J).le. 0.5 .and. &
  1951. TOPO_IN(I+IHE(J),J+1) .le. 0.5 .and. &
  1952. TOPO_IN(I+IHW(J),J+1) .le. 0.5 .and. &
  1953. TOPO_IN(I+IHE(J),J-1) .le. 0.5 .and. &
  1954. TOPO_IN(I+IHW(J),J-1) .lt. 0.5) then
  1955. rincr(I,J)=0.125*( 4.*PSFC_IN(I,J)+ &
  1956. PSFC_IN(I+IHE(J),J+1)+PSFC_IN(I+IHE(J),J-1)+ &
  1957. PSFC_IN(I+IHW(J),J+1)+PSFC_IN(I+IHW(J),J-1) ) &
  1958. - PSFC_IN(I,J)
  1959. ! if (rincr(I,J) .ne. 0 .and. abs(rincr(I,J)) .gt. 20.) then
  1960. ! write(message,*) 'II, I,J,rincr: ', II, I,J,rincr(I,J)
  1961. ! CALL wrf_message(message)
  1962. ! endif
  1963. endif
  1964. endif
  1965. ENDDO
  1966. ENDDO
  1967. DO J=JTS+1,min(JTE,JDE-1)-1
  1968. DO I=ITS+1,min(ITE,IDE-1)-1
  1969. PSFC_IN(I,J)=PSFC_IN(I,J) + rincr(I,J)
  1970. ENDDO
  1971. ENDDO
  1972. ! write(message,*) ' -------------------------------------------------- '
  1973. ! CALL wrf_message(message)
  1974. end do II_loop
  1975. ALLOCATE(DUM2D(IMS:IME,JMS:JME))
  1976. DO J=JMS,JME
  1977. DO I=IMS,IME
  1978. DUM2D(I,J)=-9.
  1979. END DO
  1980. END DO
  1981. DO J=JTS,min(JTE,JDE-1)
  1982. I_loop: DO I=ITS,min(ITE,IDE-1)
  1983. IF (PSFC_IN(I,J) .lt. 0.1) THEN
  1984. write(message,*) 'QUITTING BECAUSE I,J, PSFC_IN: ', I,J,PSFC_IN(I,J)
  1985. call wrf_error_fatal(message)
  1986. ENDIF
  1987. BOT_INPUT_PRESS=PSFC_IN(I,J)
  1988. BOT_INPUT_HGT=TOPO_IN(I,J)
  1989. IF (I .eq. Ilook .AND. J .eq. Jlook) THEN
  1990. write(message,*) ' TERRAIN_HGT_T: ', I,J, TERRAIN_HGT_T(I,J)
  1991. CALL wrf_message(message)
  1992. write(message,*) ' PSFC_IN, TOPO_IN: ', &
  1993. I, J, PSFC_IN(I,J),TOPO_IN(I,J)
  1994. CALL wrf_message(message)
  1995. DO L=1,generic
  1996. write(message,*) ' L,PRESS3D_IN, Z3D_IN: ', &
  1997. I,J,L, PRESS3D_IN(I,J,L),Z3D_IN(I,J,L)
  1998. CALL wrf_debug(10,message)
  1999. END DO
  2000. ENDIF
  2001. DO L=2,generic-1
  2002. IF ( PRESS3D_IN(i,j,L) .gt. PSFC_IN(I,J) .AND. &
  2003. Z3D_IN(I,J,L) .lt. TERRAIN_HGT_T(I,J) .AND. &
  2004. Z3D_IN(I,J,L+1) .gt. TERRAIN_HGT_T(I,J) ) THEN
  2005. BOT_INPUT_PRESS=PRESS3D_IN(i,j,L)
  2006. BOT_INPUT_HGT=Z3D_IN(I,J,L)
  2007. ! IF (I .eq. Ilook .and. J .eq. Jlook) THEN
  2008. ! write(message,*) 'BOT_INPUT_PRESS, BOT_INPUT_HGT NOW : ', &
  2009. ! Ilook,Jlook, BOT_INPUT_PRESS, BOT_INPUT_HGT
  2010. ! CALL wrf_message(message)
  2011. ! ENDIF
  2012. ENDIF
  2013. END DO
  2014. !!!!!!!!!!!!!!!!!!!!!! START HYDRO CHECK
  2015. IF ( PRESS3D_IN(i,j,1) .ne. 200100. .AND. &
  2016. (PSFC_IN(I,J) .gt. PRESS3D_IN(i,j,2) .OR. &
  2017. TOPO_IN(I,J) .lt. Z3D_IN(I,J,2)) ) THEN ! extrapolate downward
  2018. IF (J .eq. JTS .AND. I .eq. ITS) THEN
  2019. write(message,*) 'hydro check - should only be for isobaric input'
  2020. CALL wrf_message(message)
  2021. ENDIF
  2022. IF (Z3D_IN(I,J,2) .ne. TOPO_IN(I,J)) THEN
  2023. dpdz=(PRESS3D_IN(i,j,2)-PSFC_IN(I,J))/(Z3D_IN(I,J,2)-TOPO_IN(I,J))
  2024. rhs=-9.81*((PRESS3D_IN(i,j,2)+ PSFC_IN(I,J))/2.)/(287.04* T3D_IN(I,J,2))
  2025. IF ( abs(PRESS3D_IN(i,j,2)-PSFC_IN(I,J)) .gt. 290.) THEN
  2026. IF (dpdz .lt. 1.05*rhs .OR. dpdz .gt. 0.95*rhs) THEN
  2027. write(message,*) 'I,J,P(2),Psfc,Z(2),Zsfc: ', &
  2028. I,J,PRESS3D_IN(i,j,2),PSFC_IN(I,J),Z3D_IN(I,J,2),TOPO_IN(I,J)
  2029. IF (mod(I,5).eq.0 .AND. mod(J,5).eq.0) CALL wrf_debug(50,message)
  2030. CYCLE I_loop
  2031. ENDIF
  2032. ENDIF
  2033. ELSE ! z(2) equals TOPO_IN
  2034. IF (PRESS3D_IN(i,j,2) .eq. PSFC_IN(I,J)) THEN
  2035. ! write(message,*) 'all equal at I,J: ', I,J
  2036. ! CALL wrf_message(message)
  2037. ELSE
  2038. ! write(message,*) 'heights equal, pressures not: ', &
  2039. ! PRESS3D_IN(i,j,2), PSFC_IN(I,J)
  2040. ! CALL wrf_message(message)
  2041. CYCLE I_loop
  2042. ENDIF
  2043. ENDIF
  2044. IF ( abs(PRESS3D_IN(i,j,2)-PSFC_IN(I,J)) .gt. 290.) THEN
  2045. IF (PRESS3D_IN(i,j,2) .lt. PSFC_IN(I,J) .and. &
  2046. Z3D_IN(I,J,2) .lt. TOPO_IN(I,J)) THEN
  2047. ! write(message,*) 'surface data mismatch(a) at I,J: ', I,J
  2048. ! CALL wrf_message(message)
  2049. CYCLE I_loop
  2050. ELSEIF (PRESS3D_IN(i,j,2) .gt. PSFC_IN(I,J) .AND. &
  2051. Z3D_IN(I,J,2) .gt. TOPO_IN(I,J)) THEN
  2052. ! write(message,*) 'surface data mismatch(b) at I,J: ', I,J
  2053. ! CALL wrf_message(message)
  2054. CYCLE I_loop
  2055. ENDIF
  2056. ENDIF
  2057. ENDIF
  2058. !!!!!!! loop over a few more levels
  2059. DO L=3,6
  2060. IF ( PRESS3D_IN(i,j,1) .ne. 200100. .AND. &
  2061. (((PSFC_IN(I,J)-PRESS3D_IN(i,j,L)) .lt. 400.) .OR. &
  2062. TOPO_IN(I,J) .lt. Z3D_IN(I,J,L))) then
  2063. IF (Z3D_IN(I,J,L) .ne. TOPO_IN(I,J)) THEN
  2064. dpdz=(PRESS3D_IN(i,j,L)-PSFC_IN(I,J))/ &
  2065. (Z3D_IN(I,J,L)-TOPO_IN(I,J))
  2066. rhs=-9.81*((PRESS3D_IN(i,j,L)+ PSFC_IN(I,J))/2.)/ &
  2067. (287.04*T3D_IN(I,J,L))
  2068. IF ( abs(PRESS3D_IN(i,j,L)-PSFC_IN(I,J)) .gt. 290.) THEN
  2069. IF (dpdz .lt. 1.05*rhs .or. dpdz .gt. 0.95*rhs) THEN
  2070. write(message,*) 'I,J,L,Piso,Psfc,Ziso,Zsfc: ', &
  2071. I,J,L,PRESS3D_IN(i,j,L),PSFC_IN(I,J),&
  2072. Z3D_IN(I,J,L),TOPO_IN(I,J)
  2073. IF (mod(I,5).eq.0 .AND. mod(J,5).eq.0) &
  2074. CALL wrf_debug(50,message)
  2075. CYCLE I_loop
  2076. ENDIF
  2077. ENDIF
  2078. ELSE
  2079. IF (PRESS3D_IN(i,j,2) .eq. PSFC_IN(I,J)) THEN
  2080. ! write(message,*) 'all equal at I,J: ', I,J
  2081. ! CALL wrf_message(message)
  2082. ELSE
  2083. CYCLE I_loop
  2084. ENDIF
  2085. ENDIF
  2086. ENDIF
  2087. IF ( abs(PRESS3D_IN(i,j,L)-PSFC_IN(I,J)) .gt. 290.) THEN
  2088. IF (PRESS3D_IN(i,j,L) .lt. PSFC_IN(I,J) .AND. &
  2089. Z3D_IN(I,J,L) .lt. TOPO_IN(I,J)) THEN
  2090. CYCLE I_loop
  2091. ELSEIF (PRESS3D_IN(i,j,L) .gt. PSFC_IN(I,J) .AND. &
  2092. Z3D_IN(I,J,L) .gt. TOPO_IN(I,J)) THEN
  2093. CYCLE I_loop
  2094. ENDIF
  2095. ENDIF
  2096. END DO
  2097. !!!!!!!!!!!!!!!!!!!!!! END HYDRO CHECK
  2098. IF (TERRAIN_HGT_T(I,J) .eq. BOT_INPUT_HGT ) THEN
  2099. dum2d(I,J)=BOT_INPUT_PRESS
  2100. IF (BOT_INPUT_HGT .ne. 0. .and. (BOT_INPUT_HGT-INT(BOT_INPUT_HGT) .ne. 0.) ) THEN
  2101. write(message,*) 'with BOT_INPUT_HGT: ', BOT_INPUT_HGT, &
  2102. 'set dum2d to bot_input_pres: ', I,J,dum2d(I,J)
  2103. CALL wrf_message(message)
  2104. ENDIF
  2105. IF (dum2d(I,J) .lt. 50000. .OR. dum2d(I,J) .gt. 109000.) THEN
  2106. write(message,*) 'bad dum2d(a): ', I,J,DUM2D(I,J)
  2107. CALL wrf_message(message)
  2108. ENDIF
  2109. ELSEIF (TERRAIN_HGT_T(I,J) .lt. BOT_INPUT_HGT ) THEN
  2110. ! target is below lowest possible input...extrapolate
  2111. IF ( BOT_INPUT_PRESS-PRESS3D_IN(I,J,2) .gt. 500. ) THEN
  2112. dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,2)) ) / &
  2113. (BOT_INPUT_HGT-Z3D_IN(i,j,2))
  2114. IF (I .eq. Ilook .and. J .eq. Jlook) THEN
  2115. write(message,*) 'I,J,dlnpdz(a): ', I,J,dlnpdz
  2116. CALL wrf_message(message)
  2117. ENDIF
  2118. ELSE
  2119. !! thin layer and/or just have lowest level - difference with 3rd level data
  2120. IF ( abs(BOT_INPUT_PRESS - PRESS3D_IN(i,j,3)) .gt. 290. ) THEN
  2121. dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,3)) ) / &
  2122. (BOT_INPUT_HGT-Z3D_IN(i,j,3))
  2123. IF (I .eq. Ilook .and. J .eq. Jlook) then
  2124. write(message,*) 'p diff: ', BOT_INPUT_PRESS, PRESS3D_IN(i,j,3)
  2125. CALL wrf_message(message)
  2126. write(message,*) 'z diff: ', BOT_INPUT_HGT, Z3D_IN(i,j,3)
  2127. CALL wrf_message(message)
  2128. ENDIF
  2129. ELSE
  2130. !! Loop up to level 7 looking for a sufficiently thick layer
  2131. FIND_THICK: DO LL=4,7
  2132. IF( abs(BOT_INPUT_PRESS - PRESS3D_IN(i,j,LL)) .gt. 290.) THEN
  2133. dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,LL)) ) / &
  2134. (BOT_INPUT_HGT-Z3D_IN(i,j,LL))
  2135. EXIT FIND_THICK
  2136. ENDIF
  2137. END DO FIND_THICK
  2138. ENDIF
  2139. ENDIF
  2140. dum2d(I,J)= exp(log(BOT_INPUT_PRESS) + dlnpdz * &
  2141. (TERRAIN_HGT_T(I,J) - BOT_INPUT_HGT) )
  2142. IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
  2143. write(message,*) 'bad dum2d(b): ', I,J,DUM2D(I,J)
  2144. CALL wrf_message(message)
  2145. write(message,*) 'BOT_INPUT_PRESS, dlnpdz, TERRAIN_HGT_T, BOT_INPUT_HGT: ', &
  2146. BOT_INPUT_PRESS, dlnpdz, TERRAIN_HGT_T(I,J), BOT_INPUT_HGT
  2147. CALL wrf_message(message)
  2148. write(message,*) 'Z3D_IN: ', Z3D_IN(I,J,1:10)
  2149. CALL wrf_message(message)
  2150. write(message,*) 'PRESS3D_IN: ', PRESS3D_IN(I,J,1:10)
  2151. CALL wrf_message(message)
  2152. ENDIF
  2153. ELSE ! target level bounded by input levels
  2154. DO L=2,generic-1
  2155. IF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,L) .AND. &
  2156. TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,L+1) ) THEN
  2157. dlnpdz= (log(PRESS3D_IN(i,j,l))-log(PRESS3D_IN(i,j,L+1)) ) / &
  2158. (Z3D_IN(i,j,l)-Z3D_IN(i,j,L+1))
  2159. dum2d(I,J)= log(PRESS3D_IN(i,j,l)) + &
  2160. dlnpdz * (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,L) )
  2161. dum2d(i,j)=exp(dum2d(i,j))
  2162. IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
  2163. write(message,*) 'bad dum2d(c): ', I,J,DUM2D(I,J)
  2164. CALL wrf_message(message)
  2165. ENDIF
  2166. ENDIF
  2167. ENDDO
  2168. !!! account for situation where BOT_INPUT_HGT < TERRAIN_HGT_T < Z3D_IN(:,2,:)
  2169. IF (dum2d(I,J) .eq. -9 .AND. BOT_INPUT_HGT .lt. TERRAIN_HGT_T(I,J) &
  2170. .AND. TERRAIN_HGT_T(I,J) .lt. Z3D_IN(I,J,2)) then
  2171. IF (mod(I,50) .eq. 0 .AND. mod(J,50) .eq. 0) THEN
  2172. write(message,*) 'I,J,BOT_INPUT_HGT, bot_pres, TERRAIN_HGT_T: ', &
  2173. I,J,BOT_INPUT_HGT, BOT_INPUT_PRESS, TERRAIN_HGT_T(I,J)
  2174. CALL wrf_message(message)
  2175. ENDIF
  2176. dlnpdz= (log(PSFC_IN(i,j))-log(PRESS3D_IN(i,j,2)) ) / &
  2177. (TOPO_IN(i,j)-Z3D_IN(i,j,2))
  2178. dum2d(I,J)= log(PSFC_IN(i,j)) + &
  2179. dlnpdz * (TERRAIN_HGT_T(I,J) - TOPO_IN(i,j) )
  2180. dum2d(i,j)= exp(dum2d(i,j))
  2181. IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
  2182. write(message,*) 'bad dum2d(d): ', I,J,DUM2D(I,J)
  2183. CALL wrf_message(message)
  2184. ENDIF
  2185. ENDIF
  2186. IF (dum2d(I,J) .eq. -9.) THEN
  2187. write(message,*) 'must have flukey situation in new ', I,J
  2188. CALL wrf_message(message)
  2189. write(message,*) 'I,J,BOT_INPUT_HGT, bot_pres, TERRAIN_HGT_T: ', &
  2190. I,J,BOT_INPUT_HGT, BOT_INPUT_PRESS, TERRAIN_HGT_T(I,J)
  2191. CALL wrf_message(message)
  2192. DO L=1,generic-1
  2193. IF ( TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,L) ) THEN
  2194. ! problematic with HGT_M substitution for "input" surface height?
  2195. dum2d(i,j)=PRESS3D_IN(I,J,L)
  2196. IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
  2197. write(message,*) 'bad dum2d(e): ', I,J,DUM2D(I,J)
  2198. CALL wrf_message(message)
  2199. ENDIF
  2200. ENDIF
  2201. ENDDO
  2202. IF ( TERRAIN_HGT_T(I,J) .eq. TOPO_IN(I,J)) THEN
  2203. dum2d(I,J)=PSFC_IN(I,J)
  2204. IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
  2205. write(message,*) 'bad dum2d(grid%f): ', I,J,DUM2D(I,J)
  2206. CALL wrf_message(message)
  2207. ENDIF
  2208. write(message,*) 'matched input topo, psfc: ', I,J,TOPO_IN(I,J),PSFC_IN(I,J)
  2209. CALL wrf_message(message)
  2210. ENDIF
  2211. IF (dum2d(I,J) .eq. -9.) THEN
  2212. CALL wrf_error_fatal("quitting due to undefined surface pressure")
  2213. ENDIF
  2214. ENDIF
  2215. DEFINED_PSFC(I,J)=.TRUE.
  2216. IF (I .eq. Ilook .AND. J .eq. Jlook) THEN
  2217. write(message,*) 'newstyle psfc: ', I,J,dum2d(I,J)
  2218. CALL wrf_message(message)
  2219. ENDIF
  2220. ENDIF
  2221. ENDDO I_loop
  2222. ENDDO
  2223. ! write(message,*) 'psfc points (new style)'
  2224. ! CALL wrf_message(message)
  2225. ! DO J=min(JTE,JDE-1),JTS,-loopinc
  2226. ! write(message,633) (dum2d(I,J)/100.,I=ITS,min(ITE,IDE-1),iloopinc)
  2227. ! CALL wrf_message(message)
  2228. ! END DO
  2229. 633 format(35(f5.0,1x))
  2230. write(message,*) 'PSFC extremes (new style)'
  2231. CALL wrf_message(message)
  2232. write(message,*) minval(dum2d,MASK=DEFINED_PSFC),maxval(dum2d,MASK=DEFINED_PSFC)
  2233. CALL wrf_message(message)
  2234. ! IF (minval(dum2d,MASK=DEFINED_PSFC) .lt. 50000. .or. maxval(dum2d,MASK=DEFINED_PSFC) .gt. 108000.) THEN
  2235. DO J=JTS,min(JTE,JDE-1)
  2236. DO I=ITS,min(ITE,IDE-1)
  2237. IF (DEFINED_PSFC(I,J) .AND. dum2d(I,J) .lt. 50000. ) THEN
  2238. IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN
  2239. WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
  2240. CALL wrf_debug(2,message)
  2241. ELSE
  2242. CALL wrf_error_fatal("quit due to unrealistic surface pressure")
  2243. ENDIF
  2244. ENDIF
  2245. IF (DEFINED_PSFC(I,J) .AND. dum2d(I,J) .gt. 108000. ) THEN
  2246. IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN
  2247. WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
  2248. CALL wrf_debug(2,message)
  2249. ELSE
  2250. CALL wrf_error_fatal("quit due to unrealistic surface pressure")
  2251. ENDIF
  2252. ENDIF
  2253. END DO
  2254. END DO
  2255. !! "traditional" isobaric only approach ------------------------------------------------
  2256. ALLOCATE (DUM2DB(IMS:IME,JMS:JME))
  2257. DO J=JMS,JME
  2258. DO I=IMS,IME
  2259. DUM2DB(I,J)=-9.
  2260. END DO
  2261. END DO
  2262. DO J=JTS,min(JTE,JDE-1)
  2263. DO I=ITS,min(ITE,IDE-1)
  2264. IF (TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,2)) THEN ! targ below lowest
  2265. IF ( abs(PRESS3D_IN(i,j,2)-PRESS3D_IN(i,j,3)) .gt. 290.) THEN
  2266. dlnpdz= (log(PRESS3D_IN(i,j,2))-log(PRESS3D_IN(i,j,3)) ) / &
  2267. (Z3D_IN(i,j,2)-Z3D_IN(i,j,3))
  2268. ELSE
  2269. dlnpdz= (log(PRESS3D_IN(i,j,2))-log(PRESS3D_IN(i,j,4)) ) / &
  2270. (Z3D_IN(i,j,2)-Z3D_IN(i,j,4))
  2271. ENDIF
  2272. DUM2DB(I,J)= exp( log(PRESS3D_IN(i,j,2)) + dlnpdz * &
  2273. (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,2)) )
  2274. IF (I .eq. Ilook .and. J .eq. Jlook) THEN
  2275. write(message,*) 'I,K, trad: dlnpdz, press_in(2), terrain_t, Z3D_IN(2): ', I,J,dlnpdz, &
  2276. PRESS3D_IN(i,j,2), TERRAIN_HGT_T(I,J), Z3D_IN(i,j,2)
  2277. CALL wrf_message(message)
  2278. ENDIF
  2279. DEFINED_PSFCB(i,j)=.true.
  2280. ELSEIF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,2)) THEN ! target level bounded by input levels
  2281. DO L=2,generic-1
  2282. IF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,L) .AND. &
  2283. TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,L+1) ) THEN
  2284. dlnpdz= (log(PRESS3D_IN(i,j,l))-log(PRESS3D_IN(i,j,L+1)) ) / &
  2285. (Z3D_IN(i,j,l)-Z3D_IN(i,j,L+1))
  2286. DUM2DB(I,J)= log(PRESS3D_IN(i,j,l)) + &
  2287. dlnpdz * (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,L) )
  2288. DUM2DB(i,j)=exp(DUM2DB(i,j))
  2289. DEFINED_PSFCB(i,j)=.true.
  2290. IF (DUM2DB(I,J) .lt. 13000.) THEN
  2291. write(message,*) 'I,J,L,terrain,Z3d(L),z3d(L+1),p3d(L),p3d(l+1): ', I,J,L, &
  2292. TERRAIN_HGT_T(I,J),Z3D_IN(I,J,L),Z3D_IN(I,J,L+1),PRESS3D_IN(I,J,L), &
  2293. PRESS3D_IN(I,J,L+1)
  2294. CALL wrf_error_fatal(message)
  2295. ENDIF
  2296. ENDIF
  2297. ENDDO
  2298. ELSEIF (TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,2)) THEN
  2299. DUM2DB(i,j)=PRESS3D_IN(I,J,2)
  2300. DEFINED_PSFCB(i,j)=.true.
  2301. ENDIF
  2302. IF (DUM2DB(I,J) .eq. -9.) THEN
  2303. write(message,*) 'must have flukey situation in trad ', I,J
  2304. CALL wrf_message(message)
  2305. DO L=1,generic-1
  2306. IF ( TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,L) ) THEN
  2307. DUM2DB(i,j)=PRESS3D_IN(I,J,L)
  2308. DEFINED_PSFCB(i,j)=.true.
  2309. ENDIF
  2310. ENDDO
  2311. ENDIF
  2312. IF (DUM2DB(I,J) .eq. -9.) THEN
  2313. write(message,*) 'HOPELESS PSFC, I QUIT'
  2314. CALL wrf_error_fatal(message)
  2315. ENDIF
  2316. if (I .eq. Ilook .and. J .eq. Jlook) THEN
  2317. write(message,*) ' traditional psfc: ', I,J,DUM2DB(I,J)
  2318. CALL wrf_message(message)
  2319. ENDIF
  2320. ENDDO
  2321. ENDDO
  2322. ! write(message,*) 'psfc points (traditional)'
  2323. ! CALL wrf_message(message)
  2324. ! DO J=min(JTE,JDE-1),JTS,-loopinc
  2325. ! write(message,633) (DUM2DB(I,J)/100.,I=its,min(ite,IDE-1),iloopinc)
  2326. ! CALL wrf_message(message)
  2327. ! ENDDO
  2328. write(message,*) 'PSFC extremes (traditional)'
  2329. CALL wrf_message(message)
  2330. write(message,*) minval(DUM2DB,MASK=DEFINED_PSFCB),maxval(DUM2DB,MASK=DEFINED_PSFCB)
  2331. CALL wrf_message(message)
  2332. DO J=JTS,min(JTE,JDE-1)
  2333. DO I=ITS,min(ITE,IDE-1)
  2334. IF (DEFINED_PSFCB(I,J) .AND. dum2db(I,J) .lt. 50000. ) THEN
  2335. IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN
  2336. WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
  2337. CALL wrf_debug(2,message)
  2338. ELSE
  2339. CALL wrf_error_fatal("quit due to unrealistic surface pressure")
  2340. ENDIF
  2341. ENDIF
  2342. IF (DEFINED_PSFCB(I,J) .AND. dum2db(I,J) .gt. 108000. ) THEN
  2343. IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN
  2344. WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
  2345. CALL wrf_debug(2,message)
  2346. ELSE
  2347. CALL wrf_error_fatal("quit due to unrealistic surface pressure")
  2348. ENDIF
  2349. ENDIF
  2350. ! IF (DEFINED_PSFCB(I,J) .AND. ( dum2db(I,J) .lt. 50000. .OR. dum2db(I,J) .gt. 108000. )) THEN
  2351. ! IF (TERRAIN_HGT_T(I,J) .gt. -10. .and. TERRAIN_HGT_T(I,J) .lt. 5000.) THEN
  2352. ! write(0,*) 'I,J, terrain_hgt_t, dum2db: ', I,J, terrain_hgt_t(I,J),dum2db(I,J)
  2353. ! CALL wrf_error_fatal("quit due to unrealistic surface pressure")
  2354. ! ELSE
  2355. ! WRITE(message,*) 'surface pressure allowed because surface height is extreme value of: ', TERRAIN_HGT_T(I,J)
  2356. ! CALL wrf_debug(2,message)
  2357. ! ENDIF
  2358. ! ENDIF
  2359. ENDDO
  2360. ENDDO
  2361. !!!!! end traditional
  2362. DO J=JTS,min(JTE,JDE-1)
  2363. DO I=ITS,min(ITE,IDE-1)
  2364. IF (DEFINED_PSFCB(I,J) .and. DEFINED_PSFC(I,J)) THEN
  2365. IF ( abs(dum2d(I,J)-DUM2DB(I,J)) .gt. 400.) THEN
  2366. write(message,*) 'BIG DIFF I,J, dum2d, DUM2DB: ', I,J,dum2d(I,J),DUM2DB(I,J)
  2367. CALL wrf_message(message)
  2368. ENDIF
  2369. !! do we have enough confidence in new style to give it more than 50% weight?
  2370. psfc_out(I,J)=0.5*(dum2d(I,J)+DUM2DB(I,J))
  2371. ELSEIF (DEFINED_PSFC(I,J)) THEN
  2372. psfc_out(I,J)=dum2d(I,J)
  2373. ELSEIF (DEFINED_PSFCB(I,J)) THEN
  2374. psfc_out(I,J)=DUM2DB(I,J)
  2375. ELSE
  2376. write(message,*) 'I,J,dum2d,DUM2DB: ', I,J,dum2d(I,J),DUM2DB(I,J)
  2377. CALL wrf_message(message)
  2378. write(message,*) 'I,J,DEFINED_PSFC(I,J),DEFINED_PSFCB(I,J): ', I,J,DEFINED_PSFC(I,J),DEFINED_PSFCB(I,J)
  2379. CALL wrf_message(message)
  2380. call wrf_error_fatal("psfc_out completely undefined")
  2381. ENDIF
  2382. IF (I .eq. Ilook .AND. J .eq. Jlook) THEN
  2383. write(message,*) ' combined psfc: ', I,J,psfc_out(I,J)
  2384. CALL wrf_message(message)
  2385. ENDIF
  2386. IF (psfc_out(I,J) .lt. 50000. ) THEN
  2387. IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN
  2388. WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
  2389. CALL wrf_debug(2,message)
  2390. ELSE
  2391. write(message,*) 'possibly bad combo on psfc_out: ', I,J, psfc_out(I,J)
  2392. CALL wrf_debug(2,message)
  2393. write(message,*) 'DEFINED_PSFC, dum2d: ', DEFINED_PSFC(I,J),dum2d(I,J)
  2394. CALL wrf_debug(2,message)
  2395. write(message,*) 'DEFINED_PSFCB, DUM2DB: ', DEFINED_PSFCB(I,J),DUM2DB(I,J)
  2396. CALL wrf_debug(2,message)
  2397. CALL wrf_error_fatal("quit due to unrealistic surface pressure")
  2398. ENDIF
  2399. ENDIF
  2400. IF (psfc_out(I,J) .gt. 108000. ) THEN
  2401. IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN
  2402. WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
  2403. CALL wrf_debug(2,message)
  2404. ELSE
  2405. write(message,*) 'possibly bad combo on psfc_out: ', I,J, psfc_out(I,J)
  2406. CALL wrf_debug(2,message)
  2407. write(message,*) 'DEFINED_PSFC, dum2d: ', DEFINED_PSFC(I,J),dum2d(I,J)
  2408. CALL wrf_debug(2,message)
  2409. write(message,*) 'DEFINED_PSFCB, DUM2DB: ', DEFINED_PSFCB(I,J),DUM2DB(I,J)
  2410. CALL wrf_debug(2,message)
  2411. CALL wrf_error_fatal("quit due to unrealistic surface pressure")
  2412. ENDIF
  2413. ENDIF
  2414. ENDDO
  2415. ENDDO
  2416. deallocate(dum2d,dum2db)
  2417. END SUBROUTINE compute_nmm_surfacep
  2418. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  2419. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  2420. SUBROUTINE compute_3d_pressure(psfc_out,SGML1,SGML2,pdtop,pt &
  2421. &, pd,p3d_out &
  2422. &, IDS,IDE,JDS,JDE,KDS,KDE &
  2423. &, IMS,IME,JMS,JME,KMS,KME &
  2424. &, ITS,ITE,JTS,JTE,KTS,KTE )
  2425. INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
  2426. INTEGER :: IMS,IME,JMS,JME,KMS,KME
  2427. INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE
  2428. REAL, INTENT(IN) :: psfc_out(IMS:IME,JMS:JME)
  2429. REAL, INTENT(IN) :: SGML1(KDE),SGML2(KDE),pdtop,pt
  2430. REAL, INTENT(OUT):: p3d_out(IMS:IME,JMS:JME,KDS:KDE-1)
  2431. REAL, INTENT(OUT):: pd(IMS:IME,JMS:JME)
  2432. CHARACTER (len=255) :: message
  2433. ! write(message,*) 'pdtop, pt, psfc_out(1,1): ', pdtop, pt, psfc_out(1,1)
  2434. ! CALL wrf_message(message)
  2435. DO J=JTS,min(JTE,JDE-1)
  2436. DO I=ITS,min(ITE,IDE-1)
  2437. pd(I,J)=psfc_out(I,J)-pdtop-pt
  2438. ENDDO
  2439. ENDDO
  2440. DO J=JTS,min(JTE,JDE-1)
  2441. DO K=KDS,KDE-1
  2442. DO I=ITS,min(ITE,IDE-1)
  2443. p3d_out(I,J,K)=pd(I,J)*SGML2(K)+pdtop*SGML1(K)+pt
  2444. IF (p3d_out(I,J,K) .ge. psfc_out(I,J) .or. p3d_out(I,J,K) .le. pt) THEN
  2445. write(message,*) 'I,K,J,p3d_out: ', I,K,J,p3d_out(I,J,K)
  2446. CALL wrf_error_fatal(message)
  2447. ENDIF
  2448. ENDDO
  2449. ENDDO
  2450. ENDDO
  2451. END SUBROUTINE compute_3d_pressure
  2452. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  2453. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  2454. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  2455. SUBROUTINE interp_press2press_lin(press_in,press_out, &
  2456. data_in, data_out,generic &
  2457. &, extrapolate,ignore_lowest,TFIELD &
  2458. &, IDS,IDE,JDS,JDE,KDS,KDE &
  2459. &, IMS,IME,JMS,JME,KMS,KME &
  2460. &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time )
  2461. ! Interpolates data from one set of pressure surfaces to
  2462. ! another set of pressures
  2463. INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
  2464. INTEGER :: IMS,IME,JMS,JME,KMS,KME
  2465. INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic
  2466. INTEGER :: internal_time
  2467. ! REAL, INTENT(IN) :: press_in(IMS:IME,generic,JMS:JME)
  2468. REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic)
  2469. REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1)
  2470. ! REAL, INTENT(IN) :: data_in(IMS:IME,generic,JMS:JME)
  2471. REAL, INTENT(IN) :: data_in(IMS:IME,JMS:JME,generic)
  2472. REAL, INTENT(OUT) :: data_out(IMS:IME,JMS:JME,KMS:KME)
  2473. LOGICAL, INTENT(IN) :: extrapolate, ignore_lowest, TFIELD
  2474. LOGICAL :: col_smooth
  2475. INTEGER :: i,j
  2476. INTEGER :: k,kk
  2477. REAL :: desired_press
  2478. REAL :: dvaldlnp,dlnp,tadiabat,tiso
  2479. REAL, PARAMETER :: ADIAFAC=9.81/1004.
  2480. REAL, PARAMETER :: TSTEXTRAPFAC=.0065
  2481. DO K=KMS,KME
  2482. DO J=JMS,JME
  2483. DO I=IMS,IME
  2484. DATA_OUT(I,J,K)=-99999.9
  2485. ENDDO
  2486. ENDDO
  2487. ENDDO
  2488. IF (ignore_lowest) then
  2489. LMIN=2
  2490. ELSE
  2491. LMIN=1
  2492. ENDIF
  2493. DO j = JTS, min(JTE,JDE-1)
  2494. test_i: DO i = ITS, min(ITE,IDE-1)
  2495. IF (internal_time_loop .gt. 1) THEN
  2496. IF (J .ne. JDS .and. J .ne. JDE-1 .and. &
  2497. I .ne. IDS .and. I .ne. IDE-1 ) THEN
  2498. !! not on boundary
  2499. CYCLE test_i
  2500. ENDIF
  2501. ENDIF
  2502. col_smooth=.false.
  2503. output_loop: DO k = KDS,KDE-1
  2504. desired_press = press_out(i,j,k)
  2505. if (K .gt. KDS) then
  2506. if (TFIELD .and. col_smooth .and. desired_press .le. press_in(i,j,LMIN) &
  2507. .and. press_out(i,j,k-1) .ge. press_in(i,j,LMIN)) then
  2508. MAX_SMOOTH=K
  2509. ! write(message,*) 'I,J, MAX_SMOOTH: ', I,J, MAX_SMOOTH
  2510. ! CALL wrf_debug(100,message)
  2511. endif
  2512. endif
  2513. ! keep track of where the extrapolation begins
  2514. IF (desired_press .GT. press_in(i,j,LMIN)) THEN
  2515. IF (TFIELD .and. K .eq. 1 .and. (desired_press - press_in(i,j,LMIN)) .gt. 3000.) then
  2516. col_smooth=.TRUE. ! due to large extrapolation distance
  2517. ENDIF
  2518. IF ((desired_press - press_in(i,j,LMIN)).LT. 50.) THEN ! 0.5 mb
  2519. data_out(i,j,k) = data_in(i,j,LMIN)
  2520. ELSE
  2521. IF (extrapolate) THEN
  2522. ! Extrapolate downward because desired P level is below
  2523. ! the lowest level in our input data. Extrapolate using simple
  2524. ! 1st derivative of value with respect to ln P for the bottom 2
  2525. ! input layers.
  2526. ! Add a check to make sure we are not using the gradient of
  2527. ! a very thin layer
  2528. if (TFIELD) then
  2529. tiso=0.5*(data_in(i,j,1)+data_in(i,j,2))
  2530. endif
  2531. IF ( (press_in(i,j,LMIN)-press_in(i,j,LMIN+1)) .GT. 500.) THEN ! likely isobaric data
  2532. dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+1))
  2533. dvaldlnp = (data_in(i,j,LMIN) - data_in(i,j,LMIN+1)) / dlnp
  2534. ELSE ! assume terrain following
  2535. dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+5))
  2536. dvaldlnp = (data_in(i,j,LMIN) - data_in(i,j,LMIN+5)) / dlnp
  2537. ENDIF
  2538. data_out(i,j,k) = data_in(i,j,LMIN) + dvaldlnp * &
  2539. ( log(desired_press)-log(press_in(i,j,LMIN)) )
  2540. if (TFIELD .and. data_out(i,j,k) .lt. tiso-0.2) then
  2541. ! restrict slope to -1K/10 hPa
  2542. dvaldlnp=max(dvaldlnp, -1.0/ &
  2543. log( press_in(i,j,LMIN) / &
  2544. ( press_in(i,j,LMIN)-1000.) ))
  2545. data_out(I,J,K)= data_in(i,j,LMIN) + dvaldlnp * &
  2546. ( log(desired_press)-log(press_in(i,j,LMIN)) )
  2547. elseif (TFIELD .and. data_out(i,j,k) .gt. tiso+0.2) then
  2548. ! restrict slope to +0.8K/10 hPa
  2549. dvaldlnp=min(dvaldlnp, 0.8/ &
  2550. log( press_in(i,j,LMIN) / &
  2551. ( press_in(i,j,LMIN)-1000.) ))
  2552. data_out(I,J,K)= data_in(i,j,LMIN) + dvaldlnp * &
  2553. ( log(desired_press)-log(press_in(i,j,LMIN)) )
  2554. endif
  2555. ELSE
  2556. data_out(i,j,k) = data_in(i,j,LMIN)
  2557. ENDIF
  2558. ENDIF
  2559. ELSE IF (desired_press .LT. press_in(i,j,generic)) THEN
  2560. IF ( (press_in(i,j,generic) - desired_press) .LT. 10.) THEN
  2561. data_out(i,j,k) = data_in(i,j,generic)
  2562. ELSE
  2563. IF (extrapolate) THEN
  2564. ! Extrapolate upward
  2565. IF ((press_in(i,j,generic-1)-press_in(i,j,generic)).GT.50.) THEN
  2566. dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-1))
  2567. dvaldlnp=(data_in(i,j,generic)-data_in(i,j,generic-1))/dlnp
  2568. ELSE
  2569. dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-2))
  2570. dvaldlnp=(data_in(i,j,generic)-data_in(i,j,generic-2))/dlnp
  2571. ENDIF
  2572. data_out(i,j,k) = data_in(i,j,generic) + &
  2573. dvaldlnp * (log(desired_press)-log(press_in(i,j,generic)))
  2574. ELSE
  2575. data_out(i,j,k) = data_in(i,j,generic)
  2576. ENDIF
  2577. ENDIF
  2578. ELSE
  2579. ! We can trap between two levels and linearly interpolate
  2580. input_loop: DO kk = LMIN, generic-1
  2581. IF (desired_press .EQ. press_in(i,j,kk) )THEN
  2582. data_out(i,j,k) = data_in(i,j,kk)
  2583. EXIT input_loop
  2584. ELSE IF ( (desired_press .LT. press_in(i,j,kk)) .AND. &
  2585. (desired_press .GT. press_in(i,j,kk+1)) ) THEN
  2586. ! do trapped in lnp
  2587. dlnp = log(press_in(i,j,kk)) - log(press_in(i,j,kk+1))
  2588. dvaldlnp = (data_in(i,j,kk)-data_in(i,j,kk+1))/dlnp
  2589. data_out(i,j,k) = data_in(i,j,kk+1)+ &
  2590. dvaldlnp*(log(desired_press)-log(press_in(i,j,kk+1)))
  2591. EXIT input_loop
  2592. ENDIF
  2593. ENDDO input_loop
  2594. ENDIF
  2595. ENDDO output_loop
  2596. if (col_smooth) then
  2597. do K=max(KDS,MAX_SMOOTH-4),MAX_SMOOTH+4
  2598. data_out(I,J,K)=0.5*(data_out(I,J,K)+data_out(I,J,K+1))
  2599. enddo
  2600. endif
  2601. ENDDO test_i
  2602. ENDDO
  2603. END SUBROUTINE interp_press2press_lin
  2604. SUBROUTINE wind_adjust(press_in,press_out, &
  2605. U_in, V_in,U_out,V_out &
  2606. &, generic,depth_replace &
  2607. &, IDS,IDE,JDS,JDE,KDS,KDE &
  2608. &, IMS,IME,JMS,JME,KMS,KME &
  2609. &, ITS,ITE,JTS,JTE,KTS,KTE )
  2610. INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
  2611. INTEGER :: IMS,IME,JMS,JME,KMS,KME
  2612. INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic
  2613. INTEGER :: MAXLIN,MAXLOUT
  2614. REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic)
  2615. REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1)
  2616. REAL, INTENT(IN) :: U_in(IMS:IME,JMS:JME,generic)
  2617. REAL, INTENT(IN) :: V_in(IMS:IME,JMS:JME,generic)
  2618. REAL, INTENT(INOUT) :: U_out(IMS:IME,KMS:KME,JMS:JME)
  2619. REAL, INTENT(INOUT) :: V_out(IMS:IME,KMS:KME,JMS:JME)
  2620. REAL :: p1d_in(generic)
  2621. REAL :: p1d_out(KDS:KDE-1)
  2622. DO j = JTS, min(JTE,JDE-1)
  2623. DO i = ITS, min(ITE,IDE-1)
  2624. ! IF (press_out(I,J,1) .lt. press_in(I,J,2)) then
  2625. IF( (press_in(I,J,2)-press_out(I,J,1)) .gt. 200.) then
  2626. U_out(I,1,J)=U_in(I,J,2)
  2627. V_out(I,1,J)=V_in(I,J,2)
  2628. INLOOP: DO L=2,generic
  2629. p1d_in(L)=-9999.
  2630. IF ( (press_in(I,J,2)-press_in(I,J,L)) .lt. depth_replace) THEN
  2631. p1d_in(L)=(press_in(I,J,2)-press_in(I,J,L))
  2632. MAXLIN=L
  2633. ELSE
  2634. p1d_in(L)=(press_in(I,J,2)-press_in(I,J,L))
  2635. EXIT INLOOP
  2636. ENDIF
  2637. END DO INLOOP
  2638. OUTLOOP: DO L=KDS,KDE-1
  2639. p1d_out(L)=-9999.
  2640. IF ( (press_out(I,J,1)-press_out(I,J,L)) .lt. depth_replace) THEN
  2641. p1d_out(L)=(press_out(I,J,1)-press_out(I,J,L))
  2642. MAXLOUT=L
  2643. ELSE
  2644. EXIT OUTLOOP
  2645. ENDIF
  2646. END DO OUTLOOP
  2647. DO L=1,MAXLOUT
  2648. ptarg=p1d_out(L)
  2649. FINDLOOP: DO LL=2,MAXLIN
  2650. if (p1d_in(LL) .lt. ptarg .and. p1d_in(LL+1) .gt. ptarg) then
  2651. dlnp=log(p1d_in(LL))-log(p1d_in(LL+1))
  2652. dudlnp=(U_in(I,J,LL)-U_in(I,J,LL+1))/dlnp
  2653. dvdlnp=(V_in(I,J,LL)-V_in(I,J,LL+1))/dlnp
  2654. U_out(I,L,J)=U_in(I,J,LL)+dudlnp*(log(ptarg)-log(p1d_in(LL)))
  2655. V_out(I,L,J)=V_in(I,J,LL)+dvdlnp*(log(ptarg)-log(p1d_in(LL)))
  2656. EXIT FINDLOOP
  2657. endif
  2658. END DO FINDLOOP
  2659. END DO ! MAXLOUT loop
  2660. ENDIF
  2661. ENDDO
  2662. ENDDO
  2663. END SUBROUTINE wind_adjust
  2664. !--------------------------------------------------------------------
  2665. SUBROUTINE interp_press2press_log(press_in,press_out, &
  2666. data_in, data_out, generic &
  2667. &, extrapolate,ignore_lowest &
  2668. &, IDS,IDE,JDS,JDE,KDS,KDE &
  2669. &, IMS,IME,JMS,JME,KMS,KME &
  2670. &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time )
  2671. ! Interpolates ln(data) from one set of pressure surfaces to
  2672. ! another set of pressures
  2673. INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
  2674. INTEGER :: IMS,IME,JMS,JME,KMS,KME
  2675. INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic
  2676. INTEGER :: internal_time
  2677. ! REAL, INTENT(IN) :: press_in(IMS:IME,generic,JMS:JME)
  2678. REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic)
  2679. REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1)
  2680. ! REAL, INTENT(IN) :: data_in(IMS:IME,generic,JMS:JME)
  2681. ! REAL, INTENT(IN) :: data_in(IMS:IME,JMS:JME,generic)
  2682. REAL :: data_in(IMS:IME,JMS:JME,generic)
  2683. REAL, INTENT(OUT) :: data_out(IMS:IME,JMS:JME,KMS:KME)
  2684. LOGICAL, INTENT(IN) :: extrapolate, ignore_lowest
  2685. INTEGER :: i,j
  2686. INTEGER :: k,kk
  2687. REAL :: desired_press
  2688. REAL :: dlnvaldlnp,dlnp
  2689. DO K=1,generic
  2690. DO j = JTS, min(JTE,JDE-1)
  2691. DO i = ITS, min(ITE,IDE-1)
  2692. DATA_IN(I,J,K)=max(DATA_in(I,J,K),1.e-13)
  2693. ENDDO
  2694. ENDDO
  2695. ENDDO
  2696. DO K=KMS,KME
  2697. DO J=JMS,JME
  2698. DO I=IMS,IME
  2699. DATA_OUT(I,J,K)=-99999.9
  2700. ENDDO
  2701. ENDDO
  2702. ENDDO
  2703. IF (ignore_lowest) then
  2704. LMIN=2
  2705. ELSE
  2706. LMIN=1
  2707. ENDIF
  2708. DO j = JTS, min(JTE,JDE-1)
  2709. test_i: DO i = ITS, min(ITE,IDE-1)
  2710. IF (internal_time .gt. 1) THEN
  2711. IF (J .ne. JDS .and. J .ne. JDE-1 .and. &
  2712. I .ne. IDS .and. I .ne. IDE-1 ) THEN
  2713. !! not on boundary
  2714. CYCLE test_i
  2715. ENDIF
  2716. ENDIF
  2717. output_loop: DO k = KDS,KDE-1
  2718. desired_press = press_out(i,j,k)
  2719. IF (desired_press .GT. press_in(i,j,LMIN)) THEN
  2720. IF ((desired_press - press_in(i,j,LMIN)).LT. 10.) THEN ! 0.1 mb
  2721. data_out(i,j,k) = data_in(i,j,LMIN)
  2722. ELSE
  2723. IF (extrapolate) THEN
  2724. ! Extrapolate downward because desired P level is below
  2725. ! the lowest level in our input data. Extrapolate using simple
  2726. ! 1st derivative of value with respect to ln P for the bottom 2
  2727. ! input layers.
  2728. ! Add a check to make sure we are not using the gradient of
  2729. ! a very thin layer
  2730. IF ( (press_in(i,j,LMIN)-press_in(i,j,LMIN+1)) .GT. 100.) THEN
  2731. dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+1))
  2732. dlnvaldlnp = ( log(data_in(i,j,LMIN)) - log(data_in(i,j,LMIN+1)) ) / dlnp
  2733. ELSE
  2734. dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+2))
  2735. dlnvaldlnp = (log(data_in(i,j,LMIN)) - log(data_in(i,j,LMIN+2))) / dlnp
  2736. ENDIF
  2737. data_out(i,j,k) = exp(log(data_in(i,j,LMIN)) + dlnvaldlnp * &
  2738. ( log(desired_press)-log(press_in(i,j,LMIN))))
  2739. ELSE
  2740. data_out(i,j,k) = data_in(i,j,LMIN)
  2741. ENDIF
  2742. ENDIF
  2743. ELSE IF (desired_press .LT. press_in(i,j,generic)) THEN
  2744. IF ( (press_in(i,j,generic) - desired_press) .LT. 10.) THEN
  2745. data_out(i,j,k) = data_in(i,j,generic)
  2746. ELSE
  2747. IF (extrapolate) THEN
  2748. ! Extrapolate upward
  2749. IF ((press_in(i,j,generic-1)-press_in(i,j,generic)).GT.50.) THEN
  2750. dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-1))
  2751. dlnvaldlnp=(log(data_in(i,j,generic))-log(data_in(i,j,generic-1)))/dlnp
  2752. ELSE
  2753. dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-2))
  2754. dlnvaldlnp=(log(data_in(i,j,generic))-log(data_in(i,j,generic-2)))/dlnp
  2755. ENDIF
  2756. data_out(i,j,k) = exp(log(data_in(i,j,generic)) + &
  2757. dlnvaldlnp * (log(desired_press)-log(press_in(i,j,generic))))
  2758. ELSE
  2759. data_out(i,j,k) = data_in(i,j,generic)
  2760. ENDIF
  2761. ENDIF
  2762. ELSE
  2763. ! We can trap between two levels and linearly interpolate
  2764. input_loop: DO kk = LMIN, generic-1
  2765. IF (desired_press .EQ. press_in(i,j,kk) )THEN
  2766. data_out(i,j,k) = data_in(i,j,kk)
  2767. EXIT input_loop
  2768. ELSE IF ( (desired_press .LT. press_in(i,j,kk)) .AND. &
  2769. (desired_press .GT. press_in(i,j,kk+1)) ) THEN
  2770. ! do trapped in lnp
  2771. dlnp = log(press_in(i,j,kk)) - log(press_in(i,j,kk+1))
  2772. dlnvaldlnp = (log(data_in(i,j,kk))-log(data_in(i,j,kk+1)))/dlnp
  2773. data_out(i,j,k) = exp(log(data_in(i,j,kk+1))+ &
  2774. dlnvaldlnp*(log(desired_press)-log(press_in(i,j,kk+1))))
  2775. EXIT input_loop
  2776. ENDIF
  2777. ENDDO input_loop
  2778. ENDIF
  2779. ENDDO output_loop
  2780. ENDDO test_i
  2781. ENDDO
  2782. END SUBROUTINE interp_press2press_log
  2783. !-------------------------------------------------------------------
  2784. SUBROUTINE rh_to_mxrat (rh, t, p, q , wrt_liquid , &
  2785. ids , ide , jds , jde , kds , kde , &
  2786. ims , ime , jms , jme , kms , kme , &
  2787. its , ite , jts , jte , kts , kte )
  2788. IMPLICIT NONE
  2789. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  2790. ims , ime , jms , jme , kms , kme , &
  2791. its , ite , jts , jte , kts , kte
  2792. LOGICAL , INTENT(IN) :: wrt_liquid
  2793. ! REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p , t
  2794. ! REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: rh
  2795. REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(IN) :: p , t
  2796. REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(INOUT) :: rh
  2797. REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(OUT) :: q
  2798. ! Local vars
  2799. INTEGER :: i , j , k
  2800. REAL :: ew , q1 , t1
  2801. REAL, PARAMETER :: T_REF = 0.0
  2802. REAL, PARAMETER :: MW_AIR = 28.966
  2803. REAL, PARAMETER :: MW_VAP = 18.0152
  2804. REAL, PARAMETER :: A0 = 6.107799961
  2805. REAL, PARAMETER :: A1 = 4.436518521e-01
  2806. REAL, PARAMETER :: A2 = 1.428945805e-02
  2807. REAL, PARAMETER :: A3 = 2.650648471e-04
  2808. REAL, PARAMETER :: A4 = 3.031240396e-06
  2809. REAL, PARAMETER :: A5 = 2.034080948e-08
  2810. REAL, PARAMETER :: A6 = 6.136820929e-11
  2811. REAL, PARAMETER :: ES0 = 6.1121
  2812. REAL, PARAMETER :: C1 = 9.09718
  2813. REAL, PARAMETER :: C2 = 3.56654
  2814. REAL, PARAMETER :: C3 = 0.876793
  2815. REAL, PARAMETER :: EIS = 6.1071
  2816. REAL :: RHS
  2817. REAL, PARAMETER :: TF = 273.16
  2818. REAL :: TK
  2819. REAL :: ES
  2820. REAL :: QS
  2821. REAL, PARAMETER :: EPS = 0.622
  2822. REAL, PARAMETER :: SVP1 = 0.6112
  2823. REAL, PARAMETER :: SVP2 = 17.67
  2824. REAL, PARAMETER :: SVP3 = 29.65
  2825. REAL, PARAMETER :: SVPT0 = 273.15
  2826. ! This subroutine computes mixing ratio (q, kg/kg) from basic variables
  2827. ! pressure (p, Pa), temperature (t, K) and relative humidity (rh, 1-100%).
  2828. ! The reference temperature (t_ref, C) is used to describe the temperature
  2829. ! at which the liquid and ice phase change occurs.
  2830. DO k = kts , kte
  2831. DO j = jts , MIN ( jde-1 , jte )
  2832. DO i = its , MIN (ide-1 , ite )
  2833. rh(i,j,k) = MIN ( MAX ( rh(i,j,k) , 1. ) , 100. )
  2834. END DO
  2835. END DO
  2836. END DO
  2837. IF ( wrt_liquid ) THEN
  2838. DO k = kts , kte
  2839. DO j = jts , MIN ( jde-1 , jte )
  2840. DO i = its , MIN (ide-1 , ite )
  2841. es=svp1*10.*EXP(svp2*(t(i,j,k)-svpt0)/(t(i,j,k)-svp3))
  2842. qs=eps*es/(p(i,j,k)/100.-es)
  2843. q(i,j,k)=MAX(.01*rh(i,j,k)*qs,0.0)
  2844. END DO
  2845. END DO
  2846. END DO
  2847. ELSE
  2848. DO k = kts , kte
  2849. DO j = jts , MIN ( jde-1 , jte )
  2850. DO i = its , MIN (ide-1 , ite )
  2851. t1 = t(i,j,k) - 273.16
  2852. ! Obviously dry.
  2853. IF ( t1 .lt. -200. ) THEN
  2854. q(i,j,k) = 0
  2855. ELSE
  2856. ! First compute the ambient vapor pressure of water
  2857. IF ( ( t1 .GE. t_ref ) .AND. ( t1 .GE. -47.) ) THEN ! liq phase ESLO
  2858. ew = a0 + t1 * (a1 + t1 * (a2 + t1 * (a3 + t1 * (a4 + t1 * (a5 + t1 * a6)))))
  2859. ELSE IF ( ( t1 .GE. t_ref ) .AND. ( t1 .LT. -47. ) ) then !liq phas poor ES
  2860. ew = es0 * exp(17.67 * t1 / ( t1 + 243.5))
  2861. ELSE
  2862. tk = t(i,j,k)
  2863. rhs = -c1 * (tf / tk - 1.) - c2 * alog10(tf / tk) + &
  2864. c3 * (1. - tk / tf) + alog10(eis)
  2865. ew = 10. ** rhs
  2866. END IF
  2867. ! Now sat vap pres obtained compute local vapor pressure
  2868. ew = MAX ( ew , 0. ) * rh(i,j,k) * 0.01
  2869. ! Now compute the specific humidity using the partial vapor
  2870. ! pressures of water vapor (ew) and dry air (p-ew). The
  2871. ! constants assume that the pressure is in hPa, so we divide
  2872. ! the pressures by 100.
  2873. q1 = mw_vap * ew
  2874. q1 = q1 / (q1 + mw_air * (p(i,j,k)/100. - ew))
  2875. q(i,j,k) = q1 / (1. - q1 )
  2876. END IF
  2877. END DO
  2878. END DO
  2879. END DO
  2880. END IF
  2881. END SUBROUTINE rh_to_mxrat
  2882. !--=------------------------------------------------------------------
  2883. SUBROUTINE boundary_smooth(h, landmask, grid, nsmth , nrow &
  2884. &, IDS,IDE,JDS,JDE,KDS,KDE &
  2885. &, IMS,IME,JMS,JME,KMS,KME &
  2886. &, ITS,ITE,JTS,JTE,KTS,KTE )
  2887. implicit none
  2888. TYPE (domain) :: grid
  2889. integer :: IDS,IDE,JDS,JDE,KDS,KDE
  2890. integer :: IMS,IME,JMS,JME,KMS,KME
  2891. integer :: ITS,ITE,JTS,JTE,KTS,KTE
  2892. integer:: ihw(JDS:JDE-1),ihe(JDS:JDE-1),nsmth,nrow
  2893. real:: h(IMS:IME,JMS:JME),landmask(IMS:IME,JMS:JME)
  2894. real :: h_old(IMS:IME,JMS:JME)
  2895. real :: hbms(IDS:IDE-1,JDS:JDE-1)
  2896. real :: hse(IDS:IDE-1,JDS:JDE-1)
  2897. real :: hne(IDS:IDE-1,JDS:JDE-1)
  2898. integer :: IPS,IPE,JPS,JPE,KPS,KPE
  2899. integer :: ihl, ihh, m2l, ibas,jmelin
  2900. integer :: I,J,KS,IOFFSET,JSTART,JEND
  2901. character (len=255) :: message
  2902. ips=its
  2903. ipe=ite
  2904. jps=jts
  2905. jpe=jte
  2906. kps=kts
  2907. kpe=kte
  2908. do j= JTS,min(JTE,JDE-1)
  2909. ihw(J)=-mod(J,2)
  2910. ihe(j)=ihw(J)+1
  2911. end do
  2912. do J=JTS,min(JTE,JDE-1)
  2913. do I=ITS,min(ITE,IDE-1)
  2914. hbms(I,J)=landmask(I,J)
  2915. enddo
  2916. enddo
  2917. jmelin=(JDE-1)-nrow+1
  2918. ibas=nrow/2
  2919. m2l=mod(nrow,2)
  2920. do j=jts,min(jte,jde-1)
  2921. ihl=ibas+mod(j,2)+m2l*mod(J+1,2)
  2922. ihh=(IDE-1)-ibas-m2l*mod(J+1,2)
  2923. do i=its,min(ite,ide-1)
  2924. if (I .ge. ihl .and. I .le. ihh .and. J .ge. nrow .and. J .le. jmelin) then
  2925. hbms(I,J)=0.
  2926. endif
  2927. end do
  2928. end do
  2929. 634 format(30(f2.0,1x))
  2930. do KS=1,nsmth
  2931. grid%ht_gc=h
  2932. #ifdef DM_PARALLEL
  2933. # include "HALO_NMM_MG.inc"
  2934. #endif
  2935. h=grid%ht_gc
  2936. h_old=grid%ht_gc
  2937. do J=JTS,min(JTE,JDE-1)
  2938. do I=ITS, min(ITE,IDE-1)
  2939. if (I .ge. (IDS+mod(J,2)) .and. J .gt. JDS .and. J .lt. JDE-1 .and. I .lt. IDE-1) then
  2940. h(i,j)= ( h_old(i+ihe(j),j+1) + h_old(i+ihw(j),j-1) + h_old(i+ihe(j),j-1) + h_old(i+ihw(j),j+1) - &
  2941. 4. *h_old(i,j) )*hbms(i,j)*0.125+h_old(i,j)
  2942. endif
  2943. enddo
  2944. enddo
  2945. ! special treatment for four corners
  2946. if (hbms(1,1) .eq. 1 .and. ITS .le. 1 .and. JTS .le. 1) then
  2947. h(1,1)=0.75*h(1,1)+0.125*h(1+ihe(1),2)+ &
  2948. 0.0625*(h(2,1)+h(1,3))
  2949. endif
  2950. if (hbms(IDE-1,1) .eq. 1 .and. ITE .ge. IDE-2 .and. JTS .le. 1) then
  2951. h(IDE-1,1)=0.75*h(IDE-1,1)+0.125*h(IDE-1+ihw(1),2)+ &
  2952. 0.0625*(h(IDE-1-1,1)+h(IDE-1,3))
  2953. endif
  2954. if (hbms(1,JDE-1) .eq. 1 .and. ITS .le. 1 .and. JTE .ge. JDE-2) then
  2955. h(1,JDE-1)=0.75*h(1,JDE-1)+0.125*h(1+ihe(JDE-1),JDE-1-1)+ &
  2956. 0.0625*(h(2,JDE-1)+h(1,JDE-1-2))
  2957. endif
  2958. if (hbms(IDE-1,JDE-1) .eq. 1 .and. ITE .ge. IDE-2 .and. JTE .ge. JDE-2) then
  2959. h(IDE-1,JDE-1)=0.75*h(IDE-1,JDE-1)+0.125*h(IDE-1+ihw(JDE-1),JDE-1-1)+ &
  2960. 0.0625*(h(IDE-1-1,JDE-1)+h(IDE-1,JDE-1-2))
  2961. endif
  2962. do J=JMS,JME
  2963. do I=IMS,IME
  2964. grid%ht_gc(I,J)=h(I,J)
  2965. enddo
  2966. enddo
  2967. #ifdef DM_PARALLEL
  2968. # include "HALO_NMM_MG.inc"
  2969. #endif
  2970. do J=JMS,JME
  2971. do I=IMS,IME
  2972. h(I,J)=grid%ht_gc(I,J)
  2973. enddo
  2974. enddo
  2975. ! S bound
  2976. if (JTS .eq. JDS) then
  2977. J=JTS
  2978. do I=ITS,ITE
  2979. if (I .ge. IDS+1 .and. I .le. IDE-2) then
  2980. if (hbms(I,J) .eq. 1) then
  2981. h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J+1)+h(I+ihe(J),J+1))
  2982. endif
  2983. endif
  2984. enddo
  2985. endif
  2986. ! N bound
  2987. if (JTE .eq. JDE) then
  2988. J=JDE-1
  2989. write(message,*) 'DOING N BOUND SMOOTHING for J= ', J
  2990. CALL wrf_debug(100,message)
  2991. do I=ITS,min(ITE,IDE-1)
  2992. if (hbms(I,J) .eq. 1 .and. I .ge. IDS+1 .and. I .le. IDE-2) then
  2993. h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J-1)+h(I+ihe(J),J-1))
  2994. endif
  2995. enddo
  2996. endif
  2997. ! W bound
  2998. if (ITS .eq. IDS) then
  2999. I=ITS
  3000. do J=JTS,min(JTE,JDE-1)
  3001. if (hbms(I,J) .eq. 1 .and. J .ge. JDS+2 .and. J .le. JDE-3 .and. mod(J,2) .eq. 1) then
  3002. h(I,J)=0.75*h(I,J)+0.125*(h(I+ihe(J),J+1)+h(I+ihe(J),J-1))
  3003. endif
  3004. enddo
  3005. endif
  3006. ! E bound
  3007. if (ITE .eq. IDE) then
  3008. write(message,*) 'DOING E BOUND SMOOTHING for I= ', min(ITE,IDE-1)
  3009. CALL wrf_debug(100,message)
  3010. I=min(ITE,IDE-1)
  3011. do J=JTS,min(JTE,JDE-1)
  3012. if (hbms(I,J) .eq. 1 .and. J .ge. JDS+2 .and. J .le. JDE-3 .and. mod(J,2) .eq. 1) then
  3013. h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J+1)+h(I+ihw(J),J-1))
  3014. endif
  3015. enddo
  3016. endif
  3017. enddo ! end ks loop
  3018. do J=JMS,JME
  3019. do I=IMS,IME
  3020. grid%ht_gc(I,J)=h(I,J)
  3021. enddo
  3022. enddo
  3023. #ifdef DM_PARALLEL
  3024. #include "HALO_NMM_MG.inc"
  3025. #endif
  3026. do J=JMS,JME
  3027. do I=IMS,IME
  3028. h(I,J)=grid%ht_gc(I,J)
  3029. enddo
  3030. enddo
  3031. ! extra smoothing along inner boundary
  3032. if (JTS .eq. JDS) then
  3033. if (ITE .eq. IDE) then
  3034. IOFFSET=1
  3035. else
  3036. IOFFSET=0
  3037. endif
  3038. ! Southern Boundary
  3039. do i=its,min(ITE,IDE-1)-IOFFSET
  3040. h(i,2)=0.25*(h(i,1)+h(i+1,1)+ &
  3041. h(i,3)+h(i+1,3))
  3042. enddo
  3043. endif
  3044. if (JTE .eq. JDE) then
  3045. if (ITE .eq. IDE) then
  3046. IOFFSET=1
  3047. else
  3048. IOFFSET=0
  3049. endif
  3050. do i=its,min(ITE,IDE-1)-IOFFSET
  3051. h(i,(JDE-1)-1)=0.25*(h(i,(JDE-1)-2)+h(i+1,(JDE-1)-2)+ &
  3052. h(i,JDE-1)+h(i+1,JDE-1))
  3053. enddo
  3054. endif
  3055. if (JTS .eq. 1) then
  3056. JSTART=4
  3057. else
  3058. JSTART=JTS+mod(JTS,2) ! needs to be even
  3059. endif
  3060. if (JTE .eq. JDE) then
  3061. JEND=(JDE-1)-3
  3062. else
  3063. JEND=JTE
  3064. endif
  3065. if (ITS .eq. IDS) then
  3066. ! Western Boundary
  3067. do j=JSTART,JEND,2
  3068. h(1,j)=0.25*(h(1,j-1)+h(2,j-1)+ &
  3069. h(1,j+1)+h(2,j+1))
  3070. enddo
  3071. endif
  3072. if (ITE .eq. IDE) then
  3073. ! Eastern Boundary
  3074. do j=JSTART,JEND,2
  3075. h((IDE-1)-1,j)=0.25*(h((IDE-1)-1,j-1)+h((IDE-1),j-1)+ &
  3076. h((IDE-1)-1,j+1)+h((IDE-1),j+1))
  3077. enddo
  3078. endif
  3079. END SUBROUTINE boundary_smooth
  3080. !--------------------------------------------------------------------
  3081. SUBROUTINE monthly_interp_to_date ( field_in , date_str , field_out , &
  3082. ids , ide , jds , jde , kds , kde , &
  3083. ims , ime , jms , jme , kms , kme , &
  3084. its , ite , jts , jte , kts , kte )
  3085. ! Linrarly in time interpolate data to a current valid time. The data is
  3086. ! assumed to come in "monthly", valid at the 15th of every month.
  3087. IMPLICIT NONE
  3088. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  3089. ims , ime , jms , jme , kms , kme , &
  3090. its , ite , jts , jte , kts , kte
  3091. CHARACTER (LEN=24) , INTENT(IN) :: date_str
  3092. REAL , DIMENSION(ims:ime,jms:jme,12) , INTENT(IN) :: field_in
  3093. REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_out
  3094. ! Local vars
  3095. INTEGER :: i , j , l
  3096. INTEGER , DIMENSION(0:13) :: middle
  3097. INTEGER :: target_julyr , target_julday , target_date
  3098. INTEGER :: julyr , julday , int_month, next_month
  3099. REAL :: gmt
  3100. CHARACTER (LEN=4) :: yr
  3101. CHARACTER (LEN=2) :: mon , day15
  3102. WRITE(day15,FMT='(I2.2)') 15
  3103. DO l = 1 , 12
  3104. WRITE(mon,FMT='(I2.2)') l
  3105. CALL get_julgmt ( date_str(1:4)//'-'//mon//'-'//day15//'_'//'00:00:00.0000' , julyr , julday , gmt )
  3106. middle(l) = julyr*1000 + julday
  3107. END DO
  3108. l = 0
  3109. middle(l) = middle( 1) - 31
  3110. l = 13
  3111. middle(l) = middle(12) + 31
  3112. CALL get_julgmt ( date_str , target_julyr , target_julday , gmt )
  3113. target_date = target_julyr * 1000 + target_julday
  3114. find_month : DO l = 0 , 12
  3115. IF ( ( middle(l) .LT. target_date ) .AND. ( middle(l+1) .GE. target_date ) ) THEN
  3116. DO j = jts , MIN ( jde-1 , jte )
  3117. DO i = its , MIN (ide-1 , ite )
  3118. int_month = MOD ( l , 12 )
  3119. IF ( int_month .EQ. 0 ) int_month = 12
  3120. IF (int_month == 12) THEN
  3121. next_month=1
  3122. ELSE
  3123. next_month=int_month+1
  3124. ENDIF
  3125. field_out(i,j) = ( field_in(i,j,next_month) * ( target_date - middle(l) ) + &
  3126. field_in(i,j,int_month ) * ( middle(l+1) - target_date ) ) / &
  3127. ( middle(l+1) - middle(l) )
  3128. END DO
  3129. END DO
  3130. EXIT find_month
  3131. END IF
  3132. END DO find_month
  3133. END SUBROUTINE monthly_interp_to_date
  3134. !---------------------------------------------------------------------
  3135. SUBROUTINE monthly_min_max ( field_in , field_min , field_max , &
  3136. ids , ide , jds , jde , kds , kde , &
  3137. ims , ime , jms , jme , kms , kme , &
  3138. its , ite , jts , jte , kts , kte )
  3139. ! Plow through each month, find the max, min values for each i,j.
  3140. IMPLICIT NONE
  3141. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  3142. ims , ime , jms , jme , kms , kme , &
  3143. its , ite , jts , jte , kts , kte
  3144. REAL , DIMENSION(ims:ime,jms:jme,12) , INTENT(IN) :: field_in
  3145. REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_min , field_max
  3146. ! Local vars
  3147. INTEGER :: i , j , l
  3148. REAL :: minner , maxxer
  3149. DO j = jts , MIN(jde-1,jte)
  3150. DO i = its , MIN(ide-1,ite)
  3151. minner = field_in(i,j,1)
  3152. maxxer = field_in(i,j,1)
  3153. DO l = 2 , 12
  3154. IF ( field_in(i,j,l) .LT. minner ) THEN
  3155. minner = field_in(i,j,l)
  3156. END IF
  3157. IF ( field_in(i,j,l) .GT. maxxer ) THEN
  3158. maxxer = field_in(i,j,l)
  3159. END IF
  3160. END DO
  3161. field_min(i,j) = minner
  3162. field_max(i,j) = maxxer
  3163. END DO
  3164. END DO
  3165. END SUBROUTINE monthly_min_max
  3166. !-----------------------------------------------------------------------
  3167. SUBROUTINE reverse_vert_coord ( field, start_z, end_z &
  3168. &, IDS,IDE,JDS,JDE,KDS,KDE &
  3169. &, IMS,IME,JMS,JME,KMS,KME &
  3170. &, ITS,ITE,JTS,JTE,KTS,KTE )
  3171. IMPLICIT NONE
  3172. INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
  3173. ims , ime , jms , jme , kms , kme , &
  3174. its , ite , jts , jte , kts , kte, &
  3175. start_z, end_z
  3176. REAL, INTENT(INOUT) :: field(IMS:IME,JMS:JME,end_z)
  3177. ! local
  3178. INTEGER :: I,J,L
  3179. REAL, ALLOCATABLE :: dum3d(:,:,:)
  3180. allocate(dum3d(IMS:IME,JMS:JME,end_z))
  3181. DO L=start_z,end_z
  3182. DO J=jts,min(jte,jde-1)
  3183. DO I=its,min(ite,ide-1)
  3184. dum3d(I,J,L)=field(I,J,end_z-L+start_z)
  3185. END DO
  3186. END DO
  3187. END DO
  3188. DO L=start_z,end_z
  3189. DO J=jts,min(jte,jde-1)
  3190. DO I=its,min(ite,ide-1)
  3191. field(I,J,L)=dum3d(I,J,L)
  3192. END DO
  3193. END DO
  3194. END DO
  3195. DEALLOCATE(dum3d)
  3196. END SUBROUTINE reverse_vert_coord
  3197. !--------------------------------------------------------------------
  3198. SUBROUTINE compute_nmm_levels(ninterface, ptop, eta_levels)
  3199. USE module_model_constants
  3200. IMPLICIT NONE
  3201. character(len=132):: message
  3202. integer :: ninterface,Lthick,L
  3203. real, parameter:: gamma=.0065
  3204. real, parameter:: t_stand=288.
  3205. real, parameter:: p_stand=101325.
  3206. real :: maxdz_compute, ptop
  3207. real :: plower,pupper,tlay, sum
  3208. real :: eta_levels(ninterface)
  3209. real, allocatable:: Z(:)
  3210. real, allocatable:: deta_levels_spline(:)
  3211. logical:: print_pbl_warn
  3212. !----------------------------------------------------
  3213. allocate(Z(ninterface))
  3214. allocate(deta_levels_spline(ninterface-1))
  3215. CALL compute_eta_spline(ninterface-1,deta_levels_spline,ptop)
  3216. sum=0.
  3217. DO L=1,ninterface-1
  3218. sum=sum+deta_levels_spline(L)
  3219. ENDDO
  3220. eta_levels(1)=1.00
  3221. DO L=2,ninterface
  3222. eta_levels(L)=eta_levels(L-1)-deta_levels_spline(L-1)
  3223. ENDDO
  3224. eta_levels(ninterface)=0.00
  3225. DO L=2,ninterface-1
  3226. eta_levels(L)=0.5*(eta_levels(L))+0.25*(eta_levels(L-1)+eta_levels(L+1))
  3227. ENDDO
  3228. Z(1)=0.
  3229. maxdz_compute=0.
  3230. print_pbl_warn=.false.
  3231. DO L=2,ninterface
  3232. tlay=max( t_stand-gamma*Z(L-1), 216.5)
  3233. plower=ptop+(p_stand-ptop)*eta_levels(L-1)
  3234. pupper=ptop+(p_stand-ptop)*eta_levels(L)
  3235. Z(L)=Z(L-1)+(tlay*r_d/g)*(log(plower)-log(pupper))
  3236. if (plower .gt. 85000. .and. pupper .lt. 85000. .and. L .lt. 10) then
  3237. print_pbl_warn=.true.
  3238. endif
  3239. write(message,*) 'L, eta(l), pupper, Z(L): ', L, eta_levels(L),pupper,Z(L)
  3240. CALL wrf_debug(100,message)
  3241. if (Z(L)-Z(L-1) .gt. maxdz_compute) then
  3242. Lthick=L
  3243. endif
  3244. maxdz_compute=max(maxdz_compute,Z(L)-Z(L-1))
  3245. ENDDO
  3246. if (print_pbl_warn) then
  3247. write(message,*) 'WARNING - PBL MAY BE POORLY RESOLVED WITH NUMBER OF VERTICAL LEVELS'
  3248. CALL wrf_message(message)
  3249. write(message,*) ' - CONSIDER INCREASING THE VERTICAL RESOLUTION'
  3250. CALL wrf_message(message)
  3251. endif
  3252. write(message,*) 'thickest layer was: ', maxdz_compute , 'meters thick at level: ', Lthick
  3253. CALL wrf_message(message)
  3254. END SUBROUTINE compute_nmm_levels
  3255. !---------------------------
  3256. SUBROUTINE compute_eta_spline(LM, dsg, ptop)
  3257. IMPLICIT NONE
  3258. real:: dsg(LM), ptop, sum, rsum
  3259. real, allocatable:: xold(:),dold(:)
  3260. real, allocatable:: xnew(:),sgm(:)
  3261. real, allocatable:: pps(:),qqs(:),y2s(:)
  3262. integer nlev,LM,L,KOLD
  3263. IF (LM .ge. 46) THEN
  3264. KOLD=9
  3265. allocate(xold(KOLD))
  3266. allocate(dold(KOLD))
  3267. xold(1)=.00
  3268. dold(1)=.006
  3269. xold(2)=.13
  3270. dold(2)=.009
  3271. xold(3)=.19
  3272. dold(3)=.012
  3273. xold(4)=.30
  3274. dold(4)=.036
  3275. xold(5)=.42
  3276. dold(5)=.041
  3277. xold(6)=.56
  3278. dold(6)=.040
  3279. xold(7)=.69
  3280. dold(7)=.018
  3281. if (ptop .ge. 2000.) then
  3282. xold(8)=.90
  3283. dold(8)=.012
  3284. xold(9)=1.0
  3285. dold(9)=.006
  3286. else
  3287. xold(8)=.90
  3288. dold(8)=.008
  3289. xold(9)=1.0
  3290. dold(9)=.003
  3291. endif
  3292. ELSE
  3293. KOLD=8
  3294. allocate(xold(KOLD))
  3295. allocate(dold(KOLD))
  3296. xold(1)=.00
  3297. dold(1)=.006
  3298. xold(2)=.18
  3299. dold(2)=.015
  3300. xold(3)=.32
  3301. dold(3)=.035
  3302. xold(4)=.50
  3303. dold(4)=.040
  3304. xold(5)=.68
  3305. dold(5)=.030
  3306. xold(6)=.75
  3307. dold(6)=.017
  3308. xold(7)=.85
  3309. dold(7)=.012
  3310. if (ptop .ge. 2000.) then
  3311. xold(8)=1.0
  3312. dold(8)=.013
  3313. else
  3314. xold(8)=1.0
  3315. dold(8)=.008
  3316. endif
  3317. ENDIF
  3318. allocate(xnew(lm))
  3319. allocate(sgm(lm+1))
  3320. allocate(pps(lm))
  3321. allocate(qqs(lm))
  3322. allocate(y2s(lm))
  3323. DO L=1,LM
  3324. xnew(l)=float(l-1)/float(lm-1)
  3325. ENDDO
  3326. y2s=0.
  3327. CALL spline(kold,xold,dold,y2s,lm,xnew,dsg,pps,qqs)
  3328. sum=0.
  3329. DO l=1,lm
  3330. sum=sum+dsg(l)
  3331. ENDDO
  3332. rsum=1./sum
  3333. sgm(1)=0.
  3334. DO L=1,lm-1
  3335. dsg(l)=dsg(l)*rsum
  3336. sgm(l+1)=sgm(l)+dsg(l)
  3337. ENDDO
  3338. sgm(lm+1)=1.
  3339. dsg(lm)=sgm(lm+1)-sgm(lm)
  3340. END SUBROUTINE compute_eta_spline
  3341. ! -------------------------------------------------------------------
  3342. subroutine spline(NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,q)
  3343. ! ********************************************************************
  3344. ! * *
  3345. ! * THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE *
  3346. ! * PROGRAMED FOR A SMALL SCALAR MACHINE. *
  3347. ! * *
  3348. ! * PROGRAMER Z. JANJIC *
  3349. ! * *
  3350. ! * NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3. *
  3351. ! * XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
  3352. ! * FUNCTION ARE GIVEN. MUST BE IN ASCENDING ORDER. *
  3353. ! * YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. *
  3354. ! * Y2 - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL *
  3355. ! * SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE *
  3356. ! * SPECIFIED. *
  3357. ! * NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED. *
  3358. ! * XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
  3359. ! * FUNCTION ARE CALCULATED. XNEW(K) MUST BE GE XOLD(1) *
  3360. ! * AND LE XOLD(NOLD). *
  3361. ! * YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. *
  3362. ! * P, q - AUXILIARY VECTORS OF THE LENGTH NOLD-2. *
  3363. ! * *
  3364. ! ********************************************************************
  3365. !
  3366. ! LOG:
  3367. !
  3368. ! JOVIC - July 2008 - fixed incorrectly dimensioned arrays,
  3369. ! PYLE and do loop leading to out of bound array
  3370. ! reference
  3371. !------
  3372. !
  3373. ! PYLE - June 2007 - eliminated use of GO TO statements.
  3374. !
  3375. !-----------------------------------------------------------------------
  3376. IMPLICIT NONE
  3377. !-----------------------------------------------------------------------
  3378. INTEGER,INTENT(IN) :: NNEW,NOLD
  3379. REAL,DIMENSION(NOLD),INTENT(IN) :: XOLD,YOLD
  3380. REAL,DIMENSION(NNEW),INTENT(IN) :: XNEW
  3381. REAL,DIMENSION(NNEW),INTENT(OUT) :: YNEW
  3382. REAL,DIMENSION(NOLD+2),INTENT(INOUT) :: P,q,Y2
  3383. !
  3384. INTEGER :: K,K1,K2,KOLD,NOLDM1, K2_hold, K_hold
  3385. REAL :: AK,BK,CK,DEN,DX,DXC,DXL,DXR,DYDXL,DYDXR &
  3386. & ,RDX,RTDXC,X,XK,XSQ,Y2K,Y2KP1
  3387. !-----------------------------------------------------------------------
  3388. NOLDM1=NOLD-1
  3389. DXL=XOLD(2)-XOLD(1)
  3390. DXR=XOLD(3)-XOLD(2)
  3391. DYDXL=(YOLD(2)-YOLD(1))/DXL
  3392. DYDXR=(YOLD(3)-YOLD(2))/DXR
  3393. RTDXC=0.5/(DXL+DXR)
  3394. P(1)= RTDXC*(6.*(DYDXR-DYDXL)-DXL*Y2(1))
  3395. q(1)=-RTDXC*DXR
  3396. K=3
  3397. first_loop: DO K=3,NOLD-1
  3398. DXL=DXR
  3399. DYDXL=DYDXR
  3400. DXR=XOLD(K+1)-XOLD(K)
  3401. DYDXR=(YOLD(K+1)-YOLD(K))/DXR
  3402. DXC=DXL+DXR
  3403. DEN=1./(DXL*q(K-2)+DXC+DXC)
  3404. P(K-1)= DEN*(6.*(DYDXR-DYDXL)-DXL*P(K-2))
  3405. q(K-1)=-DEN*DXR
  3406. END DO first_loop
  3407. DO K=NOLDM1,2,-1
  3408. Y2(K)=P(K-1)+q(K-1)*Y2(K+1)
  3409. K_hold=K
  3410. END DO
  3411. K=K_hold
  3412. !-----------------------------------------------------------------------
  3413. second_loop: DO K1=1,NNEW
  3414. XK=XNEW(K1)
  3415. third_loop: DO K2=2,NOLD
  3416. IF(XOLD(K2)>XK)THEN
  3417. KOLD=K2-1
  3418. K2_hold=K2
  3419. exit third_loop
  3420. ENDIF
  3421. K2_hold=K2
  3422. END DO third_loop
  3423. IF (XOLD(K2_hold) .le. XK) THEN
  3424. YNEW(K1)=YOLD(NOLD)
  3425. CYCLE second_loop
  3426. ENDIF
  3427. IF (K1 .eq. 1 .or. K .ne. KOLD) THEN
  3428. K=KOLD
  3429. Y2K=Y2(K)
  3430. Y2KP1=Y2(K+1)
  3431. DX=XOLD(K+1)-XOLD(K)
  3432. RDX=1./DX
  3433. AK=.1666667*RDX*(Y2KP1-Y2K)
  3434. BK=0.5*Y2K
  3435. CK=RDX*(YOLD(K+1)-YOLD(K))-.1666667*DX*(Y2KP1+Y2K+Y2K)
  3436. ENDIF
  3437. X=XK-XOLD(K)
  3438. XSQ=X*X
  3439. YNEW(K1)=AK*XSQ*X+BK*XSQ+CK*X+YOLD(K)
  3440. END DO second_loop
  3441. END SUBROUTINE SPLINE
  3442. !--------------------------------------------------------------------
  3443. SUBROUTINE NMM_SH2O(IMS,IME,JMS,JME,ISTART,IM,JSTART,JM,&
  3444. NSOIL,ISLTPK, &
  3445. sm,sice,stc,smc,sh2o)
  3446. !! INTEGER, PARAMETER:: NSOTYP=9
  3447. ! INTEGER, PARAMETER:: NSOTYP=16
  3448. INTEGER, PARAMETER:: NSOTYP=19 !!!!!!!!MAYBE???
  3449. REAL :: PSIS(NSOTYP),BETA(NSOTYP),SMCMAX(NSOTYP)
  3450. REAL :: stc(IMS:IME,NSOIL,JMS:JME), &
  3451. smc(IMS:IME,NSOIL,JMS:JME)
  3452. REAL :: sh2o(IMS:IME,NSOIL,JMS:JME),sice(IMS:IME,JMS:JME),&
  3453. sm(IMS:IME,JMS:JME)
  3454. REAL :: HLICE,GRAV,T0,BLIM
  3455. INTEGER :: ISLTPK(IMS:IME,JMS:JME)
  3456. CHARACTER(LEN=255) :: message
  3457. ! Constants used in cold start sh2o initialization
  3458. DATA HLICE/3.335E5/,GRAV/9.81/,T0/273.15/
  3459. DATA BLIM/5.5/
  3460. ! DATA PSIS /0.04,0.62,0.47,0.14,0.10,0.26,0.14,0.36,0.04/
  3461. ! DATA BETA /4.26,8.72,11.55,4.74,10.73,8.17,6.77,5.25,4.26/
  3462. ! DATA SMCMAX /0.421,0.464,0.468,0.434,0.406, &
  3463. ! 0.465,0.404,0.439,0.421/
  3464. !!! NOT SURE...PSIS=SATPSI, BETA=BB??
  3465. DATA PSIS /0.069, 0.036, 0.141, 0.759, 0.759, 0.355, &
  3466. 0.135, 0.617, 0.263, 0.098, 0.324, 0.468, &
  3467. 0.355, 0.000, 0.069, 0.036, 0.468, 0.069, 0.069 /
  3468. DATA BETA/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, &
  3469. 6.66, 8.72, 8.17, 10.73, 10.39, 11.55, &
  3470. 5.25, 0.00, 2.79, 4.26, 11.55, 2.79, 2.79 /
  3471. DATA SMCMAX/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
  3472. 0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
  3473. 0.439, 1.000, 0.200, 0.421, 0.468, 0.200, 0.339/
  3474. DO K=1,NSOIL
  3475. DO J=JSTART,JM
  3476. DO I=ISTART,IM
  3477. !tst
  3478. IF (smc(I,K,J) .gt. SMCMAX(ISLTPK(I,J))) then
  3479. if (K .eq. 1) then
  3480. write(message,*) 'I,J,reducing smc from ' ,I,J,smc(I,K,J), 'to ', SMCMAX(ISLTPK(I,J))
  3481. CALL wrf_debug(100,message)
  3482. endif
  3483. smc(I,K,J)=SMCMAX(ISLTPK(I,J))
  3484. ENDIF
  3485. !tst
  3486. IF ( (sm(I,J) .lt. 0.5) .and. (sice(I,J) .lt. 0.5) ) THEN
  3487. IF (ISLTPK(I,J) .gt. 19) THEN
  3488. WRITE(message,*) 'FORCING ISLTPK at : ', I,J
  3489. CALL wrf_message(message)
  3490. ISLTPK(I,J)=9
  3491. ELSEIF (ISLTPK(I,J) .le. 0) then
  3492. WRITE(message,*) 'FORCING ISLTPK at : ', I,J
  3493. CALL wrf_message(message)
  3494. ISLTPK(I,J)=1
  3495. ENDIF
  3496. ! cold start: determine liquid soil water content (sh2o)
  3497. ! sh2o <= smc for t < 273.149K (-0.001C)
  3498. IF (stc(I,K,J) .LT. 273.149) THEN
  3499. ! first guess following explicit solution for Flerchinger Eqn from Koren
  3500. ! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
  3501. BX = BETA(ISLTPK(I,J))
  3502. IF ( BETA(ISLTPK(I,J)) .GT. BLIM ) BX = BLIM
  3503. if ( GRAV*(-PSIS(ISLTPK(I,J))) .eq. 0 ) then
  3504. write(message,*) 'TROUBLE'
  3505. CALL wrf_message(message)
  3506. write(message,*) 'I,J: ', i,J
  3507. CALL wrf_message(message)
  3508. write(message,*) 'grav, isltpk, psis(isltpk): ', grav,isltpk(I,J),&
  3509. psis(isltpk(I,J))
  3510. CALL wrf_message(message)
  3511. endif
  3512. if (BX .eq. 0 .or. stc(I,K,J) .eq. 0) then
  3513. write(message,*) 'TROUBLE -- I,J,BX, stc: ', I,J,BX,stc(I,K,J)
  3514. CALL wrf_message(message)
  3515. endif
  3516. FK = (((HLICE/(GRAV*(-PSIS(ISLTPK(I,J)))))* &
  3517. ((stc(I,K,J)-T0)/stc(I,K,J)))** &
  3518. (-1/BX))*SMCMAX(ISLTPK(I,J))
  3519. IF (FK .LT. 0.02) FK = 0.02
  3520. sh2o(I,K,J) = MIN ( FK, smc(I,K,J) )
  3521. ! ----------------------------------------------------------------------
  3522. ! now use iterative solution for liquid soil water content using
  3523. ! FUNCTION FRH2O (from the Eta "NOAH" land-surface model) with the
  3524. ! initial guess for sh2o from above explicit first guess.
  3525. sh2o(I,K,J)=FRH2O_init(stc(I,K,J),smc(I,K,J),sh2o(I,K,J), &
  3526. SMCMAX(ISLTPK(I,J)),BETA(ISLTPK(I,J)), &
  3527. PSIS(ISLTPK(I,J)))
  3528. ELSE ! above freezing
  3529. sh2o(I,K,J)=smc(I,K,J)
  3530. ENDIF
  3531. ELSE ! water point
  3532. sh2o(I,K,J)=smc(I,K,J)
  3533. ENDIF ! test on land/ice/sea
  3534. if (sh2o(I,K,J) .gt. SMCMAX(ISLTPK(I,J))) then
  3535. write(message,*) 'sh2o > THAN SMCMAX ', I,J,sh2o(I,K,J),SMCMAX(ISLTPK(I,J)),smc(I,K,J)
  3536. CALL wrf_message(message)
  3537. endif
  3538. ENDDO
  3539. ENDDO
  3540. ENDDO
  3541. END SUBROUTINE NMM_SH2O
  3542. !-------------------------------------------------------------------
  3543. FUNCTION FRH2O_init(TKELV,smc,sh2o,SMCMAX,B,PSIS)
  3544. IMPLICIT NONE
  3545. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3546. ! PURPOSE: CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT
  3547. ! IF TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION
  3548. ! TO SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF
  3549. ! KOREN ET AL. (1999, JGR, VOL 104(D16), 19569-19585).
  3550. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3551. !
  3552. ! New version (JUNE 2001): much faster and more accurate newton iteration
  3553. ! achieved by first taking log of eqn cited above -- less than 4
  3554. ! (typically 1 or 2) iterations achieves convergence. Also, explicit
  3555. ! 1-step solution option for special case of parameter Ck=0, which reduces
  3556. ! the original implicit equation to a simpler explicit form, known as the
  3557. ! ""Flerchinger Eqn". Improved handling of solution in the limit of
  3558. ! freezing point temperature T0.
  3559. !
  3560. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3561. !
  3562. ! INPUT:
  3563. !
  3564. ! TKELV.........Temperature (Kelvin)
  3565. ! smc...........Total soil moisture content (volumetric)
  3566. ! sh2o..........Liquid soil moisture content (volumetric)
  3567. ! SMCMAX........Saturation soil moisture content (from REDPRM)
  3568. ! B.............Soil type "B" parameter (from REDPRM)
  3569. ! PSIS..........Saturated soil matric potential (from REDPRM)
  3570. !
  3571. ! OUTPUT:
  3572. ! FRH2O.........supercooled liquid water content.
  3573. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3574. REAL B
  3575. REAL BLIM
  3576. REAL BX
  3577. REAL CK
  3578. REAL DENOM
  3579. REAL DF
  3580. REAL DH2O
  3581. REAL DICE
  3582. REAL DSWL
  3583. REAL ERROR
  3584. REAL FK
  3585. REAL FRH2O_init
  3586. REAL GS
  3587. REAL HLICE
  3588. REAL PSIS
  3589. REAL sh2o
  3590. REAL smc
  3591. REAL SMCMAX
  3592. REAL SWL
  3593. REAL SWLK
  3594. REAL TKELV
  3595. REAL T0
  3596. INTEGER NLOG
  3597. INTEGER KCOUNT
  3598. PARAMETER (CK=8.0)
  3599. ! PARAMETER (CK=0.0)
  3600. PARAMETER (BLIM=5.5)
  3601. ! PARAMETER (BLIM=7.0)
  3602. PARAMETER (ERROR=0.005)
  3603. PARAMETER (HLICE=3.335E5)
  3604. PARAMETER (GS = 9.81)
  3605. PARAMETER (DICE=920.0)
  3606. PARAMETER (DH2O=1000.0)
  3607. PARAMETER (T0=273.15)
  3608. ! ### LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM) ####
  3609. ! ### SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT ####
  3610. ! ### IS NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES ####
  3611. ! ################################################################
  3612. !
  3613. BX = B
  3614. IF ( B .GT. BLIM ) BX = BLIM
  3615. ! ------------------------------------------------------------------
  3616. ! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
  3617. NLOG=0
  3618. KCOUNT=0
  3619. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3620. ! C IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), sh2o = smc
  3621. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3622. IF (TKELV .GT. (T0 - 1.E-3)) THEN
  3623. FRH2O_init=smc
  3624. ELSE
  3625. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3626. IF (CK .NE. 0.0) THEN
  3627. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3628. ! CCCCCCCCC OPTION 1: ITERATED SOLUTION FOR NONZERO CK CCCCCCCCCCC
  3629. ! CCCCCCCCCCCC IN KOREN ET AL, JGR, 1999, EQN 17 CCCCCCCCCCCCCCCCC
  3630. ! INITIAL GUESS FOR SWL (frozen content)
  3631. SWL = smc-sh2o
  3632. ! KEEP WITHIN BOUNDS.
  3633. IF (SWL .GT. (smc-0.02)) SWL=smc-0.02
  3634. IF(SWL .LT. 0.) SWL=0.
  3635. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3636. ! C START OF ITERATIONS
  3637. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3638. DO WHILE (NLOG .LT. 10 .AND. KCOUNT .EQ. 0)
  3639. NLOG = NLOG+1
  3640. DF = ALOG(( PSIS*GS/HLICE ) * ( ( 1.+CK*SWL )**2. ) * &
  3641. ( SMCMAX/(smc-SWL) )**BX) - ALOG(-(TKELV-T0)/TKELV)
  3642. DENOM = 2. * CK / ( 1.+CK*SWL ) + BX / ( smc - SWL )
  3643. SWLK = SWL - DF/DENOM
  3644. ! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
  3645. IF (SWLK .GT. (smc-0.02)) SWLK = smc - 0.02
  3646. IF(SWLK .LT. 0.) SWLK = 0.
  3647. ! MATHEMATICAL SOLUTION BOUNDS APPLIED.
  3648. DSWL=ABS(SWLK-SWL)
  3649. SWL=SWLK
  3650. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3651. ! CC IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
  3652. ! CC WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
  3653. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3654. IF ( DSWL .LE. ERROR ) THEN
  3655. KCOUNT=KCOUNT+1
  3656. END IF
  3657. END DO
  3658. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3659. ! C END OF ITERATIONS
  3660. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3661. ! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
  3662. FRH2O_init = smc - SWL
  3663. ! CCCCCCCCCCCCCCCCCCCCCCCC END OPTION 1 CCCCCCCCCCCCCCCCCCCCCCCCCCC
  3664. ENDIF
  3665. IF (KCOUNT .EQ. 0) THEN
  3666. ! Print*,'Flerchinger used in NEW version. Iterations=',NLOG
  3667. ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  3668. ! CCCCC OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0 CCCCCCCC
  3669. ! CCCCCCCCCCCCC IN KOREN ET AL., JGR, 1999, EQN 17 CCCCCCCCCCCCCCC
  3670. FK=(((HLICE/(GS*(-PSIS)))*((TKELV-T0)/TKELV))**(-1/BX))*SMCMAX
  3671. ! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
  3672. IF (FK .LT. 0.02) FK = 0.02
  3673. FRH2O_init = MIN ( FK, smc )
  3674. ! CCCCCCCCCCCCCCCCCCCCCCCCC END OPTION 2 CCCCCCCCCCCCCCCCCCCCCCCCCC
  3675. ENDIF
  3676. ENDIF
  3677. RETURN
  3678. END FUNCTION FRH2O_init
  3679. !--------------------------------------------------------------------
  3680. SUBROUTINE init_module_initialize
  3681. END SUBROUTINE init_module_initialize
  3682. !---------------------------------------------------------------------
  3683. #ifdef HWRF
  3684. ! compute earth lat-lons for before interpolations. This is gopal's doing for ocean coupling
  3685. !============================================================================================
  3686. SUBROUTINE EARTH_LATLON_hwrf ( HLAT,HLON,VLAT,VLON, & !Earth lat,lon at H and V points
  3687. DLMD1,DPHD1,WBD1,SBD1, & !input res,west & south boundaries,
  3688. CENTRAL_LAT,CENTRAL_LON, & ! central lat,lon, all in degrees
  3689. IDS,IDE,JDS,JDE,KDS,KDE, &
  3690. IMS,IME,JMS,JME,KMS,KME, &
  3691. ITS,ITE,JTS,JTE,KTS,KTE )
  3692. !============================================================================
  3693. !
  3694. IMPLICIT NONE
  3695. INTEGER, INTENT(IN ) :: IDS,IDE,JDS,JDE,KDS,KDE
  3696. INTEGER, INTENT(IN ) :: IMS,IME,JMS,JME,KMS,KME
  3697. INTEGER, INTENT(IN ) :: ITS,ITE,JTS,JTE,KTS,KTE
  3698. REAL, INTENT(IN ) :: DLMD1,DPHD1,WBD1,SBD1
  3699. REAL, INTENT(IN ) :: CENTRAL_LAT,CENTRAL_LON
  3700. REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HLAT,HLON,VLAT,VLON
  3701. ! local
  3702. INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13)
  3703. INTEGER :: I,J
  3704. REAL(KIND=KNUM) :: WB,SB,DLM,DPH,TPH0,STPH0,CTPH0
  3705. REAL(KIND=KNUM) :: TDLM,TDPH,TLMH,TLMV,TLMH0,TLMV0,TPHH,TPHV,DTR
  3706. REAL(KIND=KNUM) :: STPH,CTPH,STPV,CTPV,PI_2
  3707. REAL(KIND=KNUM) :: SPHH,CLMH,FACTH,SPHV,CLMV,FACTV
  3708. REAL(KIND=KNUM), DIMENSION(IMS:IME,JMS:JME) :: GLATH,GLONH,GLATV,GLONV
  3709. !-------------------------------------------------------------------------
  3710. !
  3711. PI_2 = ACOS(0.)
  3712. DTR = PI_2/90.
  3713. WB = WBD1 * DTR ! WB: western boundary in radians
  3714. SB = SBD1 * DTR ! SB: southern boundary in radians
  3715. DLM = DLMD1 * DTR ! DLM: dlamda in radians
  3716. DPH = DPHD1 * DTR ! DPH: dphi in radians
  3717. TDLM = DLM + DLM ! TDLM: 2.0*dlamda
  3718. TDPH = DPH + DPH ! TDPH: 2.0*DPH
  3719. ! For earth lat lon only
  3720. TPH0 = CENTRAL_LAT*DTR ! TPH0: central lat in radians
  3721. STPH0 = SIN(TPH0)
  3722. CTPH0 = COS(TPH0)
  3723. DO J = JTS,MIN(JTE,JDE) !-1) ! H./ This loop takes care of zig-zag
  3724. ! ! \.H starting points along j
  3725. TLMH0 = WB - TDLM + MOD(J+1,2) * DLM ! ./ TLMH (rotated lats at H points)
  3726. TLMV0 = WB - TDLM + MOD(J,2) * DLM ! H (//ly for V points)
  3727. TPHH = SB + (J-1)*DPH ! TPHH (rotated lons at H points) are simple trans.
  3728. TPHV = TPHH ! TPHV (rotated lons at V points) are simple trans.
  3729. STPH = SIN(TPHH)
  3730. CTPH = COS(TPHH)
  3731. STPV = SIN(TPHV)
  3732. CTPV = COS(TPHV)
  3733. ! .H
  3734. DO I = ITS,MIN(ITE,IDE) !-1) ! /
  3735. TLMH = TLMH0 + I*TDLM ! \.H .U .H
  3736. ! !H./ ----><----
  3737. SPHH = CTPH0 * STPH + STPH0 * CTPH * COS(TLMH) ! DLM + DLM
  3738. GLATH(I,J)=ASIN(SPHH) ! GLATH: Earth Lat in radians
  3739. CLMH = CTPH*COS(TLMH)/(COS(GLATH(I,J))*CTPH0) &
  3740. - TAN(GLATH(I,J))*TAN(TPH0)
  3741. IF(CLMH .GT. 1.) CLMH = 1.0
  3742. IF(CLMH .LT. -1.) CLMH = -1.0
  3743. FACTH = 1.
  3744. IF(TLMH .GT. 0.) FACTH = -1.
  3745. GLONH(I,J) = -CENTRAL_LON*DTR + FACTH*ACOS(CLMH)
  3746. ENDDO
  3747. DO I = ITS,MIN(ITE,IDE) !-1)
  3748. TLMV = TLMV0 + I*TDLM
  3749. SPHV = CTPH0 * STPV + STPH0 * CTPV * COS(TLMV)
  3750. GLATV(I,J) = ASIN(SPHV)
  3751. CLMV = CTPV*COS(TLMV)/(COS(GLATV(I,J))*CTPH0) &
  3752. - TAN(GLATV(I,J))*TAN(TPH0)
  3753. IF(CLMV .GT. 1.) CLMV = 1.
  3754. IF(CLMV .LT. -1.) CLMV = -1.
  3755. FACTV = 1.
  3756. IF(TLMV .GT. 0.) FACTV = -1.
  3757. GLONV(I,J) = -CENTRAL_LON*DTR + FACTV*ACOS(CLMV)
  3758. ENDDO
  3759. ENDDO
  3760. ! Conversion to degrees (may not be required, eventually)
  3761. DO J = JTS, MIN(JTE,JDE) !-1)
  3762. DO I = ITS, MIN(ITE,IDE) !-1)
  3763. HLAT(I,J) = GLATH(I,J) / DTR
  3764. HLON(I,J)= -GLONH(I,J)/DTR
  3765. IF(HLON(I,J) .GT. 180.) HLON(I,J) = HLON(I,J) - 360.
  3766. IF(HLON(I,J) .LT. -180.) HLON(I,J) = HLON(I,J) + 360.
  3767. !
  3768. VLAT(I,J) = GLATV(I,J) / DTR
  3769. VLON(I,J) = -GLONV(I,J) / DTR
  3770. IF(VLON(I,J) .GT. 180.) VLON(I,J) = VLON(I,J) - 360.
  3771. IF(VLON(I,J) .LT. -180.) VLON(I,J) = VLON(I,J) + 360.
  3772. ENDDO
  3773. ENDDO
  3774. END SUBROUTINE EARTH_LATLON_hwrf
  3775. SUBROUTINE G2T2H_hwrf( SM,HRES_SM, & ! output grid index and weights
  3776. HLAT,HLON, & ! target (nest) input lat lon in degrees
  3777. DLMD1,DPHD1,WBD1,SBD1, & ! parent res, west and south boundaries
  3778. CENTRAL_LAT,CENTRAL_LON, & ! parent central lat,lon, all in degrees
  3779. P_IDE,P_JDE,P_IMS,P_IME,P_JMS,P_JME, & ! parent imax and jmax
  3780. IDS,IDE,JDS,JDE,KDS,KDE, & ! target (nest) dIMEnsions
  3781. IMS,IME,JMS,JME,KMS,KME, &
  3782. ITS,ITE,JTS,JTE,KTS,KTE )
  3783. !
  3784. !*** Tom Black - Initial Version
  3785. !*** Gopal - Revised Version for WRF (includes coincident grid points)
  3786. !***
  3787. !*** GIVEN PARENT CENTRAL LAT-LONS, RESOLUTION AND WESTERN AND SOUTHERN BOUNDARY,
  3788. !*** AND THE NESTED GRID LAT-LONS AT H POINTS, THIS ROUTINE FIRST LOCATES THE
  3789. !*** INDICES,IIH,JJH, OF THE PARENT DOMAIN'S H POINTS THAT LIES CLOSEST TO THE
  3790. !*** h POINTS OF THE NESTED DOMAIN
  3791. !
  3792. !============================================================================
  3793. !
  3794. IMPLICIT NONE
  3795. INTEGER, INTENT(IN ) :: IDS,IDE,JDS,JDE,KDS,KDE
  3796. INTEGER, INTENT(IN ) :: IMS,IME,JMS,JME,KMS,KME
  3797. INTEGER, INTENT(IN ) :: ITS,ITE,JTS,JTE,KTS,KTE
  3798. INTEGER, INTENT(IN ) :: P_IDE,P_JDE,P_IMS,P_IME,P_JMS,P_JME
  3799. REAL, INTENT(IN ) :: DLMD1,DPHD1,WBD1,SBD1
  3800. REAL, INTENT(IN ) :: CENTRAL_LAT,CENTRAL_LON
  3801. REAL, DIMENSION(P_IMS:P_IME,P_JMS:P_JME), INTENT(IN) :: SM
  3802. REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(IN) :: HLAT,HLON
  3803. REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HRES_SM
  3804. ! local
  3805. INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13)
  3806. INTEGER :: IMT,JMT,N2R,MK,K,I,J,DSLP0,DSLOPE,N
  3807. INTEGER :: NROW,NCOL,KROWS
  3808. REAL(KIND=KNUM) :: X,Y,Z,TLAT,TLON
  3809. REAL(KIND=KNUM) :: PI_2,D2R,R2D,GLAT,GLON,DPH,DLM,TPH0,TLM0,WB,SB
  3810. REAL(KIND=KNUM) :: ROW,COL,SLP0,TLATHC,TLONHC,DENOM,SLOPE
  3811. REAL(KIND=KNUM) :: TLAT1,TLAT2,TLON1,TLON2,DLM1,DLM2,DLM3,DLM4,D1,D2
  3812. REAL(KIND=KNUM) :: DLA1,DLA2,DLA3,DLA4,S1,R1,DS1,AN1,AN2,AN3 ! Q
  3813. REAL(KIND=KNUM) :: DL1,DL2,DL3,DL4,DL1I,DL2I,DL3I,DL4I,SUMDL,TLONO,TLATO
  3814. REAL(KIND=KNUM) :: DTEMP
  3815. REAL , DIMENSION(IMS:IME,JMS:JME) :: TLATHX,TLONHX
  3816. INTEGER, DIMENSION(IMS:IME,JMS:JME) :: KOUTB
  3817. REAL SUM,AMAXVAL
  3818. REAL, DIMENSION (4, ims:ime, jms:jme ) :: NBWGT
  3819. LOGICAL FLIP
  3820. REAL, DIMENSION(IMS:IME,JMS:JME) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
  3821. INTEGER, DIMENSION(IMS:IME,JMS:JME) :: IIH,JJH
  3822. !-------------------------------------------------------------------------------
  3823. IMT=2*P_IDE-2 ! parent i dIMEnsions
  3824. JMT=P_JDE/2 ! parent j dIMEnsions
  3825. PI_2=ACOS(0.)
  3826. D2R=PI_2/90.
  3827. R2D=1./D2R
  3828. DPH=DPHD1*D2R
  3829. DLM=DLMD1*D2R
  3830. TPH0= CENTRAL_LAT*D2R
  3831. TLM0=-CENTRAL_LON*D2R ! NOTE THE MINUS HERE
  3832. WB=WBD1*D2R ! CONVERT NESTED GRID H POINTS FROM GEODETIC
  3833. SB=SBD1*D2R
  3834. SLP0=DPHD1/DLMD1
  3835. DSLP0=NINT(R2D*ATAN(SLP0))
  3836. DS1=SQRT(DPH*DPH+DLM*DLM) ! Q
  3837. AN1=ASIN(DLM/DS1)
  3838. AN2=ASIN(DPH/DS1)
  3839. DO J = JTS,MIN(JTE,JDE) !-1)
  3840. DO I = ITS,MIN(ITE,IDE) !-1)
  3841. !***
  3842. !*** LOCATE TARGET h POINTS (HLAT AND HLON) ON THE PARENT DOMAIN AND
  3843. !*** DETERMINE THE INDICES IN TERMS OF THE PARENT DOMAIN. FIRST
  3844. !*** CONVERT NESTED GRID h POINTS FROM GEODETIC TO TRANSFORMED
  3845. !*** COORDINATE ON THE PARENT GRID
  3846. !
  3847. GLAT=HLAT(I,J)*D2R
  3848. GLON= (360. - HLON(I,J))*D2R
  3849. X=COS(TPH0)*COS(GLAT)*COS(GLON-TLM0)+SIN(TPH0)*SIN(GLAT)
  3850. Y=-COS(GLAT)*SIN(GLON-TLM0)
  3851. Z=COS(TPH0)*SIN(GLAT)-SIN(TPH0)*COS(GLAT)*COS(GLON-TLM0)
  3852. TLAT=R2D*ATAN(Z/SQRT(X*X+Y*Y))
  3853. TLON=R2D*ATAN(Y/X)
  3854. !
  3855. ROW=TLAT/DPHD1+JMT ! JMT IS THE CENTRAL ROW OF THE PARENT DOMAIN
  3856. COL=TLON/DLMD1+P_IDE-1 ! (P_IDE-1) IS THE CENTRAL COLUMN OF THE PARENT DOMAIN
  3857. NROW=INT(ROW + 0.001) ! ROUND-OFF IS AVOIDED WITHOUT USING NINT ON PURPOSE
  3858. NCOL=INT(COL + 0.001)
  3859. TLAT=TLAT*D2R
  3860. TLON=TLON*D2R
  3861. ! WRITE(60,*)'============================================================'
  3862. ! WRITE(60,*)' ','i=',i,'j=',j
  3863. !***
  3864. !***
  3865. !*** FIRST CONSIDER THE SITUATION WHERE THE POINT h IS AT
  3866. !***
  3867. !*** V H
  3868. !***
  3869. !***
  3870. !*** h
  3871. !*** H V
  3872. !***
  3873. !*** THEN LOCATE THE NEAREST H POINT ON THE PARENT GRID
  3874. !***
  3875. IF(MOD(NROW,2).EQ.1.AND.MOD(NCOL,2).EQ.1.OR. &
  3876. MOD(NROW,2).EQ.0.AND.MOD(NCOL,2).EQ.0)THEN
  3877. TLAT1=(NROW-JMT)*DPH
  3878. TLAT2=TLAT1+DPH
  3879. TLON1=(NCOL-(P_IDE-1))*DLM
  3880. TLON2=TLON1+DLM
  3881. DLM1=TLON-TLON1
  3882. DLM2=TLON-TLON2
  3883. ! D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
  3884. ! D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
  3885. DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
  3886. D1=ACOS(DTEMP)
  3887. DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
  3888. D2=ACOS(DTEMP)
  3889. IF(D1.GT.D2)THEN
  3890. NROW=NROW+1 ! FIND THE NEAREST H ROW
  3891. NCOL=NCOL+1 ! FIND THE NEAREST H COLUMN
  3892. ENDIF
  3893. ! WRITE(60,*)'NEAREST PARENT IS:','col=',COL,'row=',ROW,'ncol=',NCOL,'nrow=',NROW
  3894. ELSE
  3895. !***
  3896. !*** NOW CONSIDER THE SITUATION WHERE THE POINT h IS AT
  3897. !***
  3898. !*** H V
  3899. !***
  3900. !***
  3901. !*** h
  3902. !*** V H
  3903. !***
  3904. !*** THEN LOCATE THE NEAREST H POINT ON THE PARENT GRID
  3905. !***
  3906. !***
  3907. TLAT1=(NROW+1-JMT)*DPH
  3908. TLAT2=TLAT1-DPH
  3909. TLON1=(NCOL-(P_IDE-1))*DLM
  3910. TLON2=TLON1+DLM
  3911. DLM1=TLON-TLON1
  3912. DLM2=TLON-TLON2
  3913. ! D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
  3914. ! D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
  3915. DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
  3916. D1=ACOS(DTEMP)
  3917. DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
  3918. D2=ACOS(DTEMP)
  3919. IF(D1.LT.D2)THEN
  3920. NROW=NROW+1 ! FIND THE NEAREST H ROW
  3921. ELSE
  3922. NCOL=NCOL+1 ! FIND THE NEAREST H COLUMN
  3923. ENDIF
  3924. ! WRITE(60,*)'NEAREST PARENT IS:','col=',COL,'row=',ROW,'ncol=',NCOL,'nrow=',NROW
  3925. ENDIF
  3926. KROWS=((NROW-1)/2)*IMT
  3927. IF(MOD(NROW,2).EQ.1)THEN
  3928. K=KROWS+(NCOL+1)/2
  3929. ELSE
  3930. K=KROWS+P_IDE-1+NCOL/2
  3931. ENDIF
  3932. !***
  3933. !*** WE NOW KNOW THAT THE INNER GRID POINT IN QUESTION IS
  3934. !*** NEAREST TO THE CENTER K AS SEEN BELOW. WE MUST FIND
  3935. !*** WHICH OF THE FOUR H-BOXES (OF WHICH THIS H POINT IS
  3936. !*** A VERTEX) SURROUNDS THE INNER GRID h POINT IN QUESTION.
  3937. !***
  3938. !**
  3939. !*** H
  3940. !***
  3941. !***
  3942. !***
  3943. !*** H V H
  3944. !***
  3945. !***
  3946. !*** h
  3947. !*** H V H V H
  3948. !***
  3949. !***
  3950. !***
  3951. !*** H V H
  3952. !***
  3953. !***
  3954. !***
  3955. !*** H
  3956. !***
  3957. !***
  3958. !*** FIND THE SLOPE OF THE LINE CONNECTING h AND THE CENTER H.
  3959. !***
  3960. N2R=K/IMT
  3961. MK=MOD(K,IMT)
  3962. !
  3963. IF(MK.EQ.0)THEN
  3964. TLATHC=SB+(2*N2R-1)*DPH
  3965. ELSE
  3966. TLATHC=SB+(2*N2R+(MK-1)/(P_IDE-1))*DPH
  3967. ENDIF
  3968. !
  3969. IF(MK.LE.(P_IDE-1))THEN
  3970. TLONHC=WB+2*(MK-1)*DLM
  3971. ELSE
  3972. TLONHC=WB+(2*(MK-(P_IDE-1))-1)*DLM
  3973. ENDIF
  3974. !
  3975. !*** EXECUTE CAUTION IF YOU NEED TO CHANGE THESE CONDITIONS. SINCE WE ARE
  3976. !*** DEALING WITH SLOPES TO GENERATE DIAMOND SHAPE H BOXES, WE NEED TO BE
  3977. !*** CAREFUL HERE
  3978. !
  3979. IF(ABS(TLON-TLONHC) .LE. 1.E-4)TLONHC=TLON
  3980. IF(ABS(TLAT-TLATHC) .LE. 1.E-4)TLATHC=TLAT
  3981. DENOM=(TLON-TLONHC)
  3982. !
  3983. !***
  3984. !***STORE THE LOCATION OF THE WESTERNMOST VERTEX OF THE H-BOX ON
  3985. !***THE OUTER GRID THAT SURROUNDS THE h POINT ON THE INNER GRID.
  3986. !***
  3987. !*** COINCIDENT CONDITIONS
  3988. IF(DENOM.EQ.0.0)THEN
  3989. IF(TLATHC.EQ.TLAT)THEN
  3990. KOUTB(I,J)=K
  3991. IIH(I,J) = NCOL
  3992. JJH(I,J) = NROW
  3993. TLATHX(I,J)=TLATHC
  3994. TLONHX(I,J)=TLONHC
  3995. HBWGT1(I,J)=1.0
  3996. HBWGT2(I,J)=0.0
  3997. HBWGT3(I,J)=0.0
  3998. HBWGT4(I,J)=0.0
  3999. ! WRITE(60,*)'TRIVIAL SOLUTION'
  4000. ELSE ! SAME LONGITUDE BUT DIFFERENT LATS
  4001. !
  4002. IF(TLATHC .GT. TLAT)THEN ! NESTED POINT SOUTH OF PARENT
  4003. KOUTB(I,J)=K-(P_IDE-1)
  4004. IIH(I,J) = NCOL-1
  4005. JJH(I,J) = NROW-1
  4006. TLATHX(I,J)=TLATHC-DPH
  4007. TLONHX(I,J)=TLONHC-DLM
  4008. ! WRITE(60,*)'VANISHING SLOPE, -ve: TLATHC-DPH, TLONHC-DLM'
  4009. ELSE ! NESTED POINT NORTH OF PARENT
  4010. KOUTB(I,J)=K+(P_IDE-1)-1
  4011. IIH(I,J) = NCOL-1
  4012. JJH(I,J) = NROW+1
  4013. TLATHX(I,J)=TLATHC+DPH
  4014. TLONHX(I,J)=TLONHC-DLM
  4015. ! WRITE(60,*)'VANISHING SLOPE, +ve: TLATHC+DPH, TLONHC-DLM'
  4016. ENDIF
  4017. !***
  4018. !***
  4019. !*** 4
  4020. !***
  4021. !*** h
  4022. !*** 1 2
  4023. !***
  4024. !*** 3
  4025. !*** DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX
  4026. TLATO=TLATHX(I,J)
  4027. TLONO=TLONHX(I,J)
  4028. DLM1=TLON-TLONO
  4029. DLA1=TLAT-TLATO ! Q
  4030. ! DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
  4031. DL1=SQRT(DLM1*DLM1+DLA1*DLA1) ! Q
  4032. !
  4033. TLATO=TLATHX(I,J)
  4034. TLONO=TLONHX(I,J)+2.*DLM
  4035. DLM2=TLON-TLONO
  4036. DLA2=TLAT-TLATO ! Q
  4037. ! DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
  4038. DL2=SQRT(DLM2*DLM2+DLA2*DLA2) ! Q
  4039. !
  4040. TLATO=TLATHX(I,J)-DPH
  4041. TLONO=TLONHX(I,J)+DLM
  4042. DLM3=TLON-TLONO
  4043. DLA3=TLAT-TLATO ! Q
  4044. ! DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
  4045. DL3=SQRT(DLM3*DLM3+DLA3*DLA3) ! Q
  4046. TLATO=TLATHX(I,J)+DPH
  4047. TLONO=TLONHX(I,J)+DLM
  4048. DLM4=TLON-TLONO
  4049. DLA4=TLAT-TLATO ! Q
  4050. ! DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
  4051. DL4=SQRT(DLM4*DLM4+DLA4*DLA4) ! Q
  4052. ! THE BILINEAR WEIGHTS
  4053. !***
  4054. !***
  4055. AN3=ATAN2(DLA1,DLM1) ! Q
  4056. R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
  4057. S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
  4058. R1=R1/DS1
  4059. S1=S1/DS1
  4060. DL1I=(1.-R1)*(1.-S1)
  4061. DL2I=R1*S1
  4062. DL3I=R1*(1.-S1)
  4063. DL4I=(1.-R1)*S1
  4064. !
  4065. HBWGT1(I,J)=DL1I
  4066. HBWGT2(I,J)=DL2I
  4067. HBWGT3(I,J)=DL3I
  4068. HBWGT4(I,J)=DL4I
  4069. !
  4070. ENDIF
  4071. ELSE
  4072. !
  4073. !*** NON-COINCIDENT POINTS
  4074. !
  4075. SLOPE=(TLAT-TLATHC)/DENOM
  4076. DSLOPE=NINT(R2D*ATAN(SLOPE))
  4077. IF(DSLOPE.LE.DSLP0.AND.DSLOPE.GE.-DSLP0)THEN
  4078. IF(TLON.GT.TLONHC)THEN
  4079. ! IF(TLONHC.GE.-WB-DLM)CALL wrf_error_fatal("1H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
  4080. KOUTB(I,J)=K
  4081. IIH(I,J) = NCOL
  4082. JJH(I,J) = NROW
  4083. TLATHX(I,J)=TLATHC
  4084. TLONHX(I,J)=TLONHC
  4085. ! WRITE(60,*)'HERE WE GO1: TLATHC, TLONHC'
  4086. ELSE
  4087. ! IF(TLONHC.LE.WB+DLM)CALL wrf_error_fatal("2H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
  4088. KOUTB(I,J)=K-1
  4089. IIH(I,J) = NCOL-2
  4090. JJH(I,J) = NROW
  4091. TLATHX(I,J)=TLATHC
  4092. TLONHX(I,J)=TLONHC -2.*DLM
  4093. ! WRITE(60,*)'HERE WE GO2: TLATHC, TLONHC -2.*DLM'
  4094. ENDIF
  4095. !
  4096. ELSEIF(DSLOPE.GT.DSLP0)THEN
  4097. IF(TLON.GT.TLONHC)THEN
  4098. ! IF(TLATHC.GE.-SB-DPH)CALL wrf_error_fatal("3H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
  4099. KOUTB(I,J)=K+(P_IDE-1)-1
  4100. IIH(I,J) = NCOL-1
  4101. JJH(I,J) = NROW+1
  4102. TLATHX(I,J)=TLATHC+DPH
  4103. TLONHX(I,J)=TLONHC-DLM
  4104. ! WRITE(60,*)'HERE WE GO3: TLATHC+DPH, TLONHC-DLM'
  4105. ELSE
  4106. ! IF(TLATHC.LE.SB+DPH)CALL wrf_error_fatal("4H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
  4107. KOUTB(I,J)=K-(P_IDE-1)
  4108. IIH(I,J) = NCOL-1
  4109. JJH(I,J) = NROW-1
  4110. TLATHX(I,J)=TLATHC-DPH
  4111. TLONHX(I,J)=TLONHC-DLM
  4112. ! WRITE(60,*)'HERE WE GO4: TLATHC-DPH, TLONHC-DLM'
  4113. ENDIF
  4114. !
  4115. ELSEIF(DSLOPE.LT.-DSLP0)THEN
  4116. IF(TLON.GT.TLONHC)THEN
  4117. ! IF(TLATHC.LE.SB+DPH)CALL wrf_error_fatal("5H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
  4118. KOUTB(I,J)=K-(P_IDE-1)
  4119. IIH(I,J) = NCOL-1
  4120. JJH(I,J) = NROW-1
  4121. TLATHX(I,J)=TLATHC-DPH
  4122. TLONHX(I,J)=TLONHC-DLM
  4123. ! WRITE(60,*)'HERE WE GO5: TLATHC-DPH, TLONHC-DLM'
  4124. ELSE
  4125. ! IF(TLATHC.GE.-SB-DPH)CALL wrf_error_fatal("6H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
  4126. KOUTB(I,J)=K+(P_IDE-1)-1
  4127. IIH(I,J) = NCOL-1
  4128. JJH(I,J) = NROW+1
  4129. TLATHX(I,J)=TLATHC+DPH
  4130. TLONHX(I,J)=TLONHC-DLM
  4131. ! WRITE(60,*)'HERE WE GO6: TLATHC+DPH, TLONHC-DLM'
  4132. ENDIF
  4133. ENDIF
  4134. !
  4135. !*** NOW WE WILL MOVE AS FOLLOWS:
  4136. !***
  4137. !***
  4138. !*** 4
  4139. !***
  4140. !***
  4141. !***
  4142. !*** h
  4143. !*** 1 2
  4144. !***
  4145. !***
  4146. !***
  4147. !***
  4148. !*** 3
  4149. !***
  4150. !***
  4151. !***
  4152. !*** DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX
  4153. TLATO=TLATHX(I,J)
  4154. TLONO=TLONHX(I,J)
  4155. DLM1=TLON-TLONO
  4156. DLA1=TLAT-TLATO ! Q
  4157. ! DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
  4158. DL1=SQRT(DLM1*DLM1+DLA1*DLA1) ! Q
  4159. !
  4160. TLATO=TLATHX(I,J) ! redundant computations
  4161. TLONO=TLONHX(I,J)+2.*DLM
  4162. DLM2=TLON-TLONO
  4163. DLA2=TLAT-TLATO ! Q
  4164. ! DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
  4165. DL2=SQRT(DLM2*DLM2+DLA2*DLA2) ! Q
  4166. !
  4167. TLATO=TLATHX(I,J)-DPH
  4168. TLONO=TLONHX(I,J)+DLM
  4169. DLM3=TLON-TLONO
  4170. DLA3=TLAT-TLATO ! Q
  4171. ! DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
  4172. DL3=SQRT(DLM3*DLM3+DLA3*DLA3) ! Q
  4173. !
  4174. TLATO=TLATHX(I,J)+DPH
  4175. TLONO=TLONHX(I,J)+DLM
  4176. DLM4=TLON-TLONO
  4177. DLA4=TLAT-TLATO ! Q
  4178. ! DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
  4179. DL4=SQRT(DLM4*DLM4+DLA4*DLA4) ! Q
  4180. ! THE BILINEAR WEIGHTS
  4181. !***
  4182. AN3=ATAN2(DLA1,DLM1) ! Q
  4183. R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
  4184. S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
  4185. R1=R1/DS1
  4186. S1=S1/DS1
  4187. DL1I=(1.-R1)*(1.-S1)
  4188. DL2I=R1*S1
  4189. DL3I=R1*(1.-S1)
  4190. DL4I=(1.-R1)*S1
  4191. !
  4192. HBWGT1(I,J)=DL1I
  4193. HBWGT2(I,J)=DL2I
  4194. HBWGT3(I,J)=DL3I
  4195. HBWGT4(I,J)=DL4I
  4196. !
  4197. ENDIF
  4198. !
  4199. !*** FINALLY STORE IIH IN TERMS OF E-GRID INDEX
  4200. !
  4201. IIH(I,J)=NINT(0.5*IIH(I,J))
  4202. ENDDO
  4203. ENDDO
  4204. !
  4205. !*** EXTENSION TO NEAREST NEIGHBOR
  4206. !
  4207. DO J = JTS,MIN(JTE,JDE) !-1)
  4208. DO I = ITS,MIN(ITE,IDE) !-1)
  4209. NBWGT(1,I,J)=HBWGT1(I,J)
  4210. NBWGT(2,I,J)=HBWGT2(I,J)
  4211. NBWGT(3,I,J)=HBWGT3(I,J)
  4212. NBWGT(4,I,J)=HBWGT4(I,J)
  4213. ENDDO
  4214. ENDDO
  4215. DO J = JTS,MIN(JTE,JDE) !-1)
  4216. DO I = ITS,MIN(ITE,IDE) !-1)
  4217. AMAXVAL=0.
  4218. DO N=1,4
  4219. AMAXVAL=amax1(NBWGT(N,I,J),AMAXVAL)
  4220. ENDDO
  4221. !
  4222. FLIP=.TRUE.
  4223. SUM=0.0
  4224. DO N=1,4
  4225. IF(AMAXVAL .EQ. NBWGT(N,I,J) .AND. FLIP)THEN
  4226. NBWGT(N,I,J)=1.0
  4227. FLIP=.FALSE.
  4228. ELSE
  4229. NBWGT(N,I,J)=0.0
  4230. ENDIF
  4231. SUM=SUM+NBWGT(N,I,J)
  4232. IF(SUM .GT. 1.0)CALL wrf_error_fatal ( "horizontal interp error - interp_hnear_nmm" )
  4233. ENDDO
  4234. IF((NBWGT(1,I,J)+NBWGT(2,I,J)+NBWGT(3,I,J)+NBWGT(4,I,J)) .NE. 1)THEN
  4235. WRITE(0,*)'------------------------------------------------------------------------'
  4236. WRITE(0,*)'FATAL: SOMETHING IS WRONG WITH THE WEIGHTS IN module_initialize_real.F'
  4237. WRITE(0,*)'------------------------------------------------------------------------'
  4238. STOP
  4239. ENDIF
  4240. ! WRITE(66,*)I,J,NBWGT(1,I,J),NBWGT(2,I,J),NBWGT(3,I,J),NBWGT(4,I,J)
  4241. ENDDO
  4242. ENDDO
  4243. DO J=MAX(3,JTS),MIN(JTE,JDE) !-1)
  4244. DO I=MAX(3,ITS),MIN(ITE,IDE) !-1)
  4245. IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
  4246. HRES_SM(I,J) = NBWGT(1,I,J)*SM(IIH(I,J),JJH(I,J) ) &
  4247. + NBWGT(2,I,J)*SM(IIH(I,J)+1, JJH(I,J) ) &
  4248. + NBWGT(3,I,J)*SM(IIH(I,J), JJH(I,J)-1) &
  4249. + NBWGT(4,I,J)*SM(IIH(I,J), JJH(I,J)+1)
  4250. ! WRITE(68,*)I,J,SM(IIH(I,J),JJH(I,J)),SM(IIH(I,J)+1, JJH(I,J)), &
  4251. ! SM(IIH(I,J), JJH(I,J)-1),SM(IIH(I,J), JJH(I,J)+1),HRES_SM(I,J)
  4252. ELSE
  4253. HRES_SM(I,J) = NBWGT(1,I,J)*SM(IIH(I,J), JJH(I,J) ) &
  4254. + NBWGT(2,I,J)*SM(IIH(I,J)+1, JJH(I,J) ) &
  4255. + NBWGT(3,I,J)*SM(IIH(I,J)+1, JJH(I,J)-1) &
  4256. + NBWGT(4,I,J)*SM(IIH(I,J)+1, JJH(I,J)+1)
  4257. ! WRITE(68,*)I,J,SM(IIH(I,J),JJH(I,J)),SM(IIH(I,J)+1, JJH(I,J)), &
  4258. ! SM(IIH(I,J)+1, JJH(I,J)-1),SM(IIH(I,J)+1, JJH(I,J)+1),HRES_SM(I,J)
  4259. ENDIF
  4260. ENDDO
  4261. ENDDO
  4262. ! Boundary treatment in J direction
  4263. DO J=MAX(3,JTS),MIN(JTE,JDE)
  4264. HRES_SM(2,J)=HRES_SM(3,J)
  4265. HRES_SM(1,J)=HRES_SM(2,J)
  4266. END DO
  4267. ! Boundary treatment in J direction and 4 corners
  4268. DO I=ITS,MIN(ITE,IDE)
  4269. HRES_SM(I,2)=HRES_SM(I,3)
  4270. HRES_SM(I,1)=HRES_SM(I,2)
  4271. END DO
  4272. RETURN
  4273. END SUBROUTINE G2T2H_hwrf
  4274. !========================================================================================
  4275. ! end gopal's doing for ocean coupling
  4276. !============================================================================================
  4277. #endif
  4278. END MODULE module_initialize_real