| blob | 03e512cd6a2beb5ead22b0e668105b5791796b39 |
1 !REAL:MODEL_LAYER:INITIALIZATION
3 ! This MODULE holds the routines which are used to perform various initializations
4 ! for individual domains utilizing the NMM dynamical core.
6 !-----------------------------------------------------------------------
8 MODULE module_initialize_real
10 USE module_bc
11 USE module_configure
12 USE module_domain
13 USE module_io_domain
14 USE module_model_constants
15 ! USE module_si_io_nmm
16 USE module_state_description
17 USE module_timing
18 USE module_soil_pre
19 #ifdef DM_PARALLEL
20 USE module_dm, ONLY : LOCAL_COMMUNICATOR &
21 ,MYTASK,NTASKS,NTASKS_X &
22 ,NTASKS_Y
23 USE module_comm_dm
24 USE module_ext_internal
25 #endif
27 INTEGER :: internal_time_loop
28 INTEGER:: MPI_COMM_COMP,MYPE
29 INTEGER:: loopinc, iloopinc
31 CONTAINS
33 !-------------------------------------------------------------------
35 SUBROUTINE init_domain ( grid )
37 IMPLICIT NONE
39 ! Input space and data. No gridded meteorological data has been stored, though.
41 ! TYPE (domain), POINTER :: grid
42 TYPE (domain) :: grid
44 ! Local data.
46 INTEGER :: idum1, idum2
48 CALL set_scalar_indices_from_config ( head_grid%id , idum1, idum2 )
50 CALL init_domain_nmm (grid &
51 !
52 #include <actual_new_args.inc>
53 !
54 )
56 END SUBROUTINE init_domain
58 !-------------------------------------------------------------------
59 !---------------------------------------------------------------------
60 SUBROUTINE init_domain_nmm ( grid &
61 !
62 # include <dummy_new_args.inc>
63 !
64 )
66 USE module_optional_input
67 IMPLICIT NONE
69 ! Input space and data. No gridded meteorological data has been stored, though.
71 ! TYPE (domain), POINTER :: grid
72 TYPE (domain) :: grid
74 # include <dummy_new_decl.inc>
76 TYPE (grid_config_rec_type) :: config_flags
78 ! Local domain indices and counters.
80 INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
82 INTEGER :: &
83 ids, ide, jds, jde, kds, kde, &
84 ims, ime, jms, jme, kms, kme, &
85 its, ite, jts, jte, kts, kte, &
86 ips, ipe, jps, jpe, kps, kpe, &
87 i, j, k, NNXP, NNYP
89 ! Local data
91 CHARACTER(LEN=19):: start_date
93 #ifdef DM_PARALLEL
95 LOGICAL,EXTERNAL :: WRF_DM_ON_MONITOR
96 logical :: test
98 ! INTEGER :: DOMDESC
99 REAL,ALLOCATABLE :: SICE_G(:,:), SM_G(:,:)
100 INTEGER, ALLOCATABLE:: IHE_G(:),IHW_G(:)
101 INTEGER, ALLOCATABLE:: ITEMP(:,:)
102 INTEGER :: my_e,my_n,my_s,my_w,my_ne,my_nw,my_se,my_sw,myi,myj,npe
103 INTEGER :: istat,inpes,jnpes
104 #else
105 integer, allocatable:: ihw(:),ihe(:)
106 #endif
109 CHARACTER (LEN=255) :: message
111 INTEGER :: error
112 REAL :: p_surf, p_level
113 REAL :: cof1, cof2
114 REAL :: qvf , qvf1 , qvf2 , pd_surf
115 REAL :: p00 , t00 , a
116 REAL :: hold_znw, rmin,rmax
118 REAL :: p_top_requested , ptsgm
119 INTEGER :: num_metgrid_levels, ICOUNT
120 REAL , DIMENSION(max_eta) :: eta_levels
122 LOGICAL :: stretch_grid, dry_sounding, debug, log_flag_sst, hyb_coor
124 REAL, ALLOCATABLE,DIMENSION(:,:):: ADUM2D,SNOWC,HT,TG_ALT, &
125 PDVP,PSFC_OUTV
127 REAL, ALLOCATABLE,DIMENSION(:,:,:):: P3D_OUT,P3DV_OUT,P3DV_IN, &
128 QTMP,QTMP2
130 INTEGER, ALLOCATABLE, DIMENSION(:):: KHL2,KVL2,KHH2,KVH2, &
131 KHLA,KHHA,KVLA,KVHA
133 ! INTEGER, ALLOCATABLE, DIMENSION(:,:):: grid%lu_index
135 REAL, ALLOCATABLE, DIMENSION(:):: DXJ,WPDARJ,CPGFUJ,CURVJ, &
136 FCPJ,FDIVJ,EMJ,EMTJ,FADJ, &
137 HDACJ,DDMPUJ,DDMPVJ
139 REAL, ALLOCATABLE,DIMENSION(:),SAVE:: SG1,SG2,DSG1,DSG2, &
140 SGML1,SGML2
142 !-- Carsel and Parrish [1988]
143 REAL , DIMENSION(100) :: lqmi
144 integer iicount
146 REAL:: TPH0D,TLM0D
147 REAL:: TPH0,WB,SB,TDLM,TDPH
148 REAL:: WBI,SBI,EBI,ANBI,STPH0,CTPH0
149 REAL:: TSPH,DTAD,DTCF
150 REAL:: ACDT,CDDAMP,DXP,FP
151 REAL:: WBD,SBD
152 REAL:: RSNOW,SNOFAC
153 REAL, PARAMETER:: SALP=2.60
154 REAL, PARAMETER:: SNUP=0.040
155 REAL:: SMCSUM,STCSUM,SEAICESUM,FISX
156 REAL:: cur_smc, aposs_smc
158 REAL:: COAC, CODAMP
160 INTEGER,PARAMETER:: DOUBLE=SELECTED_REAL_KIND(15,300)
161 REAL(KIND=DOUBLE):: TERM1,APH,TLM,TPH,DLM,DPH,STPH,CTPH
163 INTEGER:: KHH,KVH,JAM,JA, IHL, IHH, L
164 INTEGER:: II,JJ,ISRCH,ISUM,ITER,Ilook,Jlook
166 REAL, PARAMETER:: DTR=0.01745329
167 REAL, PARAMETER:: W_NMM=0.08
168 #if defined(HWRF)
169 REAL, PARAMETER:: DDFC=1.0
170 #else
171 REAL, PARAMETER:: DDFC=8.0
172 #endif
173 REAL, PARAMETER:: TWOM=.00014584
174 REAL, PARAMETER:: CP=1004.6
175 REAL, PARAMETER:: DFC=1.0
176 REAL, PARAMETER:: ROI=916.6
177 REAL, PARAMETER:: R=287.04
178 REAL, PARAMETER:: CI=2060.0
179 REAL, PARAMETER:: ROS=1500.
180 REAL, PARAMETER:: CS=1339.2
181 REAL, PARAMETER:: DS=0.050
182 REAL, PARAMETER:: AKS=.0000005
183 REAL, PARAMETER:: DZG=2.85
184 REAL, PARAMETER:: DI=.1000
185 REAL, PARAMETER:: AKI=0.000001075
186 REAL, PARAMETER:: DZI=2.0
187 REAL, PARAMETER:: THL=210.
188 REAL, PARAMETER:: PLQ=70000.
189 REAL, PARAMETER:: ERAD=6371200.
190 REAL, PARAMETER:: TG0=258.16
191 REAL, PARAMETER:: TGA=30.0
192 integer :: numzero,numexamined
193 #ifdef HWRF
194 !============================================================================
195 ! gopal's doing for ocean coupling
196 !============================================================================
198 REAL, DIMENSION(:,:), ALLOCATABLE :: NHLAT,NHLON,NVLAT,NVLON,HRES_SM
199 REAL :: NDLMD,NDPHD,NWBD,NSBD
200 INTEGER :: NIDE,NJDE,ILOC,JLOC
202 INTEGER fid, ierr, nprocs
203 CHARACTER*255 f65name, SysString
205 !============================================================================
206 ! end gopal's doing for ocean coupling
207 !============================================================================
208 #endif
210 if (ALLOCATED(ADUM2D)) DEALLOCATE(ADUM2D)
211 if (ALLOCATED(TG_ALT)) DEALLOCATE(TG_ALT)
213 !#define COPY_IN
214 !#include <scalar_derefs.inc>
215 #ifdef DM_PARALLEL
216 # include <data_calls.inc>
217 #endif
219 SELECT CASE ( model_data_order )
220 CASE ( DATA_ORDER_ZXY )
221 kds = grid%sd31 ; kde = grid%ed31 ;
222 ids = grid%sd32 ; ide = grid%ed32 ;
223 jds = grid%sd33 ; jde = grid%ed33 ;
225 kms = grid%sm31 ; kme = grid%em31 ;
226 ims = grid%sm32 ; ime = grid%em32 ;
227 jms = grid%sm33 ; jme = grid%em33 ;
229 kts = grid%sp31 ; kte = grid%ep31 ; ! tile is entire patch
230 its = grid%sp32 ; ite = grid%ep32 ; ! tile is entire patch
231 jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
233 CASE ( DATA_ORDER_XYZ )
234 ids = grid%sd31 ; ide = grid%ed31 ;
235 jds = grid%sd32 ; jde = grid%ed32 ;
236 kds = grid%sd33 ; kde = grid%ed33 ;
238 ims = grid%sm31 ; ime = grid%em31 ;
239 jms = grid%sm32 ; jme = grid%em32 ;
240 kms = grid%sm33 ; kme = grid%em33 ;
242 its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
243 jts = grid%sp32 ; jte = grid%ep32 ; ! tile is entire patch
244 kts = grid%sp33 ; kte = grid%ep33 ; ! tile is entire patch
246 CASE ( DATA_ORDER_XZY )
247 ids = grid%sd31 ; ide = grid%ed31 ;
248 kds = grid%sd32 ; kde = grid%ed32 ;
249 jds = grid%sd33 ; jde = grid%ed33 ;
251 ims = grid%sm31 ; ime = grid%em31 ;
252 kms = grid%sm32 ; kme = grid%em32 ;
253 jms = grid%sm33 ; jme = grid%em33 ;
255 its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
256 kts = grid%sp32 ; kte = grid%ep32 ; ! tile is entire patch
257 jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
259 END SELECT
261 #ifdef DM_PARALLEL
262 CALL WRF_GET_MYPROC(MYPE)
263 CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP)
264 call wrf_get_nprocx(inpes)
265 call wrf_get_nprocy(jnpes)
266 !
267 allocate(itemp(inpes,jnpes),stat=istat)
268 npe=0
269 !
270 do j=1,jnpes
271 do i=1,inpes
272 itemp(i,j)=npe
273 if(npe==mype)then
274 myi=i
275 myj=j
276 endif
277 npe=npe+1
278 enddo
279 enddo
280 !
281 my_n=-1
282 if(myj+1<=jnpes)my_n=itemp(myi,myj+1)
283 !
284 my_e=-1
285 if(myi+1<=inpes)my_e=itemp(myi+1,myj)
286 !
287 my_s=-1
288 if(myj-1>=1)my_s=itemp(myi,myj-1)
289 !
290 my_w=-1
291 if(myi-1>=1)my_w=itemp(myi-1,myj)
292 !
293 my_ne=-1
294 if((myi+1<=inpes).and.(myj+1<=jnpes)) &
295 my_ne=itemp(myi+1,myj+1)
296 !
297 my_se=-1
298 if((myi+1<=inpes).and.(myj-1>=1)) &
299 my_se=itemp(myi+1,myj-1)
300 !
301 my_sw=-1
302 if((myi-1>=1).and.(myj-1>=1)) &
303 my_sw=itemp(myi-1,myj-1)
304 !
305 my_nw=-1
306 if((myi-1>=1).and.(myj+1<=jnpes)) &
307 my_nw=itemp(myi-1,myj+1)
308 !
309 ! my_neb(1)=my_n
310 ! my_neb(2)=my_e
311 ! my_neb(3)=my_s
312 ! my_neb(4)=my_w
313 ! my_neb(5)=my_ne
314 ! my_neb(6)=my_se
315 ! my_neb(7)=my_sw
316 ! my_neb(8)=my_nw
317 !
318 deallocate(itemp)
319 #endif
321 NNXP=min(ITE,IDE-1)
322 NNYP=min(JTE,JDE-1)
324 write(message,*) 'IDE, JDE: ', IDE,JDE
325 write(message,*) 'NNXP, NNYP: ', NNXP,NNYP
326 CALL wrf_message(message)
328 JAM=6+2*(JDE-JDS-10)
330 if (internal_time_loop .eq. 1) then
331 ALLOCATE(ADUM2D(grid%sm31:grid%em31,jms:jme))
332 ALLOCATE(KHL2(JTS:NNYP),KVL2(JTS:NNYP),KHH2(JTS:NNYP),KVH2(JTS:NNYP))
333 ALLOCATE(DXJ(JTS:NNYP),WPDARJ(JTS:NNYP),CPGFUJ(JTS:NNYP),CURVJ(JTS:NNYP))
334 ALLOCATE(FCPJ(JTS:NNYP),FDIVJ(JTS:NNYP),&
335 FADJ(JTS:NNYP))
336 ALLOCATE(HDACJ(JTS:NNYP),DDMPUJ(JTS:NNYP),DDMPVJ(JTS:NNYP))
337 ALLOCATE(KHLA(JAM),KHHA(JAM))
338 ALLOCATE(KVLA(JAM),KVHA(JAM))
339 endif
342 CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
344 IF ( CONFIG_FLAGS%FRACTIONAL_SEAICE == 1 ) THEN
345 CALL WRF_ERROR_FATAL("NMM cannot use FRACTIONAL_SEAICE = 1")
346 ENDIF
348 if ( config_flags%bl_pbl_physics == BOULACSCHEME ) then
349 call wrf_error_fatal("Cannot use BOULAC PBL with NMM")
350 endif
352 write(message,*) 'cen_lat: ', config_flags%cen_lat
353 CALL wrf_debug(100,message)
354 write(message,*) 'cen_lon: ', config_flags%cen_lon
355 CALL wrf_debug(100,message)
356 write(message,*) 'dx: ', config_flags%dx
357 CALL wrf_debug(100,message)
358 write(message,*) 'dy: ', config_flags%dy
359 CALL wrf_debug(100,message)
360 write(message,*) 'config_flags%start_year: ', config_flags%start_year
361 CALL wrf_debug(100,message)
362 write(message,*) 'config_flags%start_month: ', config_flags%start_month
363 CALL wrf_debug(100,message)
364 write(message,*) 'config_flags%start_day: ', config_flags%start_day
365 CALL wrf_debug(100,message)
366 write(message,*) 'config_flags%start_hour: ', config_flags%start_hour
367 CALL wrf_debug(100,message)
369 write(start_date,435) config_flags%start_year, config_flags%start_month, &
370 config_flags%start_day, config_flags%start_hour
371 435 format(I4,'-',I2.2,'-',I2.2,'_',I2.2,':00:00')
373 grid%dlmd=config_flags%dx
374 grid%dphd=config_flags%dy
375 tph0d=config_flags%cen_lat
376 tlm0d=config_flags%cen_lon
378 !==========================================================================
380 !!
382 ! Check to see if the boundary conditions are set
383 ! properly in the namelist file.
384 ! This checks for sufficiency and redundancy.
386 CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
388 ! Some sort of "this is the first time" initialization. Who knows.
390 grid%itimestep=0
392 ! Pull in the info in the namelist to compare it to the input data.
394 grid%real_data_init_type = model_config_rec%real_data_init_type
395 write(message,*) 'what is flag_metgrid: ', flag_metgrid
396 CALL wrf_message(message)
398 IF ( flag_metgrid .EQ. 1 ) THEN ! <----- START OF VERTICAL INTERPOLATION PART ---->
400 num_metgrid_levels = grid%num_metgrid_levels
403 IF (grid%ght_gc(its,jts,num_metgrid_levels/2) .lt. grid%ght_gc(its,jts,num_metgrid_levels/2+1)) THEN
405 write(message,*) 'normal ground up file order'
406 hyb_coor=.false.
407 CALL wrf_message(message)
409 ELSE
411 hyb_coor=.true.
412 write(message,*) 'reverse the order of coordinate'
413 CALL wrf_message(message)
415 CALL reverse_vert_coord(grid%ght_gc, 2, num_metgrid_levels &
416 &, IDS,IDE,JDS,JDE,KDS,KDE &
417 &, IMS,IME,JMS,JME,KMS,KME &
418 &, ITS,ITE,JTS,JTE,KTS,KTE )
420 #if defined(HWRF)
421 if(.not. grid%use_prep_hybrid) then
422 #endif
424 CALL reverse_vert_coord(grid%p_gc, 2, num_metgrid_levels &
425 &, IDS,IDE,JDS,JDE,KDS,KDE &
426 &, IMS,IME,JMS,JME,KMS,KME &
427 &, ITS,ITE,JTS,JTE,KTS,KTE )
429 CALL reverse_vert_coord(grid%t_gc, 2, num_metgrid_levels &
430 &, IDS,IDE,JDS,JDE,KDS,KDE &
431 &, IMS,IME,JMS,JME,KMS,KME &
432 &, ITS,ITE,JTS,JTE,KTS,KTE )
434 CALL reverse_vert_coord(grid%u_gc, 2, num_metgrid_levels &
435 &, IDS,IDE,JDS,JDE,KDS,KDE &
436 &, IMS,IME,JMS,JME,KMS,KME &
437 &, ITS,ITE,JTS,JTE,KTS,KTE )
439 CALL reverse_vert_coord(grid%v_gc, 2, num_metgrid_levels &
440 &, IDS,IDE,JDS,JDE,KDS,KDE &
441 &, IMS,IME,JMS,JME,KMS,KME &
442 &, ITS,ITE,JTS,JTE,KTS,KTE )
444 CALL reverse_vert_coord(grid%rh_gc, 2, num_metgrid_levels &
445 &, IDS,IDE,JDS,JDE,KDS,KDE &
446 &, IMS,IME,JMS,JME,KMS,KME &
447 &, ITS,ITE,JTS,JTE,KTS,KTE )
449 #if defined(HWRF)
450 endif
451 #endif
453 endif
456 IF (hyb_coor) THEN
457 ! limit extreme deviations from source model topography
458 ! due to potential for nasty extrapolation/interpolation issues
459 !
460 write(message,*) 'min, max of grid%ht_gc before adjust: ', minval(grid%ht_gc), maxval(grid%ht_gc)
461 CALL wrf_debug(100,message)
462 ICOUNT=0
463 DO J=JTS,min(JTE,JDE-1)
464 DO I=ITS,min(ITE,IDE-1)
465 IF ((grid%ht_gc(I,J) - grid%ght_gc(I,J,2)) .LT. -150.) THEN
466 grid%ht_gc(I,J)=grid%ght_gc(I,J,2)-150.
467 IF (ICOUNT .LT. 20) THEN
468 write(message,*) 'increasing NMM topo toward RUC ', I,J
469 CALL wrf_debug(100,message)
470 ICOUNT=ICOUNT+1
471 ENDIF
472 ELSEIF ((grid%ht_gc(I,J) - grid%ght_gc(I,J,2)) .GT. 150.) THEN
473 grid%ht_gc(I,J)=grid%ght_gc(I,J,2)+150.
474 IF (ICOUNT .LT. 20) THEN
475 write(message,*) 'decreasing NMM topo toward RUC ', I,J
476 CALL wrf_debug(100,message)
477 ICOUNT=ICOUNT+1
478 ENDIF
479 ENDIF
480 END DO
481 END DO
483 write(message,*) 'min, max of ht_gc after correct: ', minval(grid%ht_gc), maxval(grid%ht_gc)
484 CALL wrf_debug(100,message)
485 ENDIF
487 CALL boundary_smooth(grid%ht_gc,grid%landmask, grid, 12 , 12 &
488 &, IDS,IDE,JDS,JDE,KDS,KDE &
489 &, IMS,IME,JMS,JME,KMS,KME &
490 &, ITS,ITE,JTS,JTE,KTS,KTE )
492 DO j = jts, MIN(jte,jde-1)
493 DO i = its, MIN(ite,ide-1)
494 if (grid%landmask(I,J) .gt. 0.5) grid%sm(I,J)=0.
495 if (grid%landmask(I,J) .le. 0.5) grid%sm(I,J)=1.
496 if (grid%tsk_gc(I,J) .gt. 0.) then
497 grid%nmm_tsk(I,J)=grid%tsk_gc(I,J)
498 else
499 #if defined(HWRF)
500 if(grid%use_prep_hybrid) then
501 if(grid%t(I,J,1)<100) then
502 write(*,*) 'NO VALID SURFACE TEMPERATURE: I,J,TSK_GC(I,J),T(I,J,1) = ', &
503 I,J,grid%TSK_GC(I,J),grid%T(I,J,1)
504 else
505 grid%nmm_tsk(I,J)=grid%t(I,J,1) ! stopgap measure
506 end if
507 else
508 #endif
509 grid%nmm_tsk(I,J)=grid%t_gc(I,J,1) ! stopgap measure
510 #if defined(HWRF)
511 endif
512 #endif
513 endif
514 !
515 grid%glat(I,J)=grid%hlat_gc(I,J)*DEGRAD
516 grid%glon(I,J)=grid%hlon_gc(I,J)*DEGRAD
517 grid%weasd(I,J)=grid%snow(I,J)
518 grid%xice(I,J)=grid%xice_gc(I,J)
519 ENDDO
520 ENDDO
521 ! First item is to define the target vertical coordinate
523 num_metgrid_levels = grid%num_metgrid_levels
524 eta_levels(1:kde) = model_config_rec%eta_levels(1:kde)
525 ptsgm = model_config_rec%ptsgm
526 p_top_requested = grid%p_top_requested
527 grid%pt=p_top_requested
529 if (internal_time_loop .eq. 1) then
531 if (eta_levels(1) .ne. 1.0) then
532 #if defined(HWRF)
533 if(grid%use_prep_hybrid) then
534 call wrf_error_fatal('PREP_HYBRID ERROR: eta_levels is not specified, but use_prep_hybrid=.true.')
535 end if
536 #endif
538 write(message,*) '********************************************************************* '
539 CALL wrf_message(message)
540 write(message,*) '** eta_levels appears not to be specified in the namelist'
541 CALL wrf_message(message)
542 write(message,*) '** We will call compute_nmm_levels to define layer thicknesses.'
543 CALL wrf_message(message)
544 write(message,*) '** These levels should be reasonable for running the model, '
545 CALL wrf_message(message)
546 write(message,*) '** but may not be ideal for the simulation being made. Consider '
547 CALL wrf_message(message)
548 write(message,*) '** defining your own levels by specifying eta_levels in the model '
549 CALL wrf_message(message)
550 write(message,*) '** namelist. '
551 CALL wrf_message(message)
552 write(message,*) '********************************************************************** '
553 CALL wrf_message(message)
555 CALL compute_nmm_levels(KDE,p_top_requested,eta_levels)
557 DO L=1,KDE
558 write(message,*) 'L, eta_levels(L) returned :: ', L,eta_levels(L)
559 CALL wrf_message(message)
560 ENDDO
562 endif
564 write(message,*) 'KDE-1: ', KDE-1
565 CALL wrf_debug(1,message)
566 allocate(SG1(1:KDE-1))
567 allocate(SG2(1:KDE-1))
568 allocate(DSG1(1:KDE-1))
569 allocate(DSG2(1:KDE-1))
570 allocate(SGML1(1:KDE))
571 allocate(SGML2(1:KDE))
573 CALL define_nmm_vertical_coord (kde-1, ptsgm, grid%pt,grid%pdtop, eta_levels, &
574 grid%eta1,grid%deta1,grid%aeta1, &
575 grid%eta2,grid%deta2,grid%aeta2, grid%dfl, grid%dfrlg )
577 DO L=KDS,KDE-1
578 grid%deta(L)=eta_levels(L)-eta_levels(L+1)
579 ENDDO
580 endif
582 write(message,*) 'num_metgrid_levels: ', num_metgrid_levels
583 CALL wrf_message(message)
585 DO j = jts, MIN(jte,jde-1)
586 DO i = its, MIN(ite,ide-1)
587 grid%fis(I,J)=grid%ht_gc(I,J)*g
588 !
589 ! 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
590 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
591 IF (mod(I,10) .eq. 0 .and. mod(J,10) .eq. 0) THEN
592 write(message,*) 'grid%ht_gc and grid%toposoil to swap, flag_soilhgt ::: ', &
593 I,J, grid%ht_gc(I,J),grid%toposoil(I,J),flag_soilhgt
594 CALL wrf_debug(10,message)
595 ENDIF
596 IF ( ( flag_soilhgt.EQ. 1 ) ) THEN
597 grid%ght_gc(I,J,1)=grid%toposoil(I,J)
598 ENDIF
599 ENDIF
601 ENDDO
602 ENDDO
604 numzero=0
605 numexamined=0
606 DO j = jts, MIN(jte,jde-1)
607 DO i = its, MIN(ite,ide-1)
608 numexamined=numexamined+1
609 if(grid%fis(i,j)<1e-5 .and. grid%fis(i,j)>-1e5 ) then
610 numzero=numzero+1
611 end if
612 enddo
613 enddo
614 write(message,*) 'TOTAL NEAR-ZERO FIS POINTS: ',numzero,' OF ',numexamined
615 call wrf_debug(10,message)
616 #if defined(HWRF)
617 interp_notph: if(.not. grid%use_prep_hybrid) then
618 #endif
619 if (.NOT. allocated(PDVP)) allocate(PDVP(IMS:IME,JMS:JME))
620 if (.NOT. allocated(P3D_OUT)) allocate(P3D_OUT(IMS:IME,JMS:JME,KDS:KDE-1))
621 if (.NOT. allocated(PSFC_OUTV)) allocate(PSFC_OUTV(IMS:IME,JMS:JME))
622 if (.NOT. allocated(P3DV_OUT)) allocate(P3DV_OUT(IMS:IME,JMS:JME,KDS:KDE-1))
623 if (.NOT. allocated(P3DV_IN)) allocate(P3DV_IN(IMS:IME,JMS:JME,num_metgrid_levels))
625 CALL compute_nmm_surfacep (grid%ht_gc, grid%ght_gc, grid%p_gc , grid%t_gc &
626 &, grid%psfc_out, num_metgrid_levels &
627 &, IDS,IDE,JDS,JDE,KDS,KDE &
628 &, IMS,IME,JMS,JME,KMS,KME &
629 &, ITS,ITE,JTS,JTE,KTS,KTE ) ! H points
631 if (internal_time_loop .eq. 1) then
633 write(message,*) 'psfc points (final combined)'
634 loopinc=max( (JTE-JTS)/20,1)
635 iloopinc=max( (ITE-ITS)/10,1)
636 CALL wrf_message(message)
637 DO J=min(JTE,JDE-1),JTS,-loopinc
638 write(message,633) (grid%psfc_out(I,J)/100.,I=its,min(ite,IDE-1),iloopinc)
639 CALL wrf_message(message)
640 ENDDO
642 endif
644 633 format(35(f5.0,1x))
646 CALL compute_3d_pressure (grid%psfc_out,grid%aeta1,grid%aeta2 &
647 &, grid%pdtop,grid%pt,grid%pd,p3d_out &
648 &, IDS,IDE,JDS,JDE,KDS,KDE &
649 &, IMS,IME,JMS,JME,KMS,KME &
650 &, ITS,ITE,JTS,JTE,KTS,KTE )
652 #ifdef DM_PARALLEL
653 ips=its ; ipe=ite ; jps=jts ; jpe=jte ; kps=kts ; kpe=kte
654 # include "HALO_NMM_MG2.inc"
655 #endif
657 #ifdef DM_PARALLEL
658 # include "HALO_NMM_MG3.inc"
659 #endif
661 do K=1,num_metgrid_levels
662 do J=JTS,min(JTE,JDE-1)
663 do I=ITS,min(ITE,IDE-1)
665 IF (K .eq. KTS) THEN
666 IF (J .eq. JDS .and. I .lt. IDE-1) THEN ! S boundary
667 PDVP(I,J)=0.5*(grid%pd(I,J)+grid%pd(I+1,J))
668 PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J)+grid%psfc_out(I+1,J))
669 ELSEIF (J .eq. JDE-1 .and. I .lt. IDE-1) THEN ! N boundary
670 PDVP(I,J)=0.5*(grid%pd(I,J)+grid%pd(I+1,J))
671 PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J)+grid%psfc_out(I+1,J))
672 ELSEIF (I .eq. IDS .and. mod(J,2) .eq. 0) THEN ! W boundary
673 PDVP(I,J)=0.5*(grid%pd(I,J-1)+grid%pd(I,J+1))
674 PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J-1)+grid%psfc_out(I,J+1))
675 ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 0) THEN ! E boundary
676 PDVP(I,J)=0.5*(grid%pd(I,J-1)+grid%pd(I,J+1))
677 PSFC_OUTV(I,J)=0.5*(grid%psfc_out(I,J-1)+grid%psfc_out(I,J+1))
678 ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 1) THEN ! phantom E boundary
679 PDVP(I,J)=grid%pd(I,J)
680 PSFC_OUTV(I,J)=grid%psfc_out(I,J)
681 ELSEIF (mod(J,2) .eq. 0) THEN ! interior even row
682 PDVP(I,J)=0.25*(grid%pd(I,J)+grid%pd(I-1,J)+grid%pd(I,J+1)+grid%pd(I,J-1))
683 PSFC_OUTV(I,J)=0.25*(grid%psfc_out(I,J)+grid%psfc_out(I-1,J)+ &
684 grid%psfc_out(I,J+1)+grid%psfc_out(I,J-1))
685 ELSE ! interior odd row
686 PDVP(I,J)=0.25*(grid%pd(I,J)+grid%pd(I+1,J)+grid%pd(I,J+1)+grid%pd(I,J-1))
687 PSFC_OUTV(I,J)=0.25*(grid%psfc_out(I,J)+grid%psfc_out(I+1,J)+ &
688 grid%psfc_out(I,J+1)+grid%psfc_out(I,J-1))
689 ENDIF
690 ENDIF
692 IF (J .eq. JDS .and. I .lt. IDE-1) THEN ! S boundary
693 P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J,K)+grid%p_gc(I+1,J,K))
694 ELSEIF (J .eq. JDE-1 .and. I .lt. IDE-1) THEN ! N boundary
695 P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J,K)+grid%p_gc(I+1,J,K))
696 ELSEIF (I .eq. IDS .and. mod(J,2) .eq. 0) THEN ! W boundary
697 P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J-1,K)+grid%p_gc(I,J+1,K))
698 ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 0) THEN ! E boundary
699 P3DV_IN(I,J,K)=0.5*(grid%p_gc(I,J-1,K)+grid%p_gc(I,J+1,K))
700 ELSEIF (I .eq. IDE-1 .and. mod(J,2) .eq. 1) THEN ! phantom E boundary
701 P3DV_IN(I,J,K)=grid%p_gc(I,J,K)
702 ELSEIF (mod(J,2) .eq. 0) THEN ! interior even row
703 P3DV_IN(I,J,K)=0.25*(grid%p_gc(I,J,K)+grid%p_gc(I-1,J,K) + &
704 grid%p_gc(I,J+1,K)+grid%p_gc(I,J-1,K))
705 ELSE ! interior odd row
706 P3DV_IN(I,J,K)=0.25*(grid%p_gc(I,J,K)+grid%p_gc(I+1,J,K) + &
707 grid%p_gc(I,J+1,K)+grid%p_gc(I,J-1,K))
708 ENDIF
710 enddo
711 enddo
712 enddo
714 CALL compute_3d_pressure (psfc_outv,grid%aeta1,grid%aeta2 &
715 &, grid%pdtop,grid%pt,pdvp,p3dv_out &
716 &, IDS,IDE,JDS,JDE,KDS,KDE &
717 &, IMS,IME,JMS,JME,KMS,KME &
718 &, ITS,ITE,JTS,JTE,KTS,KTE )
720 CALL interp_press2press_lin(grid%p_gc, p3d_out &
721 &, grid%t_gc, grid%t,num_metgrid_levels &
722 &, .TRUE.,.TRUE.,.TRUE. & ! extrap, ignore_lowest, t_field
723 &, IDS,IDE,JDS,JDE,KDS,KDE &
724 &, IMS,IME,JMS,JME,KMS,KME &
725 &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
728 CALL interp_press2press_lin(p3dv_in, p3dv_out &
729 &, grid%u_gc, grid%u,num_metgrid_levels &
730 &, .FALSE.,.TRUE.,.FALSE. &
731 &, IDS,IDE,JDS,JDE,KDS,KDE &
732 &, IMS,IME,JMS,JME,KMS,KME &
733 &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
735 CALL interp_press2press_lin(p3dv_in, p3dv_out &
736 &, grid%v_gc, grid%v,num_metgrid_levels &
737 &, .FALSE.,.TRUE.,.FALSE. &
738 &, IDS,IDE,JDS,JDE,KDS,KDE &
739 &, IMS,IME,JMS,JME,KMS,KME &
740 &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
742 IF (hyb_coor) THEN
743 CALL wind_adjust(p3dv_in,p3dv_out,grid%u_gc,grid%v_gc,grid%u,grid%v &
744 &, num_metgrid_levels,5000. &
745 &, IDS,IDE,JDS,JDE,KDS,KDE &
746 &, IMS,IME,JMS,JME,KMS,KME &
747 &, ITS,ITE,JTS,JTE,KTS,KTE )
748 ENDIF
751 ALLOCATE(qtmp(IMS:IME,JMS:JME,num_metgrid_levels))
752 ALLOCATE(qtmp2(IMS:IME,JMS:JME,num_metgrid_levels))
754 CALL rh_to_mxrat (grid%rh_gc, grid%t_gc, grid%p_gc, qtmp , .TRUE. , &
755 ids , ide , jds , jde , 1 , num_metgrid_levels , &
756 ims , ime , jms , jme , 1 , num_metgrid_levels , &
757 its , ite , jts , jte , 1 , num_metgrid_levels )
759 do K=1,num_metgrid_levels
760 do J=JTS,min(JTE,JDE-1)
761 do I=ITS,min(ITE,IDE-1)
762 QTMP2(I,J,K)=QTMP(I,J,K)/(1.0+QTMP(I,J,K))
763 end do
764 end do
765 end do
767 CALL interp_press2press_log(grid%p_gc, p3d_out &
768 &, QTMP2, grid%q,num_metgrid_levels &
769 &, .FALSE.,.TRUE. &
770 &, IDS,IDE,JDS,JDE,KDS,KDE &
771 &, IMS,IME,JMS,JME,KMS,KME &
772 &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time_loop )
774 IF (ALLOCATED(QTMP)) DEALLOCATE(QTMP)
775 IF (ALLOCATED(QTMP)) DEALLOCATE(QTMP2)
776 #if defined(HWRF)
777 else ! we are using prep_hybrid
778 ! Compute surface pressure:
779 grid%psfc_out=grid%pdtop+grid%pd
780 end if interp_notph
781 #endif
783 ! Get the monthly values interpolated to the current date
784 ! for the traditional monthly
785 ! fields of green-ness fraction and background grid%albedo.
787 if (internal_time_loop .eq. 1 .or. config_flags%sst_update .eq. 1) then
789 CALL monthly_interp_to_date ( grid%greenfrac_gc , current_date , grid%vegfra , &
790 ids , ide , jds , jde , kds , kde , &
791 ims , ime , jms , jme , kms , kme , &
792 its , ite , jts , jte , kts , kte )
794 CALL monthly_interp_to_date ( grid%albedo12m_gc , current_date , grid%albbck , &
795 ids , ide , jds , jde , kds , kde , &
796 ims , ime , jms , jme , kms , kme , &
797 its , ite , jts , jte , kts , kte )
799 ! Get the min/max of each i,j for the monthly green-ness fraction.
801 CALL monthly_min_max ( grid%greenfrac_gc , grid%shdmin , grid%shdmax , &
802 ids , ide , jds , jde , kds , kde , &
803 ims , ime , jms , jme , kms , kme , &
804 its , ite , jts , jte , kts , kte )
806 ! The model expects the green-ness values in percent, not fraction.
808 DO j = jts, MIN(jte,jde-1)
809 DO i = its, MIN(ite,ide-1)
810 !! grid%vegfra(i,j) = grid%vegfra(i,j) * 100.
811 grid%shdmax(i,j) = grid%shdmax(i,j) * 100.
812 grid%shdmin(i,j) = grid%shdmin(i,j) * 100.
813 grid%vegfrc(I,J)=grid%vegfra(I,J)
814 END DO
815 END DO
817 ! The model expects the albedo fields as
818 ! a fraction, not a percent. Set the water values to 8%.
820 DO j = jts, MIN(jte,jde-1)
821 DO i = its, MIN(ite,ide-1)
822 if (grid%albbck(i,j) .lt. 5.) then
823 write(message,*) 'reset albedo to 8%... I,J,albbck:: ', I,J,grid%albbck(I,J)
824 CALL wrf_debug(10,message)
825 grid%albbck(I,J)=8.
826 endif
827 grid%albbck(i,j) = grid%albbck(i,j) / 100.
828 grid%snoalb(i,j) = grid%snoalb(i,j) / 100.
829 IF ( grid%landmask(i,j) .LT. 0.5 ) THEN
830 grid%albbck(i,j) = 0.08
831 grid%snoalb(i,j) = 0.08
832 END IF
833 grid%albase(i,j)=grid%albbck(i,j)
834 grid%mxsnal(i,j)=grid%snoalb(i,j)
835 END DO
836 END DO
838 endif
840 #if defined(HWRF)
841 if(.not.grid%use_prep_hybrid) then
842 #endif
843 ! new deallocs
844 DEALLOCATE(p3d_out,p3dv_out,p3dv_in)
845 #if defined(HWRF)
846 end if
847 #endif
849 END IF ! <----- END OF VERTICAL INTERPOLATION PART ---->
852 !! compute SST at each time if updating SST
853 if (config_flags%sst_update == 1) then
855 DO j = jts, MIN(jde-1,jte)
856 DO i = its, MIN(ide-1,ite)
858 if (grid%SM(I,J) .lt. 0.5) then
859 grid%SST(I,J)=0.
860 endif
862 if (grid%SM(I,J) .gt. 0.5) then
863 grid%SST(I,J)=grid%NMM_TSK(I,J)
864 grid%NMM_TSK(I,J)=0.
865 endif
867 IF ( (grid%NMM_TSK(I,J)+grid%SST(I,J)) .lt. 200. .or. &
868 (grid%NMM_TSK(I,J)+grid%SST(I,J)) .gt. 350. ) THEN
869 write(message,*) 'TSK, SST trouble at : ', I,J
870 CALL wrf_message(message)
871 write(message,*) 'SM, NMM_TSK,SST ', grid%SM(I,J),grid%NMM_TSK(I,J),grid%SST(I,J)
872 CALL wrf_message(message)
873 ENDIF
875 ENDDO
876 ENDDO
878 endif ! sst_update test
880 if (internal_time_loop .eq. 1) then
882 !!! weasd has "snow water equivalent" in mm
884 DO j = jts, MIN(jte,jde-1)
885 DO i = its, MIN(ite,ide-1)
887 IF(grid%sm(I,J).GT.0.9) THEN
889 IF (grid%xice(I,J) .gt. 0) then
890 grid%si(I,J)=1.0
891 ENDIF
893 ! SEA
894 grid%epsr(I,J)=.97
895 grid%embck(I,J)=.97
896 grid%gffc(I,J)=0.
897 grid%albedo(I,J)=.06
898 grid%albase(I,J)=.06
899 IF(grid%si (I,J).GT.0. ) THEN
900 ! SEA-ICE
901 grid%sm(I,J)=0.
902 grid%si(I,J)=0.
903 grid%sice(I,J)=1.
904 grid%gffc(I,J)=0. ! just leave zero as irrelevant
905 grid%albedo(I,J)=.60
906 grid%albase(I,J)=.60
907 ENDIF
908 ELSE
910 grid%si(I,J)=5.0*grid%weasd(I,J)/1000.
911 ! LAND
912 grid%epsr(I,J)=1.0
913 grid%embck(I,J)=1.0
914 grid%gffc(I,J)=0.0 ! just leave zero as irrelevant
915 grid%sice(I,J)=0.
916 grid%sno(I,J)=grid%si(I,J)*.20
917 ENDIF
918 ENDDO
919 ENDDO
921 ! DETERMINE grid%albedo OVER LAND
922 DO j = jts, MIN(jte,jde-1)
923 DO i = its, MIN(ite,ide-1)
924 IF(grid%sm(I,J).LT.0.9.AND.grid%sice(I,J).LT.0.9) THEN
925 ! SNOWFREE albedo
926 IF ( (grid%sno(I,J) .EQ. 0.0) .OR. &
927 (grid%albase(I,J) .GE. grid%mxsnal(I,J) ) ) THEN
928 grid%albedo(I,J) = grid%albase(I,J)
929 ELSE
930 ! MODIFY albedo IF SNOWCOVER:
931 ! BELOW SNOWDEPTH THRESHOLD...
932 IF (grid%sno(I,J) .LT. SNUP) THEN
933 RSNOW = grid%sno(I,J)/SNUP
934 SNOFAC = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP))
935 ! ABOVE SNOWDEPTH THRESHOLD...
936 ELSE
937 SNOFAC = 1.0
938 ENDIF
939 ! CALCULATE grid%albedo ACCOUNTING FOR SNOWDEPTH AND VGFRCK
940 grid%albedo(I,J) = grid%albase(I,J) &
941 + (1.0-grid%vegfra(I,J))*SNOFAC*(grid%mxsnal(I,J)-grid%albase(I,J))
942 ENDIF
943 END IF
944 grid%si(I,J)=5.0*grid%weasd(I,J)
945 grid%sno(I,J)=grid%weasd(I,J)
947 !! convert vegfra
948 grid%vegfra(I,J)=grid%vegfra(I,J)*100.
949 !
950 ENDDO
951 ENDDO
953 #ifdef DM_PARALLEL
955 ALLOCATE(SM_G(IDS:IDE,JDS:JDE),SICE_G(IDS:IDE,JDS:JDE))
957 CALL WRF_PATCH_TO_GLOBAL_REAL( grid%sice(IMS,JMS) &
958 &, SICE_G,grid%DOMDESC &
959 &, 'z','xy' &
960 &, IDS,IDE-1,JDS,JDE-1,1,1 &
961 &, IMS,IME,JMS,JME,1,1 &
962 &, ITS,ITE,JTS,JTE,1,1 )
964 CALL WRF_PATCH_TO_GLOBAL_REAL( grid%sm(IMS,JMS) &
965 &, SM_G,grid%DOMDESC &
966 &, 'z','xy' &
967 &, IDS,IDE-1,JDS,JDE-1,1,1 &
968 &, IMS,IME,JMS,JME,1,1 &
969 &, ITS,ITE,JTS,JTE,1,1 )
972 IF (WRF_DM_ON_MONITOR()) THEN
974 637 format(40(f3.0,1x))
976 allocate(IHE_G(JDS:JDE-1),IHW_G(JDS:JDE-1))
977 DO j = JDS, JDE-1
978 IHE_G(J)=MOD(J+1,2)
979 IHW_G(J)=IHE_G(J)-1
980 ENDDO
982 DO ITER=1,10
983 DO j = jds+1, (jde-1)-1
984 DO i = ids+1, (ide-1)-1
986 ! any sea ice around point in question?
988 IF (SM_G(I,J) .ge. 0.9) THEN
989 SEAICESUM=SICE_G(I+IHE_G(J),J+1)+SICE_G(I+IHW_G(J),J+1)+ &
990 SICE_G(I+IHE_G(J),J-1)+SICE_G(I+IHW_G(J),J-1)
991 IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN
993 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. &
994 (SICE_G(I+IHW_G(J),J+1).lt.0.1 .and. SM_G(I+IHW_G(J),J+1).lt.0.1) .OR. &
995 (SICE_G(I+IHE_G(J),J-1).lt.0.1 .and. SM_G(I+IHE_G(J),J-1).lt.0.1) .OR. &
996 (SICE_G(I+IHW_G(J),J-1).lt.0.1 .and. SM_G(I+IHW_G(J),J-1).lt.0.1)) THEN
998 ! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE
1000 write(message,*) 'making seaice (1): ', I,J
1001 CALL wrf_debug(100,message)
1002 SICE_G(I,J)=1.0
1003 SM_G(I,J)=0.
1005 ENDIF
1007 ELSEIF (SEAICESUM .ge. 3) THEN
1009 ! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE
1011 write(message,*) 'making seaice (2): ', I,J
1012 CALL wrf_debug(100,message)
1013 SICE_G(I,J)=1.0
1014 SM_G(I,J)=0.
1015 ENDIF
1017 ENDIF
1019 ENDDO
1020 ENDDO
1021 ENDDO
1023 ENDIF
1025 CALL WRF_GLOBAL_TO_PATCH_REAL( SICE_G, grid%sice &
1026 &, grid%DOMDESC &
1027 &, 'z','xy' &
1028 &, IDS,IDE-1,JDS,JDE-1,1,1 &
1029 &, IMS,IME,JMS,JME,1,1 &
1030 &, ITS,ITE,JTS,JTE,1,1 )
1032 CALL WRF_GLOBAL_TO_PATCH_REAL( SM_G,grid%sm &
1033 &, grid%DOMDESC &
1034 &, 'z','xy' &
1035 &, IDS,IDE-1,JDS,JDE-1,1,1 &
1036 &, IMS,IME,JMS,JME,1,1 &
1037 &, ITS,ITE,JTS,JTE,1,1 )
1039 IF (WRF_DM_ON_MONITOR()) THEN
1040 #if defined(HWRF)
1041 ! SM_G is still needed for the high-res grid
1042 #else
1043 DEALLOCATE(SM_G)
1044 #endif
1045 deallocate(SICE_G)
1046 DEALLOCATE(IHE_G,IHW_G)
1048 ENDIF
1050 ! write(message,*) 'revised sea ice on patch'
1051 ! CALL wrf_debug(100,message)
1052 ! DO J=JTE,JTS,-(((JTE-JTS)/25)+1)
1053 ! write(message,637) (grid%sice(I,J),I=ITS,ITE,ITE/20)
1054 ! CALL wrf_debug(100,message)
1055 ! END DO
1057 #else
1058 ! serial sea ice reprocessing
1059 allocate(IHE(JDS:JDE-1),IHW(JDS:JDE-1))
1061 DO j = jts, MIN(jte,jde-1)
1062 IHE(J)=MOD(J+1,2)
1063 IHW(J)=IHE(J)-1
1064 ENDDO
1066 DO ITER=1,10
1067 DO j = jts+1, MIN(jte,jde-1)-1
1068 DO i = its+1, MIN(ite,ide-1)-1
1070 ! any sea ice around point in question?
1072 IF (grid%sm(I,J) .gt. 0.9) THEN
1073 SEAICESUM=grid%sice(I+IHE(J),J+1)+grid%sice(I+IHW(J),J+1)+ &
1074 grid%sice(I+IHE(J),J-1)+grid%sice(I+IHW(J),J-1)
1075 IF (SEAICESUM .ge. 1. .and. SEAICESUM .lt. 3.) THEN
1076 IF ((grid%sice(I+IHE(J),J+1).lt.0.1 .and. grid%sm(I+IHE(J),J+1).lt.0.1) .OR. &
1077 (grid%sice(I+IHW(J),J+1).lt.0.1 .and. grid%sm(I+IHW(J),J+1).lt.0.1) .OR. &
1078 (grid%sice(I+IHE(J),J-1).lt.0.1 .and. grid%sm(I+IHE(J),J-1).lt.0.1) .OR. &
1079 (grid%sice(I+IHW(J),J-1).lt.0.1 .and. grid%sm(I+IHW(J),J-1).lt.0.1)) THEN
1081 ! HAVE SEA ICE AND A SURROUNDING LAND POINT - CONVERT TO SEA ICE
1082 grid%sice(I,J)=1.0
1083 grid%sm(I,J)=0.
1084 ENDIF
1085 ELSEIF (SEAICESUM .ge. 3) THEN
1086 ! WATER POINT SURROUNDED BY ICE - CONVERT TO SEA ICE
1087 grid%sice(I,J)=1.0
1088 grid%sm(I,J)=0.
1089 ENDIF
1090 ENDIF
1092 ENDDO
1093 ENDDO
1094 ENDDO
1096 DEALLOCATE(IHE,IHW)
1097 #endif
1099 ! this block meant to guarantee land/sea agreement between sm and landmask
1101 DO j = jts, MIN(jte,jde-1)
1102 DO i = its, MIN(ite,ide-1)
1104 IF (grid%sm(I,J) .gt. 0.5) THEN
1105 grid%landmask(I,J)=0.0
1106 ELSEIF (grid%sm(I,J) .lt. 0.5 .and. grid%sice(I,J) .gt. 0.9) then
1107 grid%landmask(I,J)=0.0
1108 ELSEIF (grid%sm(I,J) .lt. 0.5 .and. grid%sice(I,J) .lt. 0.1) then
1109 grid%landmask(I,J)=1.0
1110 ELSE
1111 write(message,*) 'missed point in grid%landmask definition ' , I,J
1112 CALL wrf_message(message)
1113 grid%landmask(I,J)=0.0
1114 ENDIF
1115 !
1116 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
1117 write(message,*) 'set grid%nmm_tsk to: ', grid%sst(I,J)
1118 CALL wrf_message(message)
1119 grid%nmm_tsk(I,J)=grid%sst(I,J)
1120 grid%sst(I,J)=0.
1121 endif
1123 ENDDO
1124 ENDDO
1126 ! For sf_surface_physics = 1, we want to use close to a 10 cm value
1127 ! for the bottom level of the soil temps.
1129 IF ( ( model_config_rec%sf_surface_physics(grid%id) .EQ. 1 ) .AND. &
1130 ( flag_st000010 .EQ. 1 ) ) THEN
1131 DO j = jts , MIN(jde-1,jte)
1132 DO i = its , MIN(ide-1,ite)
1133 grid%soiltb(i,j) = grid%st000010(i,j)
1134 END DO
1135 END DO
1136 END IF
1138 ! Adjust the various soil temperature values depending on the difference in
1139 ! in elevation between the current model's elevation and the incoming data's
1140 ! orography.
1142 IF ( ( flag_toposoil .EQ. 1 ) ) THEN
1144 ALLOCATE(HT(ims:ime,jms:jme))
1146 DO J=jms,jme
1147 DO I=ims,ime
1148 HT(I,J)=grid%fis(I,J)/9.81
1149 END DO
1150 END DO
1152 ! if (maxval(grid%toposoil) .gt. 100.) then
1153 !
1154 ! Being avoided. Something to revisit eventually.
1155 !
1156 !1219 might be simply a matter of including toposoil
1157 !
1158 ! CODE NOT TESTED AT NCEP USING THIS FUNCTIONALITY,
1159 ! SO TO BE SAFE WILL AVOID FOR RETRO RUNS.
1160 !
1161 ! CALL adjust_soil_temp_new ( grid%soiltb , 2 , &
1162 ! grid%nmm_tsk , ht , grid%toposoil , grid%landmask, flag_toposoil , &
1163 ! grid%st000010 , st010040 , st040100 , st100200 , st010200 , &
1164 ! flag_st000010 , flag_st010040 , flag_st040100 , &
1165 ! flag_st100200 , flag_st010200 , &
1166 ! soilt000 , soilt005 , soilt020 , soilt040 , &
1167 ! soilt160 , soilt300 , &
1168 ! flag_soilt000 , flag_soilt005 , flag_soilt020 , &
1169 ! flag_soilt040 , flag_soilt160 , flag_soilt300 , &
1170 ! ids , ide , jds , jde , kds , kde , &
1171 ! ims , ime , jms , jme , kms , kme , &
1172 ! its , ite , jts , jte , kts , kte )
1173 ! endif
1175 END IF
1177 ! Process the LSM data.
1179 ! surface_input_source=1 => use data from static file
1180 ! (fractional category as input)
1181 ! surface_input_source=2 => use data from grib file
1182 ! (dominant category as input)
1184 IF ( config_flags%surface_input_source .EQ. 1 ) THEN
1185 grid%vegcat (its,jts) = 0
1186 grid%soilcat(its,jts) = 0
1187 END IF
1189 ! Generate the vegetation and soil category information
1190 ! from the fractional input
1191 ! data, or use the existing dominant category fields if they exist.
1193 IF ((grid%soilcat(its,jts) .LT. 0.5) .AND. (grid%vegcat(its,jts) .LT. 0.5)) THEN
1195 num_veg_cat = SIZE ( grid%landusef_gc , DIM=3 )
1196 num_soil_top_cat = SIZE ( grid%soilctop_gc , DIM=3 )
1197 num_soil_bot_cat = SIZE ( grid%soilcbot_gc , DIM=3 )
1199 do J=JMS,JME
1200 do K=1,num_veg_cat
1201 do I=IMS,IME
1202 grid%landusef(I,K,J)=grid%landusef_gc(I,J,K)
1203 enddo
1204 enddo
1205 enddo
1207 do J=JMS,JME
1208 do K=1,num_soil_top_cat
1209 do I=IMS,IME
1210 grid%soilctop(I,K,J)=grid%soilctop_gc(I,J,K)
1211 enddo
1212 enddo
1213 enddo
1215 do J=JMS,JME
1216 do K=1,num_soil_bot_cat
1217 do I=IMS,IME
1218 grid%soilcbot(I,K,J)=grid%soilcbot_gc(I,J,K)
1219 enddo
1220 enddo
1221 enddo
1223 ! grid%sm (1=water, 0=land)
1224 ! grid%landmask(0=water, 1=land)
1227 write(message,*) 'landmask into process_percent_cat_new'
1229 CALL wrf_debug(1,message)
1230 do J=JTE,JTS,-(((JTE-JTS)/20)+1)
1231 write(message,641) (grid%landmask(I,J),I=ITS,min(ITE,IDE-1),((ITE-ITS)/15)+1)
1232 CALL wrf_debug(1,message)
1233 enddo
1234 641 format(25(f3.0,1x))
1236 CALL process_percent_cat_new ( grid%landmask , &
1237 grid%landusef , grid%soilctop , grid%soilcbot , &
1238 grid%isltyp , grid%ivgtyp , &
1239 num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1240 ids , ide , jds , jde , kds , kde , &
1241 ims , ime , jms , jme , kms , kme , &
1242 its , ite , jts , jte , kts , kte , &
1243 model_config_rec%iswater(grid%id) )
1245 DO j = jts , MIN(jde-1,jte)
1246 DO i = its , MIN(ide-1,ite)
1247 grid%vegcat(i,j) = grid%ivgtyp(i,j)
1248 grid%soilcat(i,j) = grid%isltyp(i,j)
1249 END DO
1250 END DO
1252 ELSE
1254 ! Do we have dominant soil and veg data from the input already?
1256 IF ( grid%soilcat(its,jts) .GT. 0.5 ) THEN
1257 DO j = jts, MIN(jde-1,jte)
1258 DO i = its, MIN(ide-1,ite)
1259 grid%isltyp(i,j) = NINT( grid%soilcat(i,j) )
1260 END DO
1261 END DO
1262 END IF
1263 IF ( grid%vegcat(its,jts) .GT. 0.5 ) THEN
1264 DO j = jts, MIN(jde-1,jte)
1265 DO i = its, MIN(ide-1,ite)
1266 grid%ivgtyp(i,j) = NINT( grid%vegcat(i,j) )
1267 END DO
1268 END DO
1269 END IF
1271 ENDIF
1273 DO j = jts, MIN(jde-1,jte)
1274 DO i = its, MIN(ide-1,ite)
1276 IF (grid%sice(I,J) .lt. 0.1) THEN
1277 IF (grid%landmask(I,J) .gt. 0.5 .and. grid%sm(I,J) .gt. 0.5) THEN
1278 write(message,*) 'land mask and grid%sm both > 0.5: ', &
1279 I,J,grid%landmask(I,J),grid%sm(I,J)
1280 CALL wrf_message(message)
1281 grid%sm(I,J)=0.
1282 ELSEIF (grid%landmask(I,J) .lt. 0.5 .and. grid%sm(I,J) .lt. 0.5) THEN
1283 write(message,*) 'land mask and grid%sm both < 0.5: ', &
1284 I,J, grid%landmask(I,J),grid%sm(I,J)
1285 CALL wrf_message(message)
1286 grid%sm(I,J)=1.
1287 ENDIF
1288 ELSE
1289 IF (grid%landmask(I,J) .gt. 0.5 .and. grid%sm(I,J)+grid%sice(I,J) .gt. 0.9) then
1290 write(message,*) 'landmask says LAND, sm/sice say SEAICE: ', I,J
1291 ENDIF
1292 ENDIF
1294 ENDDO
1295 ENDDO
1297 DO j = jts, MIN(jde-1,jte)
1298 DO i = its, MIN(ide-1,ite)
1300 if (grid%sice(I,J) .gt. 0.9) then
1301 grid%isltyp(I,J)=16
1302 grid%ivgtyp(I,J)=24
1303 endif
1305 ENDDO
1306 ENDDO
1308 DO j = jts, MIN(jde-1,jte)
1309 DO i = its, MIN(ide-1,ite)
1311 if (grid%sm(I,J) .lt. 0.5) then
1312 grid%sst(I,J)=0.
1313 endif
1315 if (grid%sm(I,J) .gt. 0.5) then
1316 if (grid%sst(I,J) .lt. 0.1) then
1317 grid%sst(I,J)=grid%nmm_tsk(I,J)
1318 endif
1319 grid%nmm_tsk(I,J)=0.
1320 endif
1322 IF ( (grid%nmm_tsk(I,J)+grid%sst(I,J)) .lt. 200. .or. &
1323 (grid%nmm_tsk(I,J)+grid%sst(I,J)) .gt. 350. ) THEN
1324 write(message,*) 'TSK, sst trouble at : ', I,J
1325 CALL wrf_message(message)
1326 write(message,*) 'sm, nmm_tsk,sst ', grid%sm(I,J),grid%nmm_tsk(I,J),grid%sst(I,J)
1327 CALL wrf_message(message)
1328 ENDIF
1330 ENDDO
1331 ENDDO
1333 write(message,*) 'grid%sm'
1334 CALL wrf_message(message)
1336 DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
1337 write(message,635) (grid%sm(i,J),I=its,ite,((ite-its)/10)+1)
1338 CALL wrf_message(message)
1339 END DO
1341 write(message,*) 'sst/nmm_tsk'
1342 CALL wrf_debug(10,message)
1343 DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
1344 write(message,635) (grid%sst(I,J)+grid%nmm_tsk(I,J),I=ITS,min(ide-1,ite),((ite-its)/10)+1)
1345 CALL wrf_debug(10,message)
1346 END DO
1348 635 format(20(f5.1,1x))
1350 DO j = jts, MIN(jde-1,jte)
1351 DO i = its, MIN(ide-1,ite)
1352 IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
1353 grid%soiltb(i,j) = grid%sst(i,j)
1354 ELSE IF ( grid%landmask(i,j) .GT. 0.5 ) THEN
1355 grid%soiltb(i,j) = grid%nmm_tsk(i,j)
1356 END IF
1357 END DO
1358 END DO
1360 ! END IF
1362 ! Land use categories, dominant soil and vegetation types (if available).
1364 ! allocate(grid%lu_index(ims:ime,jms:jme))
1366 DO j = jts, MIN(jde-1,jte)
1367 DO i = its, MIN(ide-1,ite)
1368 grid%lu_index(i,j) = grid%ivgtyp(i,j)
1369 END DO
1370 END DO
1372 if (flag_sst .eq. 1) log_flag_sst=.true.
1373 if (flag_sst .eq. 0) log_flag_sst=.false.
1375 write(message,*) 'st_input dimensions: ', size(st_input,dim=1), &
1376 size(st_input,dim=2),size(st_input,dim=3)
1377 CALL wrf_debug(100,message)
1379 ! write(message,*) 'maxval st_input(1): ', maxval(st_input(:,1,:))
1380 ! CALL wrf_message(message)
1381 ! write(message,*) 'maxval st_input(2): ', maxval(st_input(:,2,:))
1382 ! CALL wrf_message(message)
1383 ! write(message,*) 'maxval st_input(3): ', maxval(st_input(:,3,:))
1384 ! CALL wrf_message(message)
1385 ! write(message,*) 'maxval st_input(4): ', maxval(st_input(:,4,:))
1386 ! CALL wrf_message(message)
1388 ! =============================================================
1390 IF (.NOT. ALLOCATED(TG_ALT))ALLOCATE(TG_ALT(grid%sm31:grid%em31,jms:jme))
1392 TPH0=TPH0D*DTR
1393 WBD=-(((ide-1)-1)*grid%dlmd)
1394 WB= WBD*DTR
1395 SBD=-(((jde-1)/2)*grid%dphd)
1396 SB= SBD*DTR
1397 DLM=grid%dlmd*DTR
1398 DPH=grid%dphd*DTR
1399 TDLM=DLM+DLM
1400 TDPH=DPH+DPH
1401 WBI=WB+TDLM
1402 SBI=SB+TDPH
1403 EBI=WB+(ide-2)*TDLM
1404 ANBI=SB+(jde-2)*DPH
1405 STPH0=SIN(TPH0)
1406 CTPH0=COS(TPH0)
1407 TSPH=3600./GRID%DT
1408 DO J=JTS,min(JTE,JDE-1)
1409 TLM=WB-TDLM+MOD(J,2)*DLM !For velocity points on the E grid
1410 TPH=SB+float(J-1)*DPH
1411 STPH=SIN(TPH)
1412 CTPH=COS(TPH)
1413 DO I=ITS,MIN(ITE,IDE-1)
1415 if (I .eq. ITS) THEN
1416 TLM=TLM+TDLM*ITS
1417 else
1418 TLM=TLM+TDLM
1419 endif
1421 TERM1=(STPH0*CTPH*COS(TLM)+CTPH0*STPH)
1422 FP=TWOM*(TERM1)
1423 grid%f(I,J)=0.5*GRID%DT*FP
1424 ENDDO
1425 ENDDO
1426 DO J=JTS,min(JTE,JDE-1)
1427 TLM=WB-TDLM+MOD(J+1,2)*DLM !For mass points on the E grid
1428 TPH=SB+float(J-1)*DPH
1429 STPH=SIN(TPH)
1430 CTPH=COS(TPH)
1431 DO I=ITS,MIN(ITE,IDE-1)
1433 if (I .eq. ITS) THEN
1434 TLM=TLM+TDLM*ITS
1435 else
1436 TLM=TLM+TDLM
1437 endif
1439 TERM1=(STPH0*CTPH*COS(TLM)+CTPH0*STPH)
1440 TERM1=MIN(TERM1,1.0D0)
1441 TERM1=MAX(TERM1,-1.0D0)
1442 APH=ASIN(TERM1)
1443 TG_ALT(I,J)=TG0+TGA*COS(APH)-grid%fis(I,J)/3333.
1444 ENDDO
1445 ENDDO
1447 DO j = jts, MIN(jde-1,jte)
1448 DO i = its, MIN(ide-1,ite)
1449 ! IF ( ( grid%landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) .AND. &
1450 ! grid%sice(I,J) .eq. 0. ) THEN
1451 ! grid%tg(i,j) = grid%sst(i,j)
1452 ! ELSEIF (grid%sice(I,J) .eq. 1) THEN
1453 ! grid%tg(i,j) = 271.16
1454 ! END IF
1456 if (grid%tg(I,J) .lt. 200.) then ! only use default TG_ALT definition if
1457 ! not getting TGROUND from grid%si
1458 grid%tg(I,J)=TG_ALT(I,J)
1459 endif
1461 if (grid%tg(I,J) .lt. 200. .or. grid%tg(I,J) .gt. 320.) then
1462 write(message,*) 'problematic grid%tg point at : ', I,J
1463 CALL wrf_message( message )
1464 endif
1466 adum2d(i,j)=grid%nmm_tsk(I,J)+grid%sst(I,J)
1468 END DO
1469 END DO
1471 DEALLOCATE(TG_ALT)
1473 write(message,*) 'call process_soil_real with num_st_levels_input: ', num_st_levels_input
1474 CALL wrf_message( message )
1476 ! =============================================================
1478 CALL process_soil_real ( adum2d, grid%tg , &
1479 grid%landmask, grid%sst, &
1480 st_input, sm_input, sw_input, &
1481 st_levels_input , sm_levels_input , &
1482 sw_levels_input , &
1483 grid%sldpth , grid%dzsoil , grid%stc , grid%smc , grid%sh2o, &
1484 flag_sst , flag_soilt000, flag_soilm000, &
1485 ids , ide , jds , jde , kds , kde , &
1486 ims , ime , jms , jme , kms , kme , &
1487 its , ite , jts , jte , kts , kte , &
1488 model_config_rec%sf_surface_physics(grid%id) , &
1489 model_config_rec%num_soil_layers , &
1490 model_config_rec%real_data_init_type , &
1491 num_st_levels_input , num_sm_levels_input , &
1492 num_sw_levels_input , &
1493 num_st_levels_alloc , num_sm_levels_alloc , &
1494 num_sw_levels_alloc )
1496 ! =============================================================
1498 ! Minimum soil values, residual, from RUC LSM scheme.
1499 ! For input from Noah and using
1500 ! RUC LSM scheme, this must be subtracted from the input
1501 ! total soil moisture. For input RUC data and using the Noah LSM scheme,
1502 ! this value must be added to the soil moisture_input.
1504 lqmi(1:num_soil_top_cat) = &
1505 (/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
1506 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
1507 0.004, 0.065 /) !dusan , 0.020, 0.004, 0.008 /)
1509 ! At the initial time we care about values of soil moisture and temperature,
1510 ! other times are ignored by the model, so we ignore them, too.
1512 account_for_zero_soil_moisture : SELECT CASE ( model_config_rec%sf_surface_physics(grid%id) )
1514 CASE ( LSMSCHEME )
1515 iicount = 0
1516 IF ( FLAG_SM000010 .EQ. 1 ) THEN
1517 DO j = jts, MIN(jde-1,jte)
1518 DO i = its, MIN(ide-1,ite)
1519 IF ((grid%landmask(i,j).gt.0.5) .and. (grid%stc(i,1,j) .gt. 200) .and. &
1520 (grid%stc(i,1,j) .lt. 400) .and. (grid%smc(i,1,j) .lt. 0.005)) then
1521 write(message,*) 'Noah > Noah: bad soil moisture at i,j = ',i,j,grid%smc(i,:,j)
1522 CALL wrf_message(message)
1523 iicount = iicount + 1
1524 grid%smc(i,:,j) = 0.005
1525 END IF
1526 END DO
1527 END DO
1528 IF ( iicount .GT. 0 ) THEN
1529 write(message,*) 'Noah -> Noah: total number of small soil moisture locations= ',&
1530 iicount
1531 CALL wrf_message(message)
1532 END IF
1533 ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
1534 DO j = jts, MIN(jde-1,jte)
1535 DO i = its, MIN(ide-1,ite)
1536 grid%smc(i,:,j) = grid%smc(i,:,j) + lqmi(grid%isltyp(i,j))
1537 END DO
1538 END DO
1539 DO j = jts, MIN(jde-1,jte)
1540 DO i = its, MIN(ide-1,ite)
1541 IF ((grid%landmask(i,j).gt.0.5) .and. (grid%stc(i,1,j) .gt. 200) .and. &
1542 (grid%stc(i,1,j) .lt. 400) .and. (grid%smc(i,1,j) .lt. 0.004)) then
1543 write(message,*) 'RUC -> Noah: bad soil moisture at i,j = ' &
1544 ,i,j,grid%smc(i,:,j)
1545 CALL wrf_message(message)
1546 iicount = iicount + 1
1547 grid%smc(i,:,j) = 0.004
1548 END IF
1549 END DO
1550 END DO
1551 IF ( iicount .GT. 0 ) THEN
1552 write(message,*) 'RUC -> Noah: total number of small soil moisture locations = ',&
1553 iicount
1554 CALL wrf_message(message)
1555 END IF
1556 END IF
1557 CASE ( RUCLSMSCHEME )
1558 iicount = 0
1559 IF ( FLAG_SM000010 .EQ. 1 ) THEN
1560 DO j = jts, MIN(jde-1,jte)
1561 DO i = its, MIN(ide-1,ite)
1562 grid%smc(i,:,j) = MAX ( grid%smc(i,:,j) - lqmi(grid%isltyp(i,j)) , 0. )
1563 END DO
1564 END DO
1565 ELSE IF ( FLAG_SOILM000 .EQ. 1 ) THEN
1566 ! no op
1567 END IF
1569 END SELECT account_for_zero_soil_moisture
1571 !!! zero out grid%nmm_tsk at water points again
1573 DO j = jts, MIN(jde-1,jte)
1574 DO i = its, MIN(ide-1,ite)
1575 if (grid%sm(I,J) .gt. 0.5) then
1576 grid%nmm_tsk(I,J)=0.
1577 endif
1578 END DO
1579 END DO
1581 !! check on grid%stc
1583 DO j = jts, MIN(jde-1,jte)
1584 DO i = its, MIN(ide-1,ite)
1586 IF (grid%sice(I,J) .gt. 0.9) then
1587 DO L = 1, grid%num_soil_layers
1588 grid%stc(I,L,J)=271.16 ! grid%tg value used by Eta/NMM
1589 END DO
1590 END IF
1592 IF (grid%sm(I,J) .gt. 0.9) then
1593 DO L = 1, grid%num_soil_layers
1594 grid%stc(I,L,J)=273.16 ! grid%tg value used by Eta/NMM
1595 END DO
1596 END IF
1598 END DO
1599 END DO
1601 DO j = jts, MIN(jde-1,jte)
1602 DO i = its, MIN(ide-1,ite)
1604 if (grid%sm(I,J) .lt. 0.1 .and. grid%stc(I,1,J) .lt. 0.1) THEN
1605 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)
1606 CALL wrf_message(message)
1607 do JJ=J-1,J+1
1608 do L=1, grid%num_soil_layers
1609 do II=I-1,I+1
1611 if (II .ge. its .and. II .le. MIN(ide-1,ite) .and. &
1612 JJ .ge. jts .and. JJ .le. MIN(jde-1,jte)) then
1614 grid%stc(I,L,J)=amax1(grid%stc(I,L,J),grid%stc(II,L,JJ))
1615 cur_smc=grid%smc(I,L,J)
1617 if ( grid%smc(II,L,JJ) .gt. 0.005 .and. grid%smc(II,L,JJ) .lt. 1.0) then
1618 aposs_smc=grid%smc(II,L,JJ)
1620 if ( cur_smc .eq. 0 ) then
1621 cur_smc=aposs_smc
1622 grid%smc(I,L,J)=cur_smc
1623 else
1624 cur_smc=amin1(cur_smc,aposs_smc)
1625 cur_smc=amin1(cur_smc,aposs_smc)
1626 grid%smc(I,L,J)=cur_smc
1627 endif
1628 endif
1630 endif ! bounds check
1632 enddo
1633 enddo
1634 enddo
1635 write(message,*) 'grid%stc, grid%smc(1) now: ', grid%stc(I,1,J),grid%smc(I,1,J)
1636 CALL wrf_message(message)
1637 endif
1639 if (grid%stc(I,1,J) .lt. 0.1) then
1640 write(message,*) 'QUITTING DUE TO STILL troublesome grid%stc value: ', I,J, grid%stc(I,1,J),grid%smc(I,1,J)
1641 call wrf_error_fatal(message)
1642 endif
1644 ENDDO
1645 ENDDO
1647 !hardwire soil stuff for time being
1649 ! RTDPTH=0.
1650 ! RTDPTH(1)=0.1
1651 ! RTDPTH(2)=0.3
1652 ! RTDPTH(3)=0.6
1654 ! grid%sldpth=0.
1655 ! grid%sldpth(1)=0.1
1656 ! grid%sldpth(2)=0.3
1657 ! grid%sldpth(3)=0.6
1658 ! grid%sldpth(4)=1.0
1660 !!! main body of nmm_specific starts here
1661 !
1662 do J=jts,min(jte,jde-1)
1663 do I=its,min(ite,ide-1)
1664 grid%res(I,J)=1.
1665 enddo
1666 enddo
1668 !! grid%hbm2
1670 grid%hbm2=0.
1672 do J=jts,min(jte,jde-1)
1673 do I=its,min(ite,ide-1)
1675 IF ( (J .ge. 3 .and. J .le. (jde-1)-2) .AND. &
1676 (I .ge. 2 .and. I .le. (ide-1)-2+mod(J,2)) ) THEN
1677 grid%hbm2(I,J)=1.
1678 ENDIF
1679 enddo
1680 enddo
1682 !! grid%hbm3
1683 grid%hbm3=0.
1685 !! LOOP OVER LOCAL DIMENSIONS
1687 do J=jts,min(jte,jde-1)
1688 grid%ihwg(J)=mod(J+1,2)-1
1689 IF (J .ge. 4 .and. J .le. (jde-1)-3) THEN
1690 IHL=(ids+1)-grid%ihwg(J)
1691 IHH=(ide-1)-2
1692 do I=its,min(ite,ide-1)
1693 IF (I .ge. IHL .and. I .le. IHH) grid%hbm3(I,J)=1.
1694 enddo
1695 ENDIF
1696 enddo
1698 !! grid%vbm2
1700 grid%vbm2=0.
1702 do J=jts,min(jte,jde-1)
1703 do I=its,min(ite,ide-1)
1705 IF ( (J .ge. 3 .and. J .le. (jde-1)-2) .AND. &
1706 (I .ge. 2 .and. I .le. (ide-1)-1-mod(J,2)) ) THEN
1708 grid%vbm2(I,J)=1.
1710 ENDIF
1712 enddo
1713 enddo
1715 !! grid%vbm3
1717 grid%vbm3=0.
1719 do J=jts,min(jte,jde-1)
1720 do I=its,min(ite,ide-1)
1722 IF ( (J .ge. 4 .and. J .le. (jde-1)-3) .AND. &
1723 (I .ge. 3-mod(J,2) .and. I .le. (ide-1)-2) ) THEN
1724 grid%vbm3(I,J)=1.
1725 ENDIF
1727 enddo
1728 enddo
1730 COAC=model_config_rec%coac(grid%id)
1731 CODAMP=model_config_rec%codamp(grid%id)
1733 DTAD=1.0
1734 ! IDTCF=DTCF, IDTCF=4
1735 DTCF=4.0 ! used?
1737 grid%dy_nmm=ERAD*DPH
1738 grid%cpgfv=-GRID%DT/(48.*grid%dy_nmm)
1739 grid%en= GRID%DT/( 4.*grid%dy_nmm)*DTAD
1740 grid%ent=GRID%DT/(16.*grid%dy_nmm)*DTAD
1742 DO J=jts,nnyp
1743 KHL2(J)=(IDE-1)*(J-1)-(J-1)/2+2
1744 KVL2(J)=(IDE-1)*(J-1)-J/2+2
1745 KHH2(J)=(IDE-1)*J-J/2-1
1746 KVH2(J)=(IDE-1)*J-(J+1)/2-1
1747 ENDDO
1749 TPH=SB-DPH
1751 DO J=jts,min(jte,jde-1)
1752 TPH=SB+float(J-1)*DPH
1753 DXP=ERAD*DLM*COS(TPH)
1754 DXJ(J)=DXP
1755 WPDARJ(J)=-W_NMM * &
1756 ((ERAD*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+grid%dy_nmm**2)/ &
1757 (GRID%DT*32.*DXP*grid%dy_nmm)
1759 CPGFUJ(J)=-GRID%DT/(48.*DXP)
1760 CURVJ(J)=.5*GRID%DT*TAN(TPH)/ERAD
1761 FCPJ(J)=GRID%DT/(CP*192.*DXP*grid%dy_nmm)
1762 FDIVJ(J)=1./(12.*DXP*grid%dy_nmm)
1763 ! EMJ(J)= GRID%DT/( 4.*DXP)*DTAD
1764 ! EMTJ(J)=GRID%DT/(16.*DXP)*DTAD
1765 FADJ(J)=-GRID%DT/(48.*DXP*grid%dy_nmm)*DTAD
1766 ACDT=GRID%DT*SQRT((ERAD*DLM*AMIN1(COS(ANBI),COS(SBI)))**2+grid%dy_nmm**2)
1767 CDDAMP=CODAMP*ACDT
1768 HDACJ(J)=COAC*ACDT/(4.*DXP*grid%dy_nmm)
1769 DDMPUJ(J)=CDDAMP/DXP
1770 DDMPVJ(J)=CDDAMP/grid%dy_nmm
1771 ENDDO
1773 DO J=JTS,min(JTE,JDE-1)
1774 TLM=WB-TDLM+MOD(J,2)*DLM
1775 TPH=SB+float(J-1)*DPH
1776 STPH=SIN(TPH)
1777 CTPH=COS(TPH)
1778 DO I=ITS,MIN(ITE,IDE-1)
1780 if (I .eq. ITS) THEN
1781 TLM=TLM+TDLM*ITS
1782 else
1783 TLM=TLM+TDLM
1784 endif
1786 FP=TWOM*(CTPH0*STPH+STPH0*CTPH*COS(TLM))
1787 grid%f(I,J)=0.5*GRID%DT*FP
1789 ENDDO
1790 ENDDO
1792 ! --------------DERIVED VERTICAL GRID CONSTANTS--------------------------
1794 grid%ef4t=.5*GRID%DT/CP
1795 grid%f4q = -GRID%DT*DTAD
1796 grid%f4d =-.5*GRID%DT*DTAD
1798 DO L=KDS,KDE-1
1799 grid%rdeta(L)=1./grid%deta(L)
1800 grid%f4q2(L)=-.25*GRID%DT*DTAD/grid%deta(L)
1801 ENDDO
1803 DO J=JTS,min(JTE,JDE-1)
1804 DO I=ITS,min(ITE,IDE-1)
1805 grid%dx_nmm(I,J)=DXJ(J)
1806 grid%wpdar(I,J)=WPDARJ(J)*grid%hbm2(I,J)
1807 grid%cpgfu(I,J)=CPGFUJ(J)*grid%vbm2(I,J)
1808 grid%curv(I,J)=CURVJ(J)*grid%vbm2(I,J)
1809 grid%fcp(I,J)=FCPJ(J)*grid%hbm2(I,J)
1810 grid%fdiv(I,J)=FDIVJ(J)*grid%hbm2(I,J)
1811 grid%fad(I,J)=FADJ(J)
1812 grid%hdacv(I,J)=HDACJ(J)*grid%vbm2(I,J)
1813 grid%hdac(I,J)=HDACJ(J)*1.25*grid%hbm2(I,J)
1814 ENDDO
1815 ENDDO
1817 DO J=JTS, MIN(JDE-1,JTE)
1819 IF (J.LE.5.OR.J.GE.(JDE-1)-4) THEN
1821 KHH=(IDE-1)-2+MOD(J,2) ! KHH is global...loop over I that have
1822 DO I=ITS,MIN(IDE-1,ITE)
1823 IF (I .ge. 2 .and. I .le. KHH) THEN
1824 grid%hdac(I,J)=grid%hdac(I,J)* DFC
1825 ENDIF
1826 ENDDO
1828 ELSE
1830 KHH=2+MOD(J,2)
1831 DO I=ITS,MIN(IDE-1,ITE)
1832 IF (I .ge. 2 .and. I .le. KHH) THEN
1833 grid%hdac(I,J)=grid%hdac(I,J)* DFC
1834 ENDIF
1835 ENDDO
1837 KHH=(IDE-1)-2+MOD(J,2)
1839 DO I=ITS,MIN(IDE-1,ITE)
1840 IF (I .ge. (IDE-1)-2 .and. I .le. KHH) THEN
1841 grid%hdac(I,J)=grid%hdac(I,J)* DFC
1842 ENDIF
1843 ENDDO
1844 ENDIF
1845 ENDDO
1847 DO J=JTS,min(JTE,JDE-1)
1848 DO I=ITS,min(ITE,IDE-1)
1849 grid%ddmpu(I,J)=DDMPUJ(J)*grid%vbm2(I,J)
1850 grid%ddmpv(I,J)=DDMPVJ(J)*grid%vbm2(I,J)
1851 grid%hdacv(I,J)=grid%hdacv(I,J)*grid%vbm2(I,J)
1852 ENDDO
1853 ENDDO
1854 ! --------------INCREASING DIFFUSION ALONG THE BOUNDARIES----------------
1856 DO J=JTS,MIN(JDE-1,JTE)
1857 IF (J.LE.5.OR.J.GE.JDE-1-4) THEN
1858 KVH=(IDE-1)-1-MOD(J,2)
1859 DO I=ITS,min(IDE-1,ITE)
1860 IF (I .ge. 2 .and. I .le. KVH) THEN
1861 grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC
1862 grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC
1863 grid%hdacv(I,J)=grid%hdacv(I,J)* DFC
1864 ENDIF
1865 ENDDO
1866 ELSE
1867 KVH=3-MOD(J,2)
1868 DO I=ITS,min(IDE-1,ITE)
1869 IF (I .ge. 2 .and. I .le. KVH) THEN
1870 grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC
1871 grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC
1872 grid%hdacv(I,J)=grid%hdacv(I,J)* DFC
1873 ENDIF
1874 ENDDO
1875 KVH=(IDE-1)-1-MOD(J,2)
1876 DO I=ITS,min(IDE-1,ITE)
1877 IF (I .ge. IDE-1-2 .and. I .le. KVH) THEN
1878 grid%ddmpu(I,J)=grid%ddmpu(I,J)*DDFC
1879 grid%ddmpv(I,J)=grid%ddmpv(I,J)*DDFC
1880 grid%hdacv(I,J)=grid%hdacv(I,J)* DFC
1881 ENDIF
1882 ENDDO
1883 ENDIF
1884 ENDDO
1886 write(message,*) 'grid%stc(1)'
1887 CALL wrf_message(message)
1888 DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
1889 write(message,635) (grid%stc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12+1)
1890 CALL wrf_message(message)
1891 ENDDO
1893 write(message,*) 'grid%smc(1)'
1894 CALL wrf_message(message)
1895 DO J=min(jde-1,jte),jts,-((jte-jts)/15+1)
1896 write(message,635) (grid%smc(I,1,J),I=its,min(ite,ide-1),(ite-its)/12+1)
1897 CALL wrf_message(message)
1898 ENDDO
1900 DO j = jts, MIN(jde-1,jte)
1901 DO i= ITS, MIN(IDE-1,ITE)
1903 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
1904 write(message,*) 'very moist on land point: ', I,J,grid%smc(I,1,J)
1905 CALL wrf_debug(10,message)
1906 endif
1908 enddo
1909 enddo
1911 !!! compute grid%emt, grid%em on global domain, and only on task 0.
1913 #ifdef DM_PARALLEL
1914 IF (wrf_dm_on_monitor()) THEN !!!! NECESSARY TO LIMIT THIS TO TASK ZERO?
1915 #else
1916 IF (JDS .eq. JTS) THEN !! set unfailable condition for serial job
1917 #endif
1919 ALLOCATE(EMJ(JDS:JDE-1),EMTJ(JDS:JDE-1))
1921 DO J=JDS,JDE-1
1922 TPH=SB+float(J-1)*DPH
1923 DXP=ERAD*DLM*COS(TPH)
1924 EMJ(J)= GRID%DT/( 4.*DXP)*DTAD
1925 EMTJ(J)=GRID%DT/(16.*DXP)*DTAD
1926 ENDDO
1928 JA=0
1929 DO 161 J=3,5
1930 JA=JA+1
1931 KHLA(JA)=2
1932 KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
1933 161 grid%emt(JA)=EMTJ(J)
1934 DO 162 J=(JDE-1)-4,(JDE-1)-2
1935 JA=JA+1
1936 KHLA(JA)=2
1937 KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
1938 162 grid%emt(JA)=EMTJ(J)
1939 DO 163 J=6,(JDE-1)-5
1940 JA=JA+1
1941 KHLA(JA)=2
1942 KHHA(JA)=2+MOD(J,2)
1943 163 grid%emt(JA)=EMTJ(J)
1944 DO 164 J=6,(JDE-1)-5
1945 JA=JA+1
1946 KHLA(JA)=(IDE-1)-2
1947 KHHA(JA)=(IDE-1)-1-MOD(J+1,2)
1948 164 grid%emt(JA)=EMTJ(J)
1950 ! --------------SPREADING OF UPSTREAM VELOCITY-POINT ADVECTION FACTOR----
1952 JA=0
1953 DO 171 J=3,5
1954 JA=JA+1
1955 KVLA(JA)=2
1956 KVHA(JA)=(IDE-1)-1-MOD(J,2)
1957 171 grid%em(JA)=EMJ(J)
1958 DO 172 J=(JDE-1)-4,(JDE-1)-2
1959 JA=JA+1
1960 KVLA(JA)=2
1961 KVHA(JA)=(IDE-1)-1-MOD(J,2)
1962 172 grid%em(JA)=EMJ(J)
1963 DO 173 J=6,(JDE-1)-5
1964 JA=JA+1
1965 KVLA(JA)=2
1966 KVHA(JA)=2+MOD(J+1,2)
1967 173 grid%em(JA)=EMJ(J)
1968 DO 174 J=6,(JDE-1)-5
1969 JA=JA+1
1970 KVLA(JA)=(IDE-1)-2
1971 KVHA(JA)=(IDE-1)-1-MOD(J,2)
1972 174 grid%em(JA)=EMJ(J)
1974 696 continue
1975 ENDIF ! wrf_dm_on_monitor/serial job
1977 call NMM_SH2O(IMS,IME,JMS,JME,ITS,NNXP,JTS,NNYP,grid%num_soil_layers,grid%isltyp, &
1978 grid%sm,grid%sice,grid%stc,grid%smc,grid%sh2o)
1980 !! must be a better place to put this, but will eliminate "phantom"
1981 !! wind points here (no wind point on eastern boundary of odd numbered rows)
1983 IF ( abs(IDE-1-ITE) .eq. 1 ) THEN ! along eastern boundary
1984 write(message,*) 'zero phantom winds'
1985 CALL wrf_message(message)
1986 DO K=1,KDE-1
1987 DO J=JDS,JDE-1,2
1988 IF (J .ge. JTS .and. J .le. JTE) THEN
1989 grid%u(IDE-1,J,K)=0.
1990 grid%v(IDE-1,J,K)=0.
1991 ENDIF
1992 ENDDO
1993 ENDDO
1994 ENDIF
1996 969 continue
1998 DO j = jms, jme
1999 DO i = ims, ime
2000 fisx=max(grid%fis(i,j),0.)
2001 grid%z0(I,J) =grid%sm(I,J)*Z0SEA+(1.-grid%sm(I,J))* &
2002 & (0.*Z0MAX+FISx *FCM+Z0LAND)
2003 ENDDO
2004 ENDDO
2006 write(message,*) 'grid%z0 over memory, leaving module_initialize_real'
2007 CALL wrf_message(message)
2008 DO J=JME,JMS,-((JME-JMS)/20+1)
2009 write(message,635) (grid%z0(I,J),I=IMS,IME,(IME-IMS)/14+1)
2010 CALL wrf_message(message)
2011 ENDDO
2014 endif ! on first_time check
2016 write(message,*) 'leaving init_domain_nmm'
2017 CALL wrf_message( TRIM(message) )
2018 !
2019 write(message,*)'STUFF MOVED TO REGISTRY:',grid%IDTAD, &
2020 & grid%NSOIL,grid%NRADL,grid%NRADS,grid%NPHS,grid%NCNVC,grid%sigma
2021 CALL wrf_message( TRIM(message) )
2022 #ifdef HWRF
2023 !=========================================================================================
2024 ! gopal's doing for ocean coupling. Produce a high resolution grid for the entire domain
2025 !=========================================================================================
2027 if(internal_time_loop.eq.1) then !Kwon's doing
2029 NDLMD=grid%dlmd/3.
2030 NDPHD=grid%dphd/3.
2031 NIDE=3*(IDE-1)-2
2032 NJDE=3*(JDE-1)-2
2033 ILOC=1
2034 JLOC=1
2035 NWBD= WBD ! + (ILOC -1)*2.*grid%dlmd + MOD(JLOC+1,2)*grid%dlmd
2036 NSBD= SBD ! + (JLOC -1)*grid%dphd
2038 ALLOCATE (NHLAT(NIDE,NJDE))
2039 ALLOCATE (NHLON(NIDE,NJDE))
2040 ALLOCATE (NVLAT(NIDE,NJDE))
2041 ALLOCATE (NVLON(NIDE,NJDE))
2042 ALLOCATE (HRES_SM(NIDE,NJDE))
2044 #if defined(DM_PARALLEL)
2045 if(wrf_dm_on_monitor()) then
2046 ! Only the monitor process does the actual work (kinda
2047 ! stupid; should be parallelized, but it's better than
2048 ! writing garbage like it did before with >1 process)
2050 ! Get high-res lat & lon:
2051 CALL EARTH_LATLON_hwrf ( NHLAT,NHLON,NVLAT,NVLON, & ! rotated lat,lon at H and V points
2052 NDLMD,NDPHD,NWBD,NSBD, &
2053 tph0d,tlm0d, &
2054 1,NIDE,1,NJDE,1,1, &
2055 1,NIDE,1,NJDE,1,1, &
2056 1,NIDE,1,NJDE,1,1 )
2058 ! Interpolate landmask to high-res grid:
2059 CALL G2T2H_hwrf ( SM_G,HRES_SM, & ! output grid index and weights
2060 NHLAT,NHLON, & ! target (hres) input lat lon in degrees
2061 grid%DLMD,grid%DPHD,WBD,SBD, & ! parent res, west and south boundaries
2062 tph0d,tlm0d, & ! parent central lat,lon, all in degrees
2063 IDE,JDE,IDS,IDE,JDS,JDE, & ! parent imax and jmax, ime,jme
2064 1,NIDE,1,NJDE,1,1, &
2065 1,NIDE,1,NJDE,1,1, &
2066 1,NIDE,1,NJDE,1,1 )
2068 ! We're done with the low-res sm grid now:
2069 deallocate(SM_G)
2071 ! Write out high-res grid for coupler:
2072 WRITE(65)NHLAT(1:NIDE,1:NJDE)
2073 WRITE(65)NHLON(1:NIDE,1:NJDE)
2074 WRITE(65)NVLAT(1:NIDE,1:NJDE)
2075 WRITE(65)NVLON(1:NIDE,1:NJDE)
2076 WRITE(65)HRES_SM(1:NIDE,1:NJDE)
2077 endif
2078 #else
2079 ! This code is the same as above, but for the non-mpi version:
2080 CALL EARTH_LATLON_hwrf ( NHLAT,NHLON,NVLAT,NVLON, & ! rotated lat,lon at H and V points
2081 NDLMD,NDPHD,NWBD,NSBD, &
2082 tph0d,tlm0d, &
2083 1,NIDE,1,NJDE,1,1, &
2084 1,NIDE,1,NJDE,1,1, &
2085 1,NIDE,1,NJDE,1,1 )
2086 CALL G2T2H_hwrf ( grid%SM,HRES_SM, & ! output grid index and weights
2087 NHLAT,NHLON, & ! target (hres) input lat lon in degrees
2088 grid%DLMD,grid%DPHD,WBD,SBD, & ! parent res, west and south boundaries
2089 tph0d,tlm0d, & ! parent central lat,lon, all in degrees
2090 IDE,JDE,IMS,IME,JMS,JME, & ! parent imax and jmax, ime,jme
2091 1,NIDE,1,NJDE,1,1, &
2092 1,NIDE,1,NJDE,1,1, &
2093 1,NIDE,1,NJDE,1,1 )
2095 WRITE(65)NHLAT(1:NIDE,1:NJDE)
2096 WRITE(65)NHLON(1:NIDE,1:NJDE)
2097 WRITE(65)NVLAT(1:NIDE,1:NJDE)
2098 WRITE(65)NVLON(1:NIDE,1:NJDE)
2099 WRITE(65)HRES_SM(1:NIDE,1:NJDE)
2100 #endif
2102 DEALLOCATE (NHLAT)
2103 DEALLOCATE (NHLON)
2104 DEALLOCATE (NVLAT)
2105 DEALLOCATE (NVLON)
2106 DEALLOCATE (HRES_SM)
2108 endif !Kwon's doing
2110 !==================================================================================
2111 ! end gopal's doing for ocean coupling.
2112 !==================================================================================
2113 #endif
2115 !#define COPY_OUT
2116 !#include <scalar_derefs.inc>
2117 RETURN
2119 END SUBROUTINE init_domain_nmm
2121 !------------------------------------------------------
2123 SUBROUTINE define_nmm_vertical_coord ( LM, PTSGM, pt, pdtop,HYBLEVS, &
2124 SG1,DSG1,SGML1, &
2125 SG2,DSG2,SGML2,dfl, dfrlg )
2127 IMPLICIT NONE
2129 ! USE module_model_constants
2131 !!! certain physical parameters here probably don't need to be defined, as defined
2132 !!! elsewhere within WRF. Done for initial testing purposes.
2134 INTEGER :: LM, LPT2, L
2135 REAL :: PTSGM, pt, PL, PT2, pdtop
2136 REAL :: RGOG, PSIG,PHYB,PHYBM
2137 REAL, PARAMETER :: Rd = 287.04 ! J deg{-1} kg{-1}
2138 REAL, PARAMETER :: CP=1004.6,GAMMA=.0065,PRF0=101325.,T0=288.
2139 REAL, PARAMETER :: g=9.81
2141 REAL, DIMENSION(LM) :: DSG,DSG1,DSG2
2142 REAL, DIMENSION(LM) :: SGML1,SGML2
2143 REAL, DIMENSION(LM+1) :: SG1,SG2,HYBLEVS,dfl,dfrlg
2145 CHARACTER(LEN=255) :: message
2147 LPT2=LM+1
2149 write(message,*) 'pt= ', pt
2150 CALL wrf_message(message)
2152 DO L=LM+1,1,-1
2153 pl=HYBLEVS(L)*(101325.-pt)+pt
2154 if(pl.lt.ptSGm) LPT2=l
2155 ENDDO
2157 IF(LPT2.lt.LM+1) THEN
2158 pt2=HYBLEVS(LPT2)*(101325.-pt)+pt
2159 ELSE
2160 pt2=pt
2161 ENDIF
2163 write(message,*) '*** Sigma system starts at ',pt2,' Pa, from level ',LPT2
2164 CALL wrf_message(message)
2166 pdtop=pt2-pt
2168 write(message,*) 'allocating DSG,DSG1,DSG2 as ', LM
2169 CALL wrf_debug(10,message)
2171 DSG=-99.
2173 DO L=1,LM
2174 DSG(L)=HYBLEVS(L)- HYBLEVS(L+1)
2175 ENDDO
2177 DSG1=0.
2178 DSG2=0.
2180 DO L=LM,1,-1
2182 IF(L.ge.LPT2) then
2183 DSG1(L)=DSG(L)
2184 ELSE
2185 DSG2(L)=DSG(L)
2186 ENDIF
2188 ENDDO
2190 SGML1=-99.
2191 SGML2=-99.
2193 IF(LPT2.le.LM+1) THEN
2195 DO L=LM+1,LPT2,-1
2196 SG2(L)=0.
2197 ENDDO
2199 DO L=LPT2,2,-1
2200 SG2(L-1)=SG2(L)+DSG2(L-1)
2201 ENDDO
2203 DO L=LPT2-1,1,-1
2204 SG2(L)=SG2(L)/SG2(1)
2205 ENDDO
2206 SG2(1)=1.
2208 DO L=LPT2-1,1,-1
2209 DSG2(L)=SG2(L)-SG2(L+1)
2210 SGML2(l)=(SG2(l)+SG2(l+1))*0.5
2211 ENDDO
2213 ENDIF
2215 DO L=LM,LPT2,-1
2216 DSG2(L)=0.
2217 SGML2(L)=0.
2218 ENDDO
2220 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2222 SG1(LM+1)=0.
2224 DO L=LM+1,LPT2,-1
2225 SG1(L-1)=SG1(L)+DSG1(L-1)
2226 ENDDO
2228 DO L=LM,LPT2,-1
2229 SG1(L)=SG1(L)/SG1(LPT2-1)
2230 ENDDO
2232 SG1(LPT2-1)=1.
2234 do l=LPT2-2,1,-1
2235 SG1(l)=1.
2236 enddo
2239 DO L=LM,LPT2,-1
2240 DSG1(L)=SG1(L)-SG1(L+1)
2241 SGML1(L)=(SG1(L)+SG1(L+1))*0.5
2242 ENDDO
2244 DO L=LPT2-1,1,-1
2245 DSG1(L)=0.
2246 SGML1(L)=1.
2247 ENDDO
2249 1000 format('l,hyblevs,psig,SG1,SG2,phyb,phybm')
2250 1100 format(' ',i4,f7.4,f10.2,2f7.4,2f10.2)
2252 write(message,1000)
2253 CALL wrf_debug(100,message)
2255 do l=1,LM+1
2256 psig=HYBLEVS(L)*(101325.-pt)+pt
2257 phyb=SG1(l)*pdtop+SG2(l)*(101325.-pdtop-pt)+pt
2258 if(l.lt.LM+1) then
2259 phybm=SGML1(l)*pdtop+SGML2(l)*(101325.-pdtop-pt)+pt
2260 else
2261 phybm=-99.
2262 endif
2264 write(message,1100) l,HYBLEVS(L),psig &
2265 ,SG1(l),SG2(l),phyb,phybm
2266 CALL wrf_debug(100,message)
2267 enddo
2270 632 format(f9.6)
2272 write(message,*) 'SG1'
2273 CALL wrf_debug(100,message)
2274 do L=LM+1,1,-1
2275 write(message,632) SG1(L)
2276 CALL wrf_debug(100,message)
2277 enddo
2279 write(message,*) 'SG2'
2280 CALL wrf_debug(100,message)
2281 do L=LM+1,1,-1
2282 write(message,632) SG2(L)
2283 CALL wrf_debug(100,message)
2284 enddo
2286 write(message,*) 'DSG1'
2287 CALL wrf_debug(100,message)
2288 do L=LM,1,-1
2289 write(message,632) DSG1(L)
2290 CALL wrf_debug(100,message)
2291 enddo
2293 write(message,*) 'DSG2'
2294 CALL wrf_debug(100,message)
2295 do L=LM,1,-1
2296 write(message,632) DSG2(L)
2297 CALL wrf_debug(100,message)
2298 enddo
2300 write(message,*) 'SGML1'
2301 CALL wrf_debug(100,message)
2302 do L=LM,1,-1
2303 write(message,632) SGML1(L)
2304 CALL wrf_debug(100,message)
2305 enddo
2307 write(message,*) 'SGML2'
2308 CALL wrf_debug(100,message)
2309 do L=LM,1,-1
2310 write(message,632) SGML2(L)
2311 CALL wrf_debug(100,message)
2312 enddo
2314 rgog=(rd*gamma)/g
2315 DO L=1,LM+1
2316 dfl(L)=g*T0*(1.-((pt+SG1(L)*pdtop+SG2(L)*(101325.-pt2)) &
2317 /101325.)**rgog)/gamma
2318 dfrlg(L)=dfl(L)/g
2319 write(message,*) 'L, dfl(L): ', L, dfl(L)
2320 CALL wrf_debug(100,message)
2321 ENDDO
2323 END SUBROUTINE define_nmm_vertical_coord
2325 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2326 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2327 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2329 SUBROUTINE compute_nmm_surfacep ( TERRAIN_HGT_T, Z3D_IN, PRESS3D_IN, T3D_IN &
2330 &, psfc_out,generic &
2331 &, IDS,IDE,JDS,JDE,KDS,KDE &
2332 &, IMS,IME,JMS,JME,KMS,KME &
2333 &, ITS,ITE,JTS,JTE,KTS,KTE )
2336 IMPLICIT NONE
2338 real, allocatable:: dum2d(:,:),DUM2DB(:,:)
2340 integer :: IDS,IDE,JDS,JDE,KDS,KDE
2341 integer :: IMS,IME,JMS,JME,KMS,KME
2342 integer :: ITS,ITE,JTS,JTE,KTS,KTE,Ilook,Jlook
2343 integer :: I,J,II,generic,L,KINSERT,K,bot_lev,LL
2344 integer :: IHE(JMS:JME),IHW(JMS:JME)
2346 real :: TERRAIN_HGT_T(IMS:IME,JMS:JME)
2347 real :: Z3D_IN(IMS:IME,JMS:JME,generic)
2348 real :: T3D_IN(IMS:IME,JMS:JME,generic)
2349 real :: PRESS3D_IN(IMS:IME,JMS:JME,generic)
2350 real :: PSFC_IN(IMS:IME,JMS:JME),TOPO_IN(IMS:IME,JMS:JME)
2351 real :: psfc_out(IMS:IME,JMS:JME),rincr(IMS:IME,JMS:JME)
2352 real :: dif1,dif2,dif3,dif4,dlnpdz,BOT_INPUT_HGT,BOT_INPUT_PRESS,dpdz,rhs
2353 real :: zin(generic),pin(generic)
2355 character (len=255) :: message
2357 logical :: DEFINED_PSFC(IMS:IME,JMS:JME), DEFINED_PSFCB(IMS:IME,JMS:JME)
2359 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2361 Ilook=25
2362 Jlook=25
2364 DO j = JMS, JME
2365 IHE(J)=MOD(J+1,2)
2366 IHW(J)=IHE(J)-1
2367 ENDDO
2369 DO J=JMS,JME
2370 DO I=IMS,IME
2371 DEFINED_PSFC(I,J)=.FALSE.
2372 DEFINED_PSFCB(I,J)=.FALSE.
2373 IF (PRESS3D_IN(I,J,1) .ne. 200100.) THEN
2374 PSFC_IN(I,J)=PRESS3D_IN(I,J,1)
2375 TOPO_IN(I,J)=Z3D_IN(I,J,1)
2376 ELSE
2377 PSFC_IN(I,J)=PRESS3D_IN(I,J,2)
2378 TOPO_IN(I,J)=Z3D_IN(I,J,2)
2379 ENDIF
2380 ENDDO
2381 ENDDO
2383 ! input surface pressure smoothing over the ocean - still needed for NAM?
2385 II_loop: do II=1,8
2387 CYCLE II_loop
2389 do J=JTS+1,min(JTE,JDE-1)-1
2390 do I=ITS+1,min(ITE,IDE-1)-1
2391 rincr(I,J)=0.
2393 if (PSFC_IN(I,J) .gt. 100000. .and. &
2394 PSFC_IN(I+IHE(J),J+1) .gt. 100000. .and. &
2395 PSFC_IN(I+IHE(J),J-1) .gt. 100000. .and. &
2396 PSFC_IN(I+IHW(J),J+1) .gt. 100000. .and. &
2397 PSFC_IN(I+IHW(J),J-1) .gt. 100000. ) then
2399 dif1=abs(PSFC_IN(I,J)-PSFC_IN(I+IHE(J),J+1))
2400 dif2=abs(PSFC_IN(I,J)-PSFC_IN(I+IHE(J),J-1))
2401 dif3=abs(PSFC_IN(I,J)-PSFC_IN(I+IHW(J),J+1))
2402 dif4=abs(PSFC_IN(I,J)-PSFC_IN(I+IHW(J),J-1))
2404 if (max(dif1,dif2,dif3,dif4) .lt. 200. .and. TOPO_IN(I,J).le. 0.5 .and. &
2405 TOPO_IN(I+IHE(J),J+1) .le. 0.5 .and. &
2406 TOPO_IN(I+IHW(J),J+1) .le. 0.5 .and. &
2407 TOPO_IN(I+IHE(J),J-1) .le. 0.5 .and. &
2408 TOPO_IN(I+IHW(J),J-1) .lt. 0.5) then
2410 rincr(I,J)=0.125*( 4.*PSFC_IN(I,J)+ &
2411 PSFC_IN(I+IHE(J),J+1)+PSFC_IN(I+IHE(J),J-1)+ &
2412 PSFC_IN(I+IHW(J),J+1)+PSFC_IN(I+IHW(J),J-1) ) &
2413 - PSFC_IN(I,J)
2415 ! if (rincr(I,J) .ne. 0 .and. abs(rincr(I,J)) .gt. 20.) then
2416 ! write(message,*) 'II, I,J,rincr: ', II, I,J,rincr(I,J)
2417 ! CALL wrf_message(message)
2418 ! endif
2420 endif
2421 endif
2423 ENDDO
2424 ENDDO
2426 DO J=JTS+1,min(JTE,JDE-1)-1
2427 DO I=ITS+1,min(ITE,IDE-1)-1
2428 PSFC_IN(I,J)=PSFC_IN(I,J) + rincr(I,J)
2429 ENDDO
2430 ENDDO
2432 ! write(message,*) ' -------------------------------------------------- '
2433 ! CALL wrf_message(message)
2435 end do II_loop
2437 ALLOCATE(DUM2D(IMS:IME,JMS:JME))
2439 DO J=JMS,JME
2440 DO I=IMS,IME
2441 DUM2D(I,J)=-9.
2442 END DO
2443 END DO
2445 DO J=JTS,min(JTE,JDE-1)
2446 I_loop: DO I=ITS,min(ITE,IDE-1)
2448 IF (PSFC_IN(I,J) .lt. 0.1) THEN
2449 write(message,*) 'QUITTING BECAUSE I,J, PSFC_IN: ', I,J,PSFC_IN(I,J)
2450 call wrf_error_fatal(message)
2451 ENDIF
2453 BOT_INPUT_PRESS=PSFC_IN(I,J)
2454 BOT_INPUT_HGT=TOPO_IN(I,J)
2456 IF (I .eq. Ilook .AND. J .eq. Jlook) THEN
2458 write(message,*) ' TERRAIN_HGT_T: ', I,J, TERRAIN_HGT_T(I,J)
2459 CALL wrf_message(message)
2460 write(message,*) ' PSFC_IN, TOPO_IN: ', &
2461 I, J, PSFC_IN(I,J),TOPO_IN(I,J)
2462 CALL wrf_message(message)
2464 DO L=1,generic
2465 write(message,*) ' L,PRESS3D_IN, Z3D_IN: ', &
2466 I,J,L, PRESS3D_IN(I,J,L),Z3D_IN(I,J,L)
2467 CALL wrf_debug(10,message)
2468 END DO
2469 ENDIF
2471 DO L=2,generic-1
2473 IF ( PRESS3D_IN(i,j,L) .gt. PSFC_IN(I,J) .AND. &
2474 Z3D_IN(I,J,L) .lt. TERRAIN_HGT_T(I,J) .AND. &
2475 Z3D_IN(I,J,L+1) .gt. TERRAIN_HGT_T(I,J) ) THEN
2477 BOT_INPUT_PRESS=PRESS3D_IN(i,j,L)
2478 BOT_INPUT_HGT=Z3D_IN(I,J,L)
2480 ! IF (I .eq. Ilook .and. J .eq. Jlook) THEN
2481 ! write(message,*) 'BOT_INPUT_PRESS, BOT_INPUT_HGT NOW : ', &
2482 ! Ilook,Jlook, BOT_INPUT_PRESS, BOT_INPUT_HGT
2483 ! CALL wrf_message(message)
2484 ! ENDIF
2486 ENDIF
2487 END DO
2489 !!!!!!!!!!!!!!!!!!!!!! START HYDRO CHECK
2491 IF ( PRESS3D_IN(i,j,1) .ne. 200100. .AND. &
2492 (PSFC_IN(I,J) .gt. PRESS3D_IN(i,j,2) .OR. &
2493 TOPO_IN(I,J) .lt. Z3D_IN(I,J,2)) ) THEN ! extrapolate downward
2495 IF (J .eq. JTS .AND. I .eq. ITS) THEN
2496 write(message,*) 'hydro check - should only be for isobaric input'
2497 CALL wrf_message(message)
2498 ENDIF
2500 IF (Z3D_IN(I,J,2) .ne. TOPO_IN(I,J)) THEN
2501 dpdz=(PRESS3D_IN(i,j,2)-PSFC_IN(I,J))/(Z3D_IN(I,J,2)-TOPO_IN(I,J))
2502 rhs=-9.81*((PRESS3D_IN(i,j,2)+ PSFC_IN(I,J))/2.)/(287.04* T3D_IN(I,J,2))
2504 IF ( abs(PRESS3D_IN(i,j,2)-PSFC_IN(I,J)) .gt. 290.) THEN
2505 IF (dpdz .lt. 1.05*rhs .OR. dpdz .gt. 0.95*rhs) THEN
2506 write(message,*) 'I,J,P(2),Psfc,Z(2),Zsfc: ', &
2507 I,J,PRESS3D_IN(i,j,2),PSFC_IN(I,J),Z3D_IN(I,J,2),TOPO_IN(I,J)
2508 IF (mod(I,5).eq.0 .AND. mod(J,5).eq.0) CALL wrf_debug(50,message)
2509 CYCLE I_loop
2510 ENDIF
2512 ENDIF
2514 ELSE ! z(2) equals TOPO_IN
2516 IF (PRESS3D_IN(i,j,2) .eq. PSFC_IN(I,J)) THEN
2517 ! write(message,*) 'all equal at I,J: ', I,J
2518 ! CALL wrf_message(message)
2519 ELSE
2520 ! write(message,*) 'heights equal, pressures not: ', &
2521 ! PRESS3D_IN(i,j,2), PSFC_IN(I,J)
2522 ! CALL wrf_message(message)
2523 CYCLE I_loop
2524 ENDIF
2526 ENDIF
2528 IF ( abs(PRESS3D_IN(i,j,2)-PSFC_IN(I,J)) .gt. 290.) THEN
2529 IF (PRESS3D_IN(i,j,2) .lt. PSFC_IN(I,J) .and. &
2530 Z3D_IN(I,J,2) .lt. TOPO_IN(I,J)) THEN
2531 ! write(message,*) 'surface data mismatch(a) at I,J: ', I,J
2532 ! CALL wrf_message(message)
2533 CYCLE I_loop
2534 ELSEIF (PRESS3D_IN(i,j,2) .gt. PSFC_IN(I,J) .AND. &
2535 Z3D_IN(I,J,2) .gt. TOPO_IN(I,J)) THEN
2536 ! write(message,*) 'surface data mismatch(b) at I,J: ', I,J
2537 ! CALL wrf_message(message)
2538 CYCLE I_loop
2539 ENDIF
2540 ENDIF
2541 ENDIF
2543 !!!!!!! loop over a few more levels
2545 DO L=3,6
2546 IF ( PRESS3D_IN(i,j,1) .ne. 200100. .AND. &
2547 (((PSFC_IN(I,J)-PRESS3D_IN(i,j,L)) .lt. 400.) .OR. &
2548 TOPO_IN(I,J) .lt. Z3D_IN(I,J,L))) then
2550 IF (Z3D_IN(I,J,L) .ne. TOPO_IN(I,J)) THEN
2551 dpdz=(PRESS3D_IN(i,j,L)-PSFC_IN(I,J))/ &
2552 (Z3D_IN(I,J,L)-TOPO_IN(I,J))
2553 rhs=-9.81*((PRESS3D_IN(i,j,L)+ PSFC_IN(I,J))/2.)/ &
2554 (287.04*T3D_IN(I,J,L))
2555 IF ( abs(PRESS3D_IN(i,j,L)-PSFC_IN(I,J)) .gt. 290.) THEN
2556 IF (dpdz .lt. 1.05*rhs .or. dpdz .gt. 0.95*rhs) THEN
2557 write(message,*) 'I,J,L,Piso,Psfc,Ziso,Zsfc: ', &
2558 I,J,L,PRESS3D_IN(i,j,L),PSFC_IN(I,J),&
2559 Z3D_IN(I,J,L),TOPO_IN(I,J)
2560 IF (mod(I,5).eq.0 .AND. mod(J,5).eq.0) &
2561 CALL wrf_debug(50,message)
2562 CYCLE I_loop
2563 ENDIF
2564 ENDIF
2565 ELSE
2566 IF (PRESS3D_IN(i,j,2) .eq. PSFC_IN(I,J)) THEN
2567 ! write(message,*) 'all equal at I,J: ', I,J
2568 ! CALL wrf_message(message)
2569 ELSE
2570 CYCLE I_loop
2571 ENDIF
2572 ENDIF
2573 ENDIF
2575 IF ( abs(PRESS3D_IN(i,j,L)-PSFC_IN(I,J)) .gt. 290.) THEN
2576 IF (PRESS3D_IN(i,j,L) .lt. PSFC_IN(I,J) .AND. &
2577 Z3D_IN(I,J,L) .lt. TOPO_IN(I,J)) THEN
2578 CYCLE I_loop
2579 ELSEIF (PRESS3D_IN(i,j,L) .gt. PSFC_IN(I,J) .AND. &
2580 Z3D_IN(I,J,L) .gt. TOPO_IN(I,J)) THEN
2581 CYCLE I_loop
2582 ENDIF
2583 ENDIF
2584 END DO
2585 !!!!!!!!!!!!!!!!!!!!!! END HYDRO CHECK
2587 IF (TERRAIN_HGT_T(I,J) .eq. BOT_INPUT_HGT ) THEN
2588 dum2d(I,J)=BOT_INPUT_PRESS
2590 IF (BOT_INPUT_HGT .ne. 0. .and. (BOT_INPUT_HGT-INT(BOT_INPUT_HGT) .ne. 0.) ) THEN
2591 write(message,*) 'with BOT_INPUT_HGT: ', BOT_INPUT_HGT, &
2592 'set dum2d to bot_input_pres: ', I,J,dum2d(I,J)
2593 CALL wrf_message(message)
2594 ENDIF
2596 IF (dum2d(I,J) .lt. 50000. .OR. dum2d(I,J) .gt. 109000.) THEN
2597 write(message,*) 'bad dum2d(a): ', I,J,DUM2D(I,J)
2598 CALL wrf_message(message)
2599 ENDIF
2601 ELSEIF (TERRAIN_HGT_T(I,J) .lt. BOT_INPUT_HGT ) THEN
2603 ! target is below lowest possible input...extrapolate
2605 IF ( BOT_INPUT_PRESS-PRESS3D_IN(I,J,2) .gt. 500. ) THEN
2606 dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,2)) ) / &
2607 (BOT_INPUT_HGT-Z3D_IN(i,j,2))
2608 IF (I .eq. Ilook .and. J .eq. Jlook) THEN
2609 write(message,*) 'I,J,dlnpdz(a): ', I,J,dlnpdz
2610 CALL wrf_message(message)
2611 ENDIF
2613 ELSE
2615 !! thin layer and/or just have lowest level - difference with 3rd level data
2616 IF ( abs(BOT_INPUT_PRESS - PRESS3D_IN(i,j,3)) .gt. 290. ) THEN
2618 dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,3)) ) / &
2619 (BOT_INPUT_HGT-Z3D_IN(i,j,3))
2621 IF (I .eq. Ilook .and. J .eq. Jlook) then
2622 write(message,*) 'p diff: ', BOT_INPUT_PRESS, PRESS3D_IN(i,j,3)
2623 CALL wrf_message(message)
2624 write(message,*) 'z diff: ', BOT_INPUT_HGT, Z3D_IN(i,j,3)
2625 CALL wrf_message(message)
2626 ENDIF
2628 ELSE
2630 !! Loop up to level 7 looking for a sufficiently thick layer
2632 FIND_THICK: DO LL=4,7
2633 IF( abs(BOT_INPUT_PRESS - PRESS3D_IN(i,j,LL)) .gt. 290.) THEN
2634 dlnpdz= (log(BOT_INPUT_PRESS)-log(PRESS3D_IN(i,j,LL)) ) / &
2635 (BOT_INPUT_HGT-Z3D_IN(i,j,LL))
2636 EXIT FIND_THICK
2637 ENDIF
2638 END DO FIND_THICK
2640 ENDIF
2642 ENDIF
2644 dum2d(I,J)= exp(log(BOT_INPUT_PRESS) + dlnpdz * &
2645 (TERRAIN_HGT_T(I,J) - BOT_INPUT_HGT) )
2647 IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
2648 write(message,*) 'bad dum2d(b): ', I,J,DUM2D(I,J)
2649 CALL wrf_message(message)
2650 write(message,*) 'BOT_INPUT_PRESS, dlnpdz, TERRAIN_HGT_T, BOT_INPUT_HGT: ', &
2651 BOT_INPUT_PRESS, dlnpdz, TERRAIN_HGT_T(I,J), BOT_INPUT_HGT
2652 CALL wrf_message(message)
2653 write(message,*) 'Z3D_IN: ', Z3D_IN(I,J,1:10)
2654 CALL wrf_message(message)
2655 write(message,*) 'PRESS3D_IN: ', PRESS3D_IN(I,J,1:10)
2656 CALL wrf_message(message)
2657 ENDIF
2659 ELSE ! target level bounded by input levels
2661 DO L=2,generic-1
2662 IF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,L) .AND. &
2663 TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,L+1) ) THEN
2664 dlnpdz= (log(PRESS3D_IN(i,j,l))-log(PRESS3D_IN(i,j,L+1)) ) / &
2665 (Z3D_IN(i,j,l)-Z3D_IN(i,j,L+1))
2666 dum2d(I,J)= log(PRESS3D_IN(i,j,l)) + &
2667 dlnpdz * (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,L) )
2668 dum2d(i,j)=exp(dum2d(i,j))
2669 IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
2670 write(message,*) 'bad dum2d(c): ', I,J,DUM2D(I,J)
2671 CALL wrf_message(message)
2672 ENDIF
2673 ENDIF
2674 ENDDO
2676 !!! account for situation where BOT_INPUT_HGT < TERRAIN_HGT_T < Z3D_IN(:,2,:)
2677 IF (dum2d(I,J) .eq. -9 .AND. BOT_INPUT_HGT .lt. TERRAIN_HGT_T(I,J) &
2678 .AND. TERRAIN_HGT_T(I,J) .lt. Z3D_IN(I,J,2)) then
2680 IF (mod(I,50) .eq. 0 .AND. mod(J,50) .eq. 0) THEN
2681 write(message,*) 'I,J,BOT_INPUT_HGT, bot_pres, TERRAIN_HGT_T: ', &
2682 I,J,BOT_INPUT_HGT, BOT_INPUT_PRESS, TERRAIN_HGT_T(I,J)
2683 CALL wrf_message(message)
2684 ENDIF
2686 dlnpdz= (log(PSFC_IN(i,j))-log(PRESS3D_IN(i,j,2)) ) / &
2687 (TOPO_IN(i,j)-Z3D_IN(i,j,2))
2688 dum2d(I,J)= log(PSFC_IN(i,j)) + &
2689 dlnpdz * (TERRAIN_HGT_T(I,J) - TOPO_IN(i,j) )
2690 dum2d(i,j)= exp(dum2d(i,j))
2691 IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
2692 write(message,*) 'bad dum2d(d): ', I,J,DUM2D(I,J)
2693 CALL wrf_message(message)
2694 ENDIF
2695 ENDIF
2697 IF (dum2d(I,J) .eq. -9.) THEN
2698 write(message,*) 'must have flukey situation in new ', I,J
2699 CALL wrf_message(message)
2700 write(message,*) 'I,J,BOT_INPUT_HGT, bot_pres, TERRAIN_HGT_T: ', &
2701 I,J,BOT_INPUT_HGT, BOT_INPUT_PRESS, TERRAIN_HGT_T(I,J)
2702 CALL wrf_message(message)
2704 DO L=1,generic-1
2705 IF ( TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,L) ) THEN
2706 ! problematic with HGT_M substitution for "input" surface height?
2707 dum2d(i,j)=PRESS3D_IN(I,J,L)
2708 IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
2709 write(message,*) 'bad dum2d(e): ', I,J,DUM2D(I,J)
2710 CALL wrf_message(message)
2711 ENDIF
2712 ENDIF
2713 ENDDO
2715 IF ( TERRAIN_HGT_T(I,J) .eq. TOPO_IN(I,J)) THEN
2716 dum2d(I,J)=PSFC_IN(I,J)
2717 IF (dum2d(I,J) .lt. 50000. .or. dum2d(I,J) .gt. 108000.) THEN
2718 write(message,*) 'bad dum2d(grid%f): ', I,J,DUM2D(I,J)
2719 CALL wrf_message(message)
2720 ENDIF
2721 write(message,*) 'matched input topo, psfc: ', I,J,TOPO_IN(I,J),PSFC_IN(I,J)
2722 CALL wrf_message(message)
2723 ENDIF
2725 IF (dum2d(I,J) .eq. -9.) THEN
2726 CALL wrf_error_fatal("quitting due to undefined surface pressure")
2727 ENDIF
2728 ENDIF
2730 DEFINED_PSFC(I,J)=.TRUE.
2732 IF (I .eq. Ilook .AND. J .eq. Jlook) THEN
2733 write(message,*) 'newstyle psfc: ', I,J,dum2d(I,J)
2734 CALL wrf_message(message)
2735 ENDIF
2737 ENDIF
2739 ENDDO I_loop
2740 ENDDO
2742 ! write(message,*) 'psfc points (new style)'
2743 ! CALL wrf_message(message)
2745 ! DO J=min(JTE,JDE-1),JTS,-loopinc
2746 ! write(message,633) (dum2d(I,J)/100.,I=ITS,min(ITE,IDE-1),iloopinc)
2747 ! CALL wrf_message(message)
2748 ! END DO
2750 633 format(35(f5.0,1x))
2752 write(message,*) 'PSFC extremes (new style)'
2753 CALL wrf_message(message)
2754 write(message,*) minval(dum2d,MASK=DEFINED_PSFC),maxval(dum2d,MASK=DEFINED_PSFC)
2755 CALL wrf_message(message)
2757 ! IF (minval(dum2d,MASK=DEFINED_PSFC) .lt. 50000. .or. maxval(dum2d,MASK=DEFINED_PSFC) .gt. 108000.) THEN
2759 DO J=JTS,min(JTE,JDE-1)
2760 DO I=ITS,min(ITE,IDE-1)
2762 IF (DEFINED_PSFC(I,J) .AND. dum2d(I,J) .lt. 50000. ) THEN
2763 IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN
2764 WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
2765 CALL wrf_debug(2,message)
2766 ELSE
2767 CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2768 ENDIF
2769 ENDIF
2771 IF (DEFINED_PSFC(I,J) .AND. dum2d(I,J) .gt. 108000. ) THEN
2772 IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN
2773 WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
2774 CALL wrf_debug(2,message)
2775 ELSE
2776 CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2777 ENDIF
2778 ENDIF
2780 END DO
2781 END DO
2785 !! "traditional" isobaric only approach ------------------------------------------------
2787 ALLOCATE (DUM2DB(IMS:IME,JMS:JME))
2788 DO J=JMS,JME
2789 DO I=IMS,IME
2790 DUM2DB(I,J)=-9.
2791 END DO
2792 END DO
2794 DO J=JTS,min(JTE,JDE-1)
2795 DO I=ITS,min(ITE,IDE-1)
2797 IF (TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,2)) THEN ! targ below lowest
2799 IF ( abs(PRESS3D_IN(i,j,2)-PRESS3D_IN(i,j,3)) .gt. 290.) THEN
2800 dlnpdz= (log(PRESS3D_IN(i,j,2))-log(PRESS3D_IN(i,j,3)) ) / &
2801 (Z3D_IN(i,j,2)-Z3D_IN(i,j,3))
2802 ELSE
2803 dlnpdz= (log(PRESS3D_IN(i,j,2))-log(PRESS3D_IN(i,j,4)) ) / &
2804 (Z3D_IN(i,j,2)-Z3D_IN(i,j,4))
2805 ENDIF
2807 DUM2DB(I,J)= exp( log(PRESS3D_IN(i,j,2)) + dlnpdz * &
2808 (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,2)) )
2810 IF (I .eq. Ilook .and. J .eq. Jlook) THEN
2811 write(message,*) 'I,K, trad: dlnpdz, press_in(2), terrain_t, Z3D_IN(2): ', I,J,dlnpdz, &
2812 PRESS3D_IN(i,j,2), TERRAIN_HGT_T(I,J), Z3D_IN(i,j,2)
2813 CALL wrf_message(message)
2814 ENDIF
2816 DEFINED_PSFCB(i,j)=.true.
2818 ELSEIF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,2)) THEN ! target level bounded by input levels
2820 DO L=2,generic-1
2821 IF (TERRAIN_HGT_T(I,J) .gt. Z3D_IN(i,j,L) .AND. &
2822 TERRAIN_HGT_T(I,J) .lt. Z3D_IN(i,j,L+1) ) THEN
2824 dlnpdz= (log(PRESS3D_IN(i,j,l))-log(PRESS3D_IN(i,j,L+1)) ) / &
2825 (Z3D_IN(i,j,l)-Z3D_IN(i,j,L+1))
2827 DUM2DB(I,J)= log(PRESS3D_IN(i,j,l)) + &
2828 dlnpdz * (TERRAIN_HGT_T(I,J) - Z3D_IN(i,j,L) )
2829 DUM2DB(i,j)=exp(DUM2DB(i,j))
2831 DEFINED_PSFCB(i,j)=.true.
2833 IF (DUM2DB(I,J) .lt. 13000.) THEN
2834 write(message,*) 'I,J,L,terrain,Z3d(L),z3d(L+1),p3d(L),p3d(l+1): ', I,J,L, &
2835 TERRAIN_HGT_T(I,J),Z3D_IN(I,J,L),Z3D_IN(I,J,L+1),PRESS3D_IN(I,J,L), &
2836 PRESS3D_IN(I,J,L+1)
2837 CALL wrf_error_fatal(message)
2838 ENDIF
2839 ENDIF
2840 ENDDO
2842 ELSEIF (TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,2)) THEN
2843 DUM2DB(i,j)=PRESS3D_IN(I,J,2)
2844 DEFINED_PSFCB(i,j)=.true.
2845 ENDIF
2847 IF (DUM2DB(I,J) .eq. -9.) THEN
2848 write(message,*) 'must have flukey situation in trad ', I,J
2849 CALL wrf_message(message)
2850 DO L=1,generic-1
2851 IF ( TERRAIN_HGT_T(I,J) .eq. Z3D_IN(i,j,L) ) THEN
2852 DUM2DB(i,j)=PRESS3D_IN(I,J,L)
2853 DEFINED_PSFCB(i,j)=.true.
2854 ENDIF
2855 ENDDO
2856 ENDIF
2858 IF (DUM2DB(I,J) .eq. -9.) THEN
2859 write(message,*) 'HOPELESS PSFC, I QUIT'
2860 CALL wrf_error_fatal(message)
2861 ENDIF
2863 if (I .eq. Ilook .and. J .eq. Jlook) THEN
2864 write(message,*) ' traditional psfc: ', I,J,DUM2DB(I,J)
2865 CALL wrf_message(message)
2866 ENDIF
2868 ENDDO
2869 ENDDO
2871 ! write(message,*) 'psfc points (traditional)'
2872 ! CALL wrf_message(message)
2873 ! DO J=min(JTE,JDE-1),JTS,-loopinc
2874 ! write(message,633) (DUM2DB(I,J)/100.,I=its,min(ite,IDE-1),iloopinc)
2875 ! CALL wrf_message(message)
2876 ! ENDDO
2878 write(message,*) 'PSFC extremes (traditional)'
2879 CALL wrf_message(message)
2880 write(message,*) minval(DUM2DB,MASK=DEFINED_PSFCB),maxval(DUM2DB,MASK=DEFINED_PSFCB)
2881 CALL wrf_message(message)
2883 DO J=JTS,min(JTE,JDE-1)
2884 DO I=ITS,min(ITE,IDE-1)
2886 IF (DEFINED_PSFCB(I,J) .AND. dum2db(I,J) .lt. 50000. ) THEN
2887 IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN
2888 WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
2889 CALL wrf_debug(2,message)
2890 ELSE
2891 CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2892 ENDIF
2893 ENDIF
2895 IF (DEFINED_PSFCB(I,J) .AND. dum2db(I,J) .gt. 108000. ) THEN
2896 IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN
2897 WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
2898 CALL wrf_debug(2,message)
2899 ELSE
2900 CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2901 ENDIF
2902 ENDIF
2904 ! IF (DEFINED_PSFCB(I,J) .AND. ( dum2db(I,J) .lt. 50000. .OR. dum2db(I,J) .gt. 108000. )) THEN
2905 ! IF (TERRAIN_HGT_T(I,J) .gt. -10. .and. TERRAIN_HGT_T(I,J) .lt. 5000.) THEN
2906 ! write(0,*) 'I,J, terrain_hgt_t, dum2db: ', I,J, terrain_hgt_t(I,J),dum2db(I,J)
2907 ! CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2908 ! ELSE
2909 ! WRITE(message,*) 'surface pressure allowed because surface height is extreme value of: ', TERRAIN_HGT_T(I,J)
2910 ! CALL wrf_debug(2,message)
2911 ! ENDIF
2912 ! ENDIF
2914 ENDDO
2915 ENDDO
2917 !!!!! end traditional
2919 DO J=JTS,min(JTE,JDE-1)
2920 DO I=ITS,min(ITE,IDE-1)
2921 IF (DEFINED_PSFCB(I,J) .and. DEFINED_PSFC(I,J)) THEN
2923 IF ( abs(dum2d(I,J)-DUM2DB(I,J)) .gt. 400.) THEN
2924 write(message,*) 'BIG DIFF I,J, dum2d, DUM2DB: ', I,J,dum2d(I,J),DUM2DB(I,J)
2925 CALL wrf_message(message)
2926 ENDIF
2928 !! do we have enough confidence in new style to give it more than 50% weight?
2929 psfc_out(I,J)=0.5*(dum2d(I,J)+DUM2DB(I,J))
2931 ELSEIF (DEFINED_PSFC(I,J)) THEN
2932 psfc_out(I,J)=dum2d(I,J)
2933 ELSEIF (DEFINED_PSFCB(I,J)) THEN
2934 psfc_out(I,J)=DUM2DB(I,J)
2935 ELSE
2936 write(message,*) 'I,J,dum2d,DUM2DB: ', I,J,dum2d(I,J),DUM2DB(I,J)
2937 CALL wrf_message(message)
2938 write(message,*) 'I,J,DEFINED_PSFC(I,J),DEFINED_PSFCB(I,J): ', I,J,DEFINED_PSFC(I,J),DEFINED_PSFCB(I,J)
2939 CALL wrf_message(message)
2940 call wrf_error_fatal("psfc_out completely undefined")
2941 ENDIF
2943 IF (I .eq. Ilook .AND. J .eq. Jlook) THEN
2944 write(message,*) ' combined psfc: ', I,J,psfc_out(I,J)
2945 CALL wrf_message(message)
2946 ENDIF
2948 IF (psfc_out(I,J) .lt. 50000. ) THEN
2949 IF (TERRAIN_HGT_T(I,J) .gt. 4500.) THEN
2950 WRITE(message,*) 'low surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
2951 CALL wrf_debug(2,message)
2952 ELSE
2953 write(message,*) 'possibly bad combo on psfc_out: ', I,J, psfc_out(I,J)
2954 CALL wrf_debug(2,message)
2955 write(message,*) 'DEFINED_PSFC, dum2d: ', DEFINED_PSFC(I,J),dum2d(I,J)
2956 CALL wrf_debug(2,message)
2957 write(message,*) 'DEFINED_PSFCB, DUM2DB: ', DEFINED_PSFCB(I,J),DUM2DB(I,J)
2958 CALL wrf_debug(2,message)
2959 CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2960 ENDIF
2961 ENDIF
2963 IF (psfc_out(I,J) .gt. 108000. ) THEN
2964 IF (TERRAIN_HGT_T(I,J) .lt. -10.) THEN
2965 WRITE(message,*) 'high surface pressure allowed because surface height is: ', TERRAIN_HGT_T(I,J)
2966 CALL wrf_debug(2,message)
2967 ELSE
2968 write(message,*) 'possibly bad combo on psfc_out: ', I,J, psfc_out(I,J)
2969 CALL wrf_debug(2,message)
2970 write(message,*) 'DEFINED_PSFC, dum2d: ', DEFINED_PSFC(I,J),dum2d(I,J)
2971 CALL wrf_debug(2,message)
2972 write(message,*) 'DEFINED_PSFCB, DUM2DB: ', DEFINED_PSFCB(I,J),DUM2DB(I,J)
2973 CALL wrf_debug(2,message)
2974 CALL wrf_error_fatal("quit due to unrealistic surface pressure")
2975 ENDIF
2976 ENDIF
2978 ENDDO
2979 ENDDO
2981 deallocate(dum2d,dum2db)
2983 END SUBROUTINE compute_nmm_surfacep
2985 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2986 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2988 SUBROUTINE compute_3d_pressure(psfc_out,SGML1,SGML2,pdtop,pt &
2989 &, pd,p3d_out &
2990 &, IDS,IDE,JDS,JDE,KDS,KDE &
2991 &, IMS,IME,JMS,JME,KMS,KME &
2992 &, ITS,ITE,JTS,JTE,KTS,KTE )
2995 INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
2996 INTEGER :: IMS,IME,JMS,JME,KMS,KME
2997 INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE
2999 REAL, INTENT(IN) :: psfc_out(IMS:IME,JMS:JME)
3000 REAL, INTENT(IN) :: SGML1(KDE),SGML2(KDE),pdtop,pt
3002 REAL, INTENT(OUT):: p3d_out(IMS:IME,JMS:JME,KDS:KDE-1)
3003 REAL, INTENT(OUT):: pd(IMS:IME,JMS:JME)
3005 CHARACTER (len=255) :: message
3007 ! write(message,*) 'pdtop, pt, psfc_out(1,1): ', pdtop, pt, psfc_out(1,1)
3008 ! CALL wrf_message(message)
3010 DO J=JTS,min(JTE,JDE-1)
3011 DO I=ITS,min(ITE,IDE-1)
3012 pd(I,J)=psfc_out(I,J)-pdtop-pt
3013 ENDDO
3014 ENDDO
3016 DO J=JTS,min(JTE,JDE-1)
3017 DO K=KDS,KDE-1
3018 DO I=ITS,min(ITE,IDE-1)
3019 p3d_out(I,J,K)=pd(I,J)*SGML2(K)+pdtop*SGML1(K)+pt
3021 IF (p3d_out(I,J,K) .ge. psfc_out(I,J) .or. p3d_out(I,J,K) .le. pt) THEN
3022 write(message,*) 'I,K,J,p3d_out: ', I,K,J,p3d_out(I,J,K)
3023 CALL wrf_error_fatal(message)
3024 ENDIF
3026 ENDDO
3027 ENDDO
3028 ENDDO
3030 END SUBROUTINE compute_3d_pressure
3032 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3033 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3034 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3036 SUBROUTINE interp_press2press_lin(press_in,press_out, &
3037 data_in, data_out,generic &
3038 &, extrapolate,ignore_lowest,TFIELD &
3039 &, IDS,IDE,JDS,JDE,KDS,KDE &
3040 &, IMS,IME,JMS,JME,KMS,KME &
3041 &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time )
3043 ! Interpolates data from one set of pressure surfaces to
3044 ! another set of pressures
3046 INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
3047 INTEGER :: IMS,IME,JMS,JME,KMS,KME
3048 INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic
3049 INTEGER :: internal_time
3051 ! REAL, INTENT(IN) :: press_in(IMS:IME,generic,JMS:JME)
3052 REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic)
3053 REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1)
3054 ! REAL, INTENT(IN) :: data_in(IMS:IME,generic,JMS:JME)
3055 REAL, INTENT(IN) :: data_in(IMS:IME,JMS:JME,generic)
3056 REAL, INTENT(OUT) :: data_out(IMS:IME,JMS:JME,KMS:KME)
3057 LOGICAL, INTENT(IN) :: extrapolate, ignore_lowest, TFIELD
3058 LOGICAL :: col_smooth
3060 INTEGER :: i,j
3061 INTEGER :: k,kk
3062 REAL :: desired_press
3063 REAL :: dvaldlnp,dlnp,tadiabat,tiso
3065 REAL, PARAMETER :: ADIAFAC=9.81/1004.
3066 REAL, PARAMETER :: TSTEXTRAPFAC=.0065
3070 DO K=KMS,KME
3071 DO J=JMS,JME
3072 DO I=IMS,IME
3073 DATA_OUT(I,J,K)=-99999.9
3074 ENDDO
3075 ENDDO
3076 ENDDO
3078 IF (ignore_lowest) then
3079 LMIN=2
3080 ELSE
3081 LMIN=1
3082 ENDIF
3084 DO j = JTS, min(JTE,JDE-1)
3085 test_i: DO i = ITS, min(ITE,IDE-1)
3087 IF (internal_time_loop .gt. 1) THEN
3088 IF (J .ne. JDS .and. J .ne. JDE-1 .and. &
3089 I .ne. IDS .and. I .ne. IDE-1 ) THEN
3090 !! not on boundary
3091 CYCLE test_i
3092 ENDIF
3093 ENDIF
3096 col_smooth=.false.
3098 output_loop: DO k = KDS,KDE-1
3100 desired_press = press_out(i,j,k)
3102 if (K .gt. KDS) then
3103 if (TFIELD .and. col_smooth .and. desired_press .le. press_in(i,j,LMIN) &
3104 .and. press_out(i,j,k-1) .ge. press_in(i,j,LMIN)) then
3105 MAX_SMOOTH=K
3106 ! write(message,*) 'I,J, MAX_SMOOTH: ', I,J, MAX_SMOOTH
3107 ! CALL wrf_debug(100,message)
3108 endif
3109 endif
3111 ! keep track of where the extrapolation begins
3113 IF (desired_press .GT. press_in(i,j,LMIN)) THEN
3114 IF (TFIELD .and. K .eq. 1 .and. (desired_press - press_in(i,j,LMIN)) .gt. 3000.) then
3115 col_smooth=.TRUE. ! due to large extrapolation distance
3116 ENDIF
3119 IF ((desired_press - press_in(i,j,LMIN)).LT. 50.) THEN ! 0.5 mb
3120 data_out(i,j,k) = data_in(i,j,LMIN)
3121 ELSE
3122 IF (extrapolate) THEN
3123 ! Extrapolate downward because desired P level is below
3124 ! the lowest level in our input data. Extrapolate using simple
3125 ! 1st derivative of value with respect to ln P for the bottom 2
3126 ! input layers.
3128 ! Add a check to make sure we are not using the gradient of
3129 ! a very thin layer
3131 if (TFIELD) then
3132 tiso=0.5*(data_in(i,j,1)+data_in(i,j,2))
3133 endif
3136 IF ( (press_in(i,j,LMIN)-press_in(i,j,LMIN+1)) .GT. 500.) THEN ! likely isobaric data
3137 dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+1))
3138 dvaldlnp = (data_in(i,j,LMIN) - data_in(i,j,LMIN+1)) / dlnp
3139 ELSE ! assume terrain following
3140 dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+5))
3141 dvaldlnp = (data_in(i,j,LMIN) - data_in(i,j,LMIN+5)) / dlnp
3142 ENDIF
3143 data_out(i,j,k) = data_in(i,j,LMIN) + dvaldlnp * &
3144 ( log(desired_press)-log(press_in(i,j,LMIN)) )
3146 if (TFIELD .and. data_out(i,j,k) .lt. tiso-0.2) then
3148 ! restrict slope to -1K/10 hPa
3149 dvaldlnp=max(dvaldlnp, -1.0/ &
3150 log( press_in(i,j,LMIN) / &
3151 ( press_in(i,j,LMIN)-1000.) ))
3153 data_out(I,J,K)= data_in(i,j,LMIN) + dvaldlnp * &
3154 ( log(desired_press)-log(press_in(i,j,LMIN)) )
3156 elseif (TFIELD .and. data_out(i,j,k) .gt. tiso+0.2) then
3158 ! restrict slope to +0.8K/10 hPa
3159 dvaldlnp=min(dvaldlnp, 0.8/ &
3160 log( press_in(i,j,LMIN) / &
3161 ( press_in(i,j,LMIN)-1000.) ))
3163 data_out(I,J,K)= data_in(i,j,LMIN) + dvaldlnp * &
3164 ( log(desired_press)-log(press_in(i,j,LMIN)) )
3166 endif
3168 ELSE
3169 data_out(i,j,k) = data_in(i,j,LMIN)
3170 ENDIF
3171 ENDIF
3172 ELSE IF (desired_press .LT. press_in(i,j,generic)) THEN
3173 IF ( (press_in(i,j,generic) - desired_press) .LT. 10.) THEN
3174 data_out(i,j,k) = data_in(i,j,generic)
3175 ELSE
3176 IF (extrapolate) THEN
3177 ! Extrapolate upward
3178 IF ((press_in(i,j,generic-1)-press_in(i,j,generic)).GT.50.) THEN
3179 dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-1))
3180 dvaldlnp=(data_in(i,j,generic)-data_in(i,j,generic-1))/dlnp
3181 ELSE
3182 dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-2))
3183 dvaldlnp=(data_in(i,j,generic)-data_in(i,j,generic-2))/dlnp
3184 ENDIF
3185 data_out(i,j,k) = data_in(i,j,generic) + &
3186 dvaldlnp * (log(desired_press)-log(press_in(i,j,generic)))
3187 ELSE
3188 data_out(i,j,k) = data_in(i,j,generic)
3189 ENDIF
3190 ENDIF
3191 ELSE
3192 ! We can trap between two levels and linearly interpolate
3194 input_loop: DO kk = LMIN, generic-1
3195 IF (desired_press .EQ. press_in(i,j,kk) )THEN
3196 data_out(i,j,k) = data_in(i,j,kk)
3197 EXIT input_loop
3198 ELSE IF ( (desired_press .LT. press_in(i,j,kk)) .AND. &
3199 (desired_press .GT. press_in(i,j,kk+1)) ) THEN
3201 ! do trapped in lnp
3203 dlnp = log(press_in(i,j,kk)) - log(press_in(i,j,kk+1))
3204 dvaldlnp = (data_in(i,j,kk)-data_in(i,j,kk+1))/dlnp
3205 data_out(i,j,k) = data_in(i,j,kk+1)+ &
3206 dvaldlnp*(log(desired_press)-log(press_in(i,j,kk+1)))
3208 EXIT input_loop
3209 ENDIF
3211 ENDDO input_loop
3212 ENDIF
3213 ENDDO output_loop
3215 if (col_smooth) then
3216 do K=max(KDS,MAX_SMOOTH-4),MAX_SMOOTH+4
3217 data_out(I,J,K)=0.5*(data_out(I,J,K)+data_out(I,J,K+1))
3218 enddo
3219 endif
3221 ENDDO test_i
3222 ENDDO
3223 END SUBROUTINE interp_press2press_lin
3225 SUBROUTINE wind_adjust(press_in,press_out, &
3226 U_in, V_in,U_out,V_out &
3227 &, generic,depth_replace &
3228 &, IDS,IDE,JDS,JDE,KDS,KDE &
3229 &, IMS,IME,JMS,JME,KMS,KME &
3230 &, ITS,ITE,JTS,JTE,KTS,KTE )
3232 INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
3233 INTEGER :: IMS,IME,JMS,JME,KMS,KME
3234 INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic
3235 INTEGER :: MAXLIN,MAXLOUT
3237 REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic)
3238 REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1)
3239 REAL, INTENT(IN) :: U_in(IMS:IME,JMS:JME,generic)
3240 REAL, INTENT(IN) :: V_in(IMS:IME,JMS:JME,generic)
3241 REAL, INTENT(INOUT) :: U_out(IMS:IME,KMS:KME,JMS:JME)
3242 REAL, INTENT(INOUT) :: V_out(IMS:IME,KMS:KME,JMS:JME)
3243 REAL :: p1d_in(generic)
3244 REAL :: p1d_out(KDS:KDE-1)
3247 DO j = JTS, min(JTE,JDE-1)
3248 DO i = ITS, min(ITE,IDE-1)
3250 ! IF (press_out(I,J,1) .lt. press_in(I,J,2)) then
3251 IF( (press_in(I,J,2)-press_out(I,J,1)) .gt. 200.) then
3253 U_out(I,1,J)=U_in(I,J,2)
3254 V_out(I,1,J)=V_in(I,J,2)
3256 INLOOP: DO L=2,generic
3257 p1d_in(L)=-9999.
3258 IF ( (press_in(I,J,2)-press_in(I,J,L)) .lt. depth_replace) THEN
3259 p1d_in(L)=(press_in(I,J,2)-press_in(I,J,L))
3260 MAXLIN=L
3261 ELSE
3262 p1d_in(L)=(press_in(I,J,2)-press_in(I,J,L))
3263 EXIT INLOOP
3264 ENDIF
3265 END DO INLOOP
3267 OUTLOOP: DO L=KDS,KDE-1
3268 p1d_out(L)=-9999.
3269 IF ( (press_out(I,J,1)-press_out(I,J,L)) .lt. depth_replace) THEN
3270 p1d_out(L)=(press_out(I,J,1)-press_out(I,J,L))
3271 MAXLOUT=L
3272 ELSE
3273 EXIT OUTLOOP
3274 ENDIF
3275 END DO OUTLOOP
3277 DO L=1,MAXLOUT
3278 ptarg=p1d_out(L)
3280 FINDLOOP: DO LL=2,MAXLIN
3282 if (p1d_in(LL) .lt. ptarg .and. p1d_in(LL+1) .gt. ptarg) then
3284 dlnp=log(p1d_in(LL))-log(p1d_in(LL+1))
3285 dudlnp=(U_in(I,J,LL)-U_in(I,J,LL+1))/dlnp
3286 dvdlnp=(V_in(I,J,LL)-V_in(I,J,LL+1))/dlnp
3287 U_out(I,L,J)=U_in(I,J,LL)+dudlnp*(log(ptarg)-log(p1d_in(LL)))
3288 V_out(I,L,J)=V_in(I,J,LL)+dvdlnp*(log(ptarg)-log(p1d_in(LL)))
3290 EXIT FINDLOOP
3291 endif
3293 END DO FINDLOOP
3294 END DO ! MAXLOUT loop
3297 ENDIF
3299 ENDDO
3300 ENDDO
3304 END SUBROUTINE wind_adjust
3305 !--------------------------------------------------------------------
3307 SUBROUTINE interp_press2press_log(press_in,press_out, &
3308 data_in, data_out, generic &
3309 &, extrapolate,ignore_lowest &
3310 &, IDS,IDE,JDS,JDE,KDS,KDE &
3311 &, IMS,IME,JMS,JME,KMS,KME &
3312 &, ITS,ITE,JTS,JTE,KTS,KTE, internal_time )
3314 ! Interpolates ln(data) from one set of pressure surfaces to
3315 ! another set of pressures
3317 INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE
3318 INTEGER :: IMS,IME,JMS,JME,KMS,KME
3319 INTEGER :: ITS,ITE,JTS,JTE,KTS,KTE,generic
3320 INTEGER :: internal_time
3322 ! REAL, INTENT(IN) :: press_in(IMS:IME,generic,JMS:JME)
3323 REAL, INTENT(IN) :: press_in(IMS:IME,JMS:JME,generic)
3324 REAL, INTENT(IN) :: press_out(IMS:IME,JMS:JME,KDS:KDE-1)
3325 ! REAL, INTENT(IN) :: data_in(IMS:IME,generic,JMS:JME)
3326 ! REAL, INTENT(IN) :: data_in(IMS:IME,JMS:JME,generic)
3327 REAL :: data_in(IMS:IME,JMS:JME,generic)
3328 REAL, INTENT(OUT) :: data_out(IMS:IME,JMS:JME,KMS:KME)
3329 LOGICAL, INTENT(IN) :: extrapolate, ignore_lowest
3331 INTEGER :: i,j
3332 INTEGER :: k,kk
3333 REAL :: desired_press
3334 REAL :: dlnvaldlnp,dlnp
3337 DO K=1,generic
3338 DO j = JTS, min(JTE,JDE-1)
3339 DO i = ITS, min(ITE,IDE-1)
3340 DATA_IN(I,J,K)=max(DATA_in(I,J,K),1.e-13)
3341 ENDDO
3342 ENDDO
3343 ENDDO
3345 DO K=KMS,KME
3346 DO J=JMS,JME
3347 DO I=IMS,IME
3348 DATA_OUT(I,J,K)=-99999.9
3349 ENDDO
3350 ENDDO
3351 ENDDO
3353 IF (ignore_lowest) then
3354 LMIN=2
3355 ELSE
3356 LMIN=1
3357 ENDIF
3359 DO j = JTS, min(JTE,JDE-1)
3360 test_i: DO i = ITS, min(ITE,IDE-1)
3362 IF (internal_time .gt. 1) THEN
3363 IF (J .ne. JDS .and. J .ne. JDE-1 .and. &
3364 I .ne. IDS .and. I .ne. IDE-1 ) THEN
3365 !! not on boundary
3366 CYCLE test_i
3367 ENDIF
3368 ENDIF
3371 output_loop: DO k = KDS,KDE-1
3373 desired_press = press_out(i,j,k)
3375 IF (desired_press .GT. press_in(i,j,LMIN)) THEN
3377 IF ((desired_press - press_in(i,j,LMIN)).LT. 10.) THEN ! 0.1 mb
3378 data_out(i,j,k) = data_in(i,j,LMIN)
3379 ELSE
3380 IF (extrapolate) THEN
3381 ! Extrapolate downward because desired P level is below
3382 ! the lowest level in our input data. Extrapolate using simple
3383 ! 1st derivative of value with respect to ln P for the bottom 2
3384 ! input layers.
3386 ! Add a check to make sure we are not using the gradient of
3387 ! a very thin layer
3389 IF ( (press_in(i,j,LMIN)-press_in(i,j,LMIN+1)) .GT. 100.) THEN
3390 dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+1))
3391 dlnvaldlnp = ( log(data_in(i,j,LMIN)) - log(data_in(i,j,LMIN+1)) ) / dlnp
3393 ELSE
3395 dlnp = log(press_in(i,j,LMIN))-log(press_in(i,j,LMIN+2))
3396 dlnvaldlnp = (log(data_in(i,j,LMIN)) - log(data_in(i,j,LMIN+2))) / dlnp
3398 ENDIF
3400 data_out(i,j,k) = exp(log(data_in(i,j,LMIN)) + dlnvaldlnp * &
3401 ( log(desired_press)-log(press_in(i,j,LMIN))))
3402 ELSE
3403 data_out(i,j,k) = data_in(i,j,LMIN)
3404 ENDIF
3405 ENDIF
3406 ELSE IF (desired_press .LT. press_in(i,j,generic)) THEN
3407 IF ( (press_in(i,j,generic) - desired_press) .LT. 10.) THEN
3408 data_out(i,j,k) = data_in(i,j,generic)
3409 ELSE
3410 IF (extrapolate) THEN
3411 ! Extrapolate upward
3412 IF ((press_in(i,j,generic-1)-press_in(i,j,generic)).GT.50.) THEN
3413 dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-1))
3414 dlnvaldlnp=(log(data_in(i,j,generic))-log(data_in(i,j,generic-1)))/dlnp
3415 ELSE
3416 dlnp =log(press_in(i,j,generic))-log(press_in(i,j,generic-2))
3417 dlnvaldlnp=(log(data_in(i,j,generic))-log(data_in(i,j,generic-2)))/dlnp
3418 ENDIF
3419 data_out(i,j,k) = exp(log(data_in(i,j,generic)) + &
3420 dlnvaldlnp * (log(desired_press)-log(press_in(i,j,generic))))
3421 ELSE
3422 data_out(i,j,k) = data_in(i,j,generic)
3423 ENDIF
3424 ENDIF
3425 ELSE
3426 ! We can trap between two levels and linearly interpolate
3428 input_loop: DO kk = LMIN, generic-1
3429 IF (desired_press .EQ. press_in(i,j,kk) )THEN
3430 data_out(i,j,k) = data_in(i,j,kk)
3431 EXIT input_loop
3432 ELSE IF ( (desired_press .LT. press_in(i,j,kk)) .AND. &
3433 (desired_press .GT. press_in(i,j,kk+1)) ) THEN
3435 ! do trapped in lnp
3437 dlnp = log(press_in(i,j,kk)) - log(press_in(i,j,kk+1))
3438 dlnvaldlnp = (log(data_in(i,j,kk))-log(data_in(i,j,kk+1)))/dlnp
3439 data_out(i,j,k) = exp(log(data_in(i,j,kk+1))+ &
3440 dlnvaldlnp*(log(desired_press)-log(press_in(i,j,kk+1))))
3442 EXIT input_loop
3444 ENDIF
3446 ENDDO input_loop
3447 ENDIF
3448 ENDDO output_loop
3449 ENDDO test_i
3450 ENDDO
3451 END SUBROUTINE interp_press2press_log
3453 !-------------------------------------------------------------------
3454 SUBROUTINE rh_to_mxrat (rh, t, p, q , wrt_liquid , &
3455 ids , ide , jds , jde , kds , kde , &
3456 ims , ime , jms , jme , kms , kme , &
3457 its , ite , jts , jte , kts , kte )
3459 IMPLICIT NONE
3461 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3462 ims , ime , jms , jme , kms , kme , &
3463 its , ite , jts , jte , kts , kte
3465 LOGICAL , INTENT(IN) :: wrt_liquid
3467 ! REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(IN) :: p , t
3468 ! REAL , DIMENSION(ims:ime,kms:kme,jms:jme) , INTENT(INOUT) :: rh
3469 REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(IN) :: p , t
3470 REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(INOUT) :: rh
3471 REAL , DIMENSION(ims:ime,jms:jme,kms:kme) , INTENT(OUT) :: q
3473 ! Local vars
3475 INTEGER :: i , j , k
3477 REAL :: ew , q1 , t1
3479 REAL, PARAMETER :: T_REF = 0.0
3480 REAL, PARAMETER :: MW_AIR = 28.966
3481 REAL, PARAMETER :: MW_VAP = 18.0152
3483 REAL, PARAMETER :: A0 = 6.107799961
3484 REAL, PARAMETER :: A1 = 4.436518521e-01
3485 REAL, PARAMETER :: A2 = 1.428945805e-02
3486 REAL, PARAMETER :: A3 = 2.650648471e-04
3487 REAL, PARAMETER :: A4 = 3.031240396e-06
3488 REAL, PARAMETER :: A5 = 2.034080948e-08
3489 REAL, PARAMETER :: A6 = 6.136820929e-11
3491 REAL, PARAMETER :: ES0 = 6.1121
3493 REAL, PARAMETER :: C1 = 9.09718
3494 REAL, PARAMETER :: C2 = 3.56654
3495 REAL, PARAMETER :: C3 = 0.876793
3496 REAL, PARAMETER :: EIS = 6.1071
3497 REAL :: RHS
3498 REAL, PARAMETER :: TF = 273.16
3499 REAL :: TK
3501 REAL :: ES
3502 REAL :: QS
3503 REAL, PARAMETER :: EPS = 0.622
3504 REAL, PARAMETER :: SVP1 = 0.6112
3505 REAL, PARAMETER :: SVP2 = 17.67
3506 REAL, PARAMETER :: SVP3 = 29.65
3507 REAL, PARAMETER :: SVPT0 = 273.15
3509 ! This subroutine computes mixing ratio (q, kg/kg) from basic variables
3510 ! pressure (p, Pa), temperature (t, K) and relative humidity (rh, 1-100%).
3511 ! The reference temperature (t_ref, C) is used to describe the temperature
3512 ! at which the liquid and ice phase change occurs.
3514 DO k = kts , kte
3515 DO j = jts , MIN ( jde-1 , jte )
3516 DO i = its , MIN (ide-1 , ite )
3517 rh(i,j,k) = MIN ( MAX ( rh(i,j,k) , 1. ) , 100. )
3518 END DO
3519 END DO
3520 END DO
3522 IF ( wrt_liquid ) THEN
3523 DO k = kts , kte
3524 DO j = jts , MIN ( jde-1 , jte )
3525 DO i = its , MIN (ide-1 , ite )
3526 es=svp1*10.*EXP(svp2*(t(i,j,k)-svpt0)/(t(i,j,k)-svp3))
3527 qs=eps*es/(p(i,j,k)/100.-es)
3528 q(i,j,k)=MAX(.01*rh(i,j,k)*qs,0.0)
3529 END DO
3530 END DO
3531 END DO
3533 ELSE
3534 DO k = kts , kte
3535 DO j = jts , MIN ( jde-1 , jte )
3536 DO i = its , MIN (ide-1 , ite )
3538 t1 = t(i,j,k) - 273.16
3540 ! Obviously dry.
3542 IF ( t1 .lt. -200. ) THEN
3543 q(i,j,k) = 0
3545 ELSE
3547 ! First compute the ambient vapor pressure of water
3549 IF ( ( t1 .GE. t_ref ) .AND. ( t1 .GE. -47.) ) THEN ! liq phase ESLO
3550 ew = a0 + t1 * (a1 + t1 * (a2 + t1 * (a3 + t1 * (a4 + t1 * (a5 + t1 * a6)))))
3552 ELSE IF ( ( t1 .GE. t_ref ) .AND. ( t1 .LT. -47. ) ) then !liq phas poor ES
3553 ew = es0 * exp(17.67 * t1 / ( t1 + 243.5))
3555 ELSE
3556 tk = t(i,j,k)
3557 rhs = -c1 * (tf / tk - 1.) - c2 * alog10(tf / tk) + &
3558 c3 * (1. - tk / tf) + alog10(eis)
3559 ew = 10. ** rhs
3561 END IF
3563 ! Now sat vap pres obtained compute local vapor pressure
3565 ew = MAX ( ew , 0. ) * rh(i,j,k) * 0.01
3567 ! Now compute the specific humidity using the partial vapor
3568 ! pressures of water vapor (ew) and dry air (p-ew). The
3569 ! constants assume that the pressure is in hPa, so we divide
3570 ! the pressures by 100.
3572 q1 = mw_vap * ew
3573 q1 = q1 / (q1 + mw_air * (p(i,j,k)/100. - ew))
3575 q(i,j,k) = q1 / (1. - q1 )
3577 END IF
3579 END DO
3580 END DO
3581 END DO
3583 END IF
3585 END SUBROUTINE rh_to_mxrat
3587 !--=------------------------------------------------------------------
3589 SUBROUTINE boundary_smooth(h, landmask, grid, nsmth , nrow &
3590 &, IDS,IDE,JDS,JDE,KDS,KDE &
3591 &, IMS,IME,JMS,JME,KMS,KME &
3592 &, ITS,ITE,JTS,JTE,KTS,KTE )
3594 implicit none
3596 TYPE (domain) :: grid
3598 integer :: IDS,IDE,JDS,JDE,KDS,KDE
3599 integer :: IMS,IME,JMS,JME,KMS,KME
3600 integer :: ITS,ITE,JTS,JTE,KTS,KTE
3601 integer:: ihw(JDS:JDE-1),ihe(JDS:JDE-1),nsmth,nrow
3602 real:: h(IMS:IME,JMS:JME),landmask(IMS:IME,JMS:JME)
3603 real :: h_old(IMS:IME,JMS:JME)
3604 real :: hbms(IDS:IDE-1,JDS:JDE-1)
3605 real :: hse(IDS:IDE-1,JDS:JDE-1)
3606 real :: hne(IDS:IDE-1,JDS:JDE-1)
3607 integer :: IPS,IPE,JPS,JPE,KPS,KPE
3608 integer :: ihl, ihh, m2l, ibas,jmelin
3609 integer :: I,J,KS,IOFFSET,JSTART,JEND
3610 character (len=255) :: message
3612 ips=its
3613 ipe=ite
3614 jps=jts
3615 jpe=jte
3616 kps=kts
3617 kpe=kte
3619 do j= JTS,min(JTE,JDE-1)
3620 ihw(J)=-mod(J,2)
3621 ihe(j)=ihw(J)+1
3622 end do
3624 do J=JTS,min(JTE,JDE-1)
3625 do I=ITS,min(ITE,IDE-1)
3626 hbms(I,J)=landmask(I,J)
3627 enddo
3628 enddo
3630 jmelin=(JDE-1)-nrow+1
3631 ibas=nrow/2
3632 m2l=mod(nrow,2)
3634 do j=jts,min(jte,jde-1)
3635 ihl=ibas+mod(j,2)+m2l*mod(J+1,2)
3636 ihh=(IDE-1)-ibas-m2l*mod(J+1,2)
3637 do i=its,min(ite,ide-1)
3638 if (I .ge. ihl .and. I .le. ihh .and. J .ge. nrow .and. J .le. jmelin) then
3639 hbms(I,J)=0.
3640 endif
3641 end do
3642 end do
3644 634 format(30(f2.0,1x))
3646 do KS=1,nsmth
3648 grid%ht_gc=h
3649 #ifdef DM_PARALLEL
3650 # include "HALO_NMM_MG.inc"
3651 #endif
3652 h=grid%ht_gc
3653 h_old=grid%ht_gc
3655 do J=JTS,min(JTE,JDE-1)
3656 do I=ITS, min(ITE,IDE-1)
3657 if (I .ge. (IDS+mod(J,2)) .and. J .gt. JDS .and. J .lt. JDE-1 .and. I .lt. IDE-1) then
3658 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) - &
3659 4. *h_old(i,j) )*hbms(i,j)*0.125+h_old(i,j)
3660 endif
3662 enddo
3663 enddo
3665 ! special treatment for four corners
3667 if (hbms(1,1) .eq. 1 .and. ITS .le. 1 .and. JTS .le. 1) then
3668 h(1,1)=0.75*h(1,1)+0.125*h(1+ihe(1),2)+ &
3669 0.0625*(h(2,1)+h(1,3))
3670 endif
3672 if (hbms(IDE-1,1) .eq. 1 .and. ITE .ge. IDE-2 .and. JTS .le. 1) then
3673 h(IDE-1,1)=0.75*h(IDE-1,1)+0.125*h(IDE-1+ihw(1),2)+ &
3674 0.0625*(h(IDE-1-1,1)+h(IDE-1,3))
3675 endif
3677 if (hbms(1,JDE-1) .eq. 1 .and. ITS .le. 1 .and. JTE .ge. JDE-2) then
3678 h(1,JDE-1)=0.75*h(1,JDE-1)+0.125*h(1+ihe(JDE-1),JDE-1-1)+ &
3679 0.0625*(h(2,JDE-1)+h(1,JDE-1-2))
3680 endif
3682 if (hbms(IDE-1,JDE-1) .eq. 1 .and. ITE .ge. IDE-2 .and. JTE .ge. JDE-2) then
3683 h(IDE-1,JDE-1)=0.75*h(IDE-1,JDE-1)+0.125*h(IDE-1+ihw(JDE-1),JDE-1-1)+ &
3684 0.0625*(h(IDE-1-1,JDE-1)+h(IDE-1,JDE-1-2))
3685 endif
3687 do J=JMS,JME
3688 do I=IMS,IME
3689 grid%ht_gc(I,J)=h(I,J)
3690 enddo
3691 enddo
3692 #ifdef DM_PARALLEL
3693 # include "HALO_NMM_MG.inc"
3694 #endif
3695 do J=JMS,JME
3696 do I=IMS,IME
3697 h(I,J)=grid%ht_gc(I,J)
3698 enddo
3699 enddo
3702 ! S bound
3703 if (JTS .eq. JDS) then
3704 J=JTS
3706 do I=ITS,ITE
3707 if (I .ge. IDS+1 .and. I .le. IDE-2) then
3708 if (hbms(I,J) .eq. 1) then
3709 h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J+1)+h(I+ihe(J),J+1))
3710 endif
3711 endif
3712 enddo
3714 endif
3716 ! N bound
3717 if (JTE .eq. JDE) then
3718 J=JDE-1
3719 write(message,*) 'DOING N BOUND SMOOTHING for J= ', J
3720 CALL wrf_debug(100,message)
3721 do I=ITS,min(ITE,IDE-1)
3722 if (hbms(I,J) .eq. 1 .and. I .ge. IDS+1 .and. I .le. IDE-2) then
3723 h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J-1)+h(I+ihe(J),J-1))
3724 endif
3725 enddo
3726 endif
3728 ! W bound
3729 if (ITS .eq. IDS) then
3730 I=ITS
3731 do J=JTS,min(JTE,JDE-1)
3732 if (hbms(I,J) .eq. 1 .and. J .ge. JDS+2 .and. J .le. JDE-3 .and. mod(J,2) .eq. 1) then
3733 h(I,J)=0.75*h(I,J)+0.125*(h(I+ihe(J),J+1)+h(I+ihe(J),J-1))
3734 endif
3735 enddo
3736 endif
3738 ! E bound
3739 if (ITE .eq. IDE) then
3740 write(message,*) 'DOING E BOUND SMOOTHING for I= ', min(ITE,IDE-1)
3741 CALL wrf_debug(100,message)
3742 I=min(ITE,IDE-1)
3743 do J=JTS,min(JTE,JDE-1)
3744 if (hbms(I,J) .eq. 1 .and. J .ge. JDS+2 .and. J .le. JDE-3 .and. mod(J,2) .eq. 1) then
3745 h(I,J)=0.75*h(I,J)+0.125*(h(I+ihw(J),J+1)+h(I+ihw(J),J-1))
3746 endif
3747 enddo
3748 endif
3750 enddo ! end ks loop
3752 do J=JMS,JME
3753 do I=IMS,IME
3754 grid%ht_gc(I,J)=h(I,J)
3755 enddo
3756 enddo
3757 #ifdef DM_PARALLEL
3758 #include "HALO_NMM_MG.inc"
3759 #endif
3760 do J=JMS,JME
3761 do I=IMS,IME
3762 h(I,J)=grid%ht_gc(I,J)
3763 enddo
3764 enddo
3766 ! extra smoothing along inner boundary
3768 if (JTS .eq. JDS) then
3769 if (ITE .eq. IDE) then
3770 IOFFSET=1
3771 else
3772 IOFFSET=0
3773 endif
3774 ! Southern Boundary
3775 do i=its,min(ITE,IDE-1)-IOFFSET
3776 h(i,2)=0.25*(h(i,1)+h(i+1,1)+ &
3777 h(i,3)+h(i+1,3))
3778 enddo
3779 endif
3782 if (JTE .eq. JDE) then
3783 if (ITE .eq. IDE) then
3784 IOFFSET=1
3785 else
3786 IOFFSET=0
3787 endif
3788 do i=its,min(ITE,IDE-1)-IOFFSET
3789 h(i,(JDE-1)-1)=0.25*(h(i,(JDE-1)-2)+h(i+1,(JDE-1)-2)+ &
3790 h(i,JDE-1)+h(i+1,JDE-1))
3791 enddo
3792 endif
3794 if (JTS .eq. 1) then
3795 JSTART=4
3796 else
3797 JSTART=JTS+mod(JTS,2) ! needs to be even
3798 endif
3800 if (JTE .eq. JDE) then
3801 JEND=(JDE-1)-3
3802 else
3803 JEND=JTE
3804 endif
3806 if (ITS .eq. IDS) then
3808 ! Western Boundary
3809 do j=JSTART,JEND,2
3810 h(1,j)=0.25*(h(1,j-1)+h(2,j-1)+ &
3811 h(1,j+1)+h(2,j+1))
3813 enddo
3814 endif
3817 if (ITE .eq. IDE) then
3818 ! Eastern Boundary
3819 do j=JSTART,JEND,2
3820 h((IDE-1)-1,j)=0.25*(h((IDE-1)-1,j-1)+h((IDE-1),j-1)+ &
3821 h((IDE-1)-1,j+1)+h((IDE-1),j+1))
3822 enddo
3823 endif
3826 END SUBROUTINE boundary_smooth
3828 !--------------------------------------------------------------------
3830 SUBROUTINE monthly_interp_to_date ( field_in , date_str , field_out , &
3831 ids , ide , jds , jde , kds , kde , &
3832 ims , ime , jms , jme , kms , kme , &
3833 its , ite , jts , jte , kts , kte )
3835 ! Linrarly in time interpolate data to a current valid time. The data is
3836 ! assumed to come in "monthly", valid at the 15th of every month.
3838 IMPLICIT NONE
3840 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3841 ims , ime , jms , jme , kms , kme , &
3842 its , ite , jts , jte , kts , kte
3844 CHARACTER (LEN=24) , INTENT(IN) :: date_str
3845 REAL , DIMENSION(ims:ime,jms:jme,12) , INTENT(IN) :: field_in
3846 REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_out
3848 ! Local vars
3850 INTEGER :: i , j , l
3851 INTEGER , DIMENSION(0:13) :: middle
3852 INTEGER :: target_julyr , target_julday , target_date
3853 INTEGER :: julyr , julday , int_month, next_month
3854 REAL :: gmt
3855 CHARACTER (LEN=4) :: yr
3856 CHARACTER (LEN=2) :: mon , day15
3859 WRITE(day15,FMT='(I2.2)') 15
3860 DO l = 1 , 12
3861 WRITE(mon,FMT='(I2.2)') l
3862 CALL get_julgmt ( date_str(1:4)//'-'//mon//'-'//day15//'_'//'00:00:00.0000' , julyr , julday , gmt )
3863 middle(l) = julyr*1000 + julday
3864 END DO
3866 l = 0
3867 middle(l) = middle( 1) - 31
3869 l = 13
3870 middle(l) = middle(12) + 31
3872 CALL get_julgmt ( date_str , target_julyr , target_julday , gmt )
3873 target_date = target_julyr * 1000 + target_julday
3874 find_month : DO l = 0 , 12
3875 IF ( ( middle(l) .LT. target_date ) .AND. ( middle(l+1) .GE. target_date ) ) THEN
3876 DO j = jts , MIN ( jde-1 , jte )
3877 DO i = its , MIN (ide-1 , ite )
3878 int_month = MOD ( l , 12 )
3879 IF ( int_month .EQ. 0 ) int_month = 12
3881 IF (int_month == 12) THEN
3882 next_month=1
3883 ELSE
3884 next_month=int_month+1
3885 ENDIF
3887 field_out(i,j) = ( field_in(i,j,next_month) * ( target_date - middle(l) ) + &
3888 field_in(i,j,int_month ) * ( middle(l+1) - target_date ) ) / &
3889 ( middle(l+1) - middle(l) )
3890 END DO
3891 END DO
3892 EXIT find_month
3893 END IF
3894 END DO find_month
3895 END SUBROUTINE monthly_interp_to_date
3897 !---------------------------------------------------------------------
3898 SUBROUTINE monthly_min_max ( field_in , field_min , field_max , &
3899 ids , ide , jds , jde , kds , kde , &
3900 ims , ime , jms , jme , kms , kme , &
3901 its , ite , jts , jte , kts , kte )
3903 ! Plow through each month, find the max, min values for each i,j.
3905 IMPLICIT NONE
3907 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3908 ims , ime , jms , jme , kms , kme , &
3909 its , ite , jts , jte , kts , kte
3911 REAL , DIMENSION(ims:ime,jms:jme,12) , INTENT(IN) :: field_in
3912 REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_min , field_max
3914 ! Local vars
3916 INTEGER :: i , j , l
3917 REAL :: minner , maxxer
3919 DO j = jts , MIN(jde-1,jte)
3920 DO i = its , MIN(ide-1,ite)
3921 minner = field_in(i,j,1)
3922 maxxer = field_in(i,j,1)
3923 DO l = 2 , 12
3924 IF ( field_in(i,j,l) .LT. minner ) THEN
3925 minner = field_in(i,j,l)
3926 END IF
3927 IF ( field_in(i,j,l) .GT. maxxer ) THEN
3928 maxxer = field_in(i,j,l)
3929 END IF
3930 END DO
3931 field_min(i,j) = minner
3932 field_max(i,j) = maxxer
3933 END DO
3934 END DO
3936 END SUBROUTINE monthly_min_max
3938 !-----------------------------------------------------------------------
3940 SUBROUTINE reverse_vert_coord ( field, start_z, end_z &
3941 &, IDS,IDE,JDS,JDE,KDS,KDE &
3942 &, IMS,IME,JMS,JME,KMS,KME &
3943 &, ITS,ITE,JTS,JTE,KTS,KTE )
3945 IMPLICIT NONE
3947 INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
3948 ims , ime , jms , jme , kms , kme , &
3949 its , ite , jts , jte , kts , kte, &
3950 start_z, end_z
3952 REAL, INTENT(INOUT) :: field(IMS:IME,JMS:JME,end_z)
3953 ! local
3955 INTEGER :: I,J,L
3956 REAL, ALLOCATABLE :: dum3d(:,:,:)
3958 allocate(dum3d(IMS:IME,JMS:JME,end_z))
3960 DO L=start_z,end_z
3961 DO J=jts,min(jte,jde-1)
3962 DO I=its,min(ite,ide-1)
3963 dum3d(I,J,L)=field(I,J,end_z-L+start_z)
3964 END DO
3965 END DO
3966 END DO
3968 DO L=start_z,end_z
3969 DO J=jts,min(jte,jde-1)
3970 DO I=its,min(ite,ide-1)
3971 field(I,J,L)=dum3d(I,J,L)
3972 END DO
3973 END DO
3974 END DO
3976 DEALLOCATE(dum3d)
3978 END SUBROUTINE reverse_vert_coord
3981 !--------------------------------------------------------------------
3983 SUBROUTINE compute_nmm_levels(ninterface, ptop, eta_levels)
3985 USE module_model_constants
3987 IMPLICIT NONE
3989 character(len=132):: message
3990 integer :: ninterface,Lthick,L
3991 real, parameter:: gamma=.0065
3992 real, parameter:: t_stand=288.
3993 real, parameter:: p_stand=101325.
3995 real :: maxdz_compute, ptop
3996 real :: plower,pupper,tlay, sum
3998 real :: eta_levels(ninterface)
3999 real, allocatable:: Z(:)
4000 real, allocatable:: deta_levels_spline(:)
4002 logical:: print_pbl_warn
4004 !----------------------------------------------------
4006 allocate(Z(ninterface))
4007 allocate(deta_levels_spline(ninterface-1))
4009 CALL compute_eta_spline(ninterface-1,deta_levels_spline,ptop)
4011 sum=0.
4012 DO L=1,ninterface-1
4013 sum=sum+deta_levels_spline(L)
4014 ENDDO
4016 eta_levels(1)=1.00
4018 DO L=2,ninterface
4019 eta_levels(L)=eta_levels(L-1)-deta_levels_spline(L-1)
4020 ENDDO
4022 eta_levels(ninterface)=0.00
4024 DO L=2,ninterface-1
4025 eta_levels(L)=0.5*(eta_levels(L))+0.25*(eta_levels(L-1)+eta_levels(L+1))
4026 ENDDO
4028 Z(1)=0.
4029 maxdz_compute=0.
4030 print_pbl_warn=.false.
4032 DO L=2,ninterface
4033 tlay=max( t_stand-gamma*Z(L-1), 216.5)
4034 plower=ptop+(p_stand-ptop)*eta_levels(L-1)
4035 pupper=ptop+(p_stand-ptop)*eta_levels(L)
4036 Z(L)=Z(L-1)+(tlay*r_d/g)*(log(plower)-log(pupper))
4038 if (plower .gt. 85000. .and. pupper .lt. 85000. .and. L .lt. 10) then
4039 print_pbl_warn=.true.
4040 endif
4042 write(message,*) 'L, eta(l), pupper, Z(L): ', L, eta_levels(L),pupper,Z(L)
4043 CALL wrf_debug(100,message)
4045 if (Z(L)-Z(L-1) .gt. maxdz_compute) then
4046 Lthick=L
4047 endif
4049 maxdz_compute=max(maxdz_compute,Z(L)-Z(L-1))
4050 ENDDO
4052 if (print_pbl_warn) then
4053 write(message,*) 'WARNING - PBL MAY BE POORLY RESOLVED WITH NUMBER OF VERTICAL LEVELS'
4054 CALL wrf_message(message)
4055 write(message,*) ' - CONSIDER INCREASING THE VERTICAL RESOLUTION'
4056 CALL wrf_message(message)
4057 endif
4059 write(message,*) 'thickest layer was: ', maxdz_compute , 'meters thick at level: ', Lthick
4060 CALL wrf_message(message)
4062 END SUBROUTINE compute_nmm_levels
4064 !---------------------------
4066 SUBROUTINE compute_eta_spline(LM, dsg, ptop)
4068 IMPLICIT NONE
4070 real:: dsg(LM), ptop, sum, rsum
4071 real, allocatable:: xold(:),dold(:)
4072 real, allocatable:: xnew(:),sgm(:)
4073 real, allocatable:: pps(:),qqs(:),y2s(:)
4074 integer nlev,LM,L,KOLD
4076 IF (LM .ge. 46) THEN
4077 KOLD=9
4078 allocate(xold(KOLD))
4079 allocate(dold(KOLD))
4081 xold(1)=.00
4082 dold(1)=.006
4083 xold(2)=.13
4084 dold(2)=.009
4085 xold(3)=.19
4086 dold(3)=.012
4087 xold(4)=.30
4088 dold(4)=.036
4089 xold(5)=.42
4090 dold(5)=.041
4091 xold(6)=.56
4092 dold(6)=.040
4093 xold(7)=.69
4094 dold(7)=.018
4096 if (ptop .ge. 2000.) then
4097 xold(8)=.90
4098 dold(8)=.012
4099 xold(9)=1.0
4100 dold(9)=.006
4101 else
4102 xold(8)=.90
4103 dold(8)=.008
4104 xold(9)=1.0
4105 dold(9)=.003
4106 endif
4108 ELSE
4110 KOLD=8
4111 allocate(xold(KOLD))
4112 allocate(dold(KOLD))
4114 xold(1)=.00
4115 dold(1)=.006
4116 xold(2)=.18
4117 dold(2)=.015
4118 xold(3)=.32
4119 dold(3)=.035
4120 xold(4)=.50
4121 dold(4)=.040
4122 xold(5)=.68
4123 dold(5)=.030
4124 xold(6)=.75
4125 dold(6)=.017
4126 xold(7)=.85
4127 dold(7)=.012
4129 if (ptop .ge. 2000.) then
4130 xold(8)=1.0
4131 dold(8)=.013
4132 else
4133 xold(8)=1.0
4134 dold(8)=.008
4135 endif
4137 ENDIF
4139 allocate(xnew(lm))
4140 allocate(sgm(lm+1))
4141 allocate(pps(lm))
4142 allocate(qqs(lm))
4143 allocate(y2s(lm))
4145 DO L=1,LM
4146 xnew(l)=float(l-1)/float(lm-1)
4147 ENDDO
4149 y2s=0.
4151 CALL spline(kold,xold,dold,y2s,lm,xnew,dsg,pps,qqs)
4153 sum=0.
4154 DO l=1,lm
4155 sum=sum+dsg(l)
4156 ENDDO
4158 rsum=1./sum
4159 sgm(1)=0.
4161 DO L=1,lm-1
4162 dsg(l)=dsg(l)*rsum
4163 sgm(l+1)=sgm(l)+dsg(l)
4164 ENDDO
4165 sgm(lm+1)=1.
4166 dsg(lm)=sgm(lm+1)-sgm(lm)
4168 END SUBROUTINE compute_eta_spline
4170 ! -------------------------------------------------------------------
4172 subroutine spline(NOLD,XOLD,YOLD,Y2,NNEW,XNEW,YNEW,P,q)
4174 ! ********************************************************************
4175 ! * *
4176 ! * THIS IS A ONE-DIMENSIONAL CUBIC SPLINE FITTING ROUTINE *
4177 ! * PROGRAMED FOR A SMALL SCALAR MACHINE. *
4178 ! * *
4179 ! * PROGRAMER Z. JANJIC *
4180 ! * *
4181 ! * NOLD - NUMBER OF GIVEN VALUES OF THE FUNCTION. MUST BE GE 3. *
4182 ! * XOLD - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
4183 ! * FUNCTION ARE GIVEN. MUST BE IN ASCENDING ORDER. *
4184 ! * YOLD - THE GIVEN VALUES OF THE FUNCTION AT THE POINTS XOLD. *
4185 ! * Y2 - THE SECOND DERIVATIVES AT THE POINTS XOLD. IF NATURAL *
4186 ! * SPLINE IS FITTED Y2(1)=0. AND Y2(NOLD)=0. MUST BE *
4187 ! * SPECIFIED. *
4188 ! * NNEW - NUMBER OF VALUES OF THE FUNCTION TO BE CALCULATED. *
4189 ! * XNEW - LOCATIONS OF THE POINTS AT WHICH THE VALUES OF THE *
4190 ! * FUNCTION ARE CALCULATED. XNEW(K) MUST BE GE XOLD(1) *
4191 ! * AND LE XOLD(NOLD). *
4192 ! * YNEW - THE VALUES OF THE FUNCTION TO BE CALCULATED. *
4193 ! * P, q - AUXILIARY VECTORS OF THE LENGTH NOLD-2. *
4194 ! * *
4195 ! ********************************************************************
4196 !
4197 ! LOG:
4198 !
4199 ! JOVIC - July 2008 - fixed incorrectly dimensioned arrays,
4200 ! PYLE and do loop leading to out of bound array
4201 ! reference
4202 !------
4203 !
4204 ! PYLE - June 2007 - eliminated use of GO TO statements.
4205 !
4206 !-----------------------------------------------------------------------
4207 IMPLICIT NONE
4208 !-----------------------------------------------------------------------
4209 INTEGER,INTENT(IN) :: NNEW,NOLD
4210 REAL,DIMENSION(NOLD),INTENT(IN) :: XOLD,YOLD
4211 REAL,DIMENSION(NNEW),INTENT(IN) :: XNEW
4212 REAL,DIMENSION(NNEW),INTENT(OUT) :: YNEW
4213 REAL,DIMENSION(NOLD+2),INTENT(INOUT) :: P,q,Y2
4214 !
4215 INTEGER :: K,K1,K2,KOLD,NOLDM1, K2_hold, K_hold
4216 REAL :: AK,BK,CK,DEN,DX,DXC,DXL,DXR,DYDXL,DYDXR &
4217 & ,RDX,RTDXC,X,XK,XSQ,Y2K,Y2KP1
4218 !-----------------------------------------------------------------------
4220 NOLDM1=NOLD-1
4222 DXL=XOLD(2)-XOLD(1)
4223 DXR=XOLD(3)-XOLD(2)
4224 DYDXL=(YOLD(2)-YOLD(1))/DXL
4225 DYDXR=(YOLD(3)-YOLD(2))/DXR
4226 RTDXC=0.5/(DXL+DXR)
4228 P(1)= RTDXC*(6.*(DYDXR-DYDXL)-DXL*Y2(1))
4229 q(1)=-RTDXC*DXR
4231 K=3
4232 first_loop: DO K=3,NOLD-1
4233 DXL=DXR
4234 DYDXL=DYDXR
4235 DXR=XOLD(K+1)-XOLD(K)
4236 DYDXR=(YOLD(K+1)-YOLD(K))/DXR
4237 DXC=DXL+DXR
4238 DEN=1./(DXL*q(K-2)+DXC+DXC)
4239 P(K-1)= DEN*(6.*(DYDXR-DYDXL)-DXL*P(K-2))
4240 q(K-1)=-DEN*DXR
4241 END DO first_loop
4243 DO K=NOLDM1,2,-1
4244 Y2(K)=P(K-1)+q(K-1)*Y2(K+1)
4245 K_hold=K
4246 END DO
4248 K=K_hold
4250 !-----------------------------------------------------------------------
4251 second_loop: DO K1=1,NNEW
4252 XK=XNEW(K1)
4253 third_loop: DO K2=2,NOLD
4255 IF(XOLD(K2)>XK)THEN
4256 KOLD=K2-1
4257 K2_hold=K2
4258 exit third_loop
4259 ENDIF
4260 K2_hold=K2
4261 END DO third_loop
4263 IF (XOLD(K2_hold) .le. XK) THEN
4264 YNEW(K1)=YOLD(NOLD)
4265 CYCLE second_loop
4266 ENDIF
4268 IF (K1 .eq. 1 .or. K .ne. KOLD) THEN
4269 K=KOLD
4270 Y2K=Y2(K)
4271 Y2KP1=Y2(K+1)
4272 DX=XOLD(K+1)-XOLD(K)
4273 RDX=1./DX
4274 AK=.1666667*RDX*(Y2KP1-Y2K)
4275 BK=0.5*Y2K
4276 CK=RDX*(YOLD(K+1)-YOLD(K))-.1666667*DX*(Y2KP1+Y2K+Y2K)
4277 ENDIF
4279 X=XK-XOLD(K)
4280 XSQ=X*X
4281 YNEW(K1)=AK*XSQ*X+BK*XSQ+CK*X+YOLD(K)
4283 END DO second_loop
4285 END SUBROUTINE SPLINE
4286 !--------------------------------------------------------------------
4287 SUBROUTINE NMM_SH2O(IMS,IME,JMS,JME,ISTART,IM,JSTART,JM,&
4288 NSOIL,ISLTPK, &
4289 sm,sice,stc,smc,sh2o)
4291 !! INTEGER, PARAMETER:: NSOTYP=9
4292 ! INTEGER, PARAMETER:: NSOTYP=16
4293 INTEGER, PARAMETER:: NSOTYP=19 !!!!!!!!MAYBE???
4295 REAL :: PSIS(NSOTYP),BETA(NSOTYP),SMCMAX(NSOTYP)
4296 REAL :: stc(IMS:IME,NSOIL,JMS:JME), &
4297 smc(IMS:IME,NSOIL,JMS:JME)
4298 REAL :: sh2o(IMS:IME,NSOIL,JMS:JME),sice(IMS:IME,JMS:JME),&
4299 sm(IMS:IME,JMS:JME)
4300 REAL :: HLICE,GRAV,T0,BLIM
4301 INTEGER :: ISLTPK(IMS:IME,JMS:JME)
4302 CHARACTER(LEN=255) :: message
4304 ! Constants used in cold start sh2o initialization
4305 DATA HLICE/3.335E5/,GRAV/9.81/,T0/273.15/
4306 DATA BLIM/5.5/
4307 ! DATA PSIS /0.04,0.62,0.47,0.14,0.10,0.26,0.14,0.36,0.04/
4308 ! DATA BETA /4.26,8.72,11.55,4.74,10.73,8.17,6.77,5.25,4.26/
4309 ! DATA SMCMAX /0.421,0.464,0.468,0.434,0.406, &
4310 ! 0.465,0.404,0.439,0.421/
4313 !!! NOT SURE...PSIS=SATPSI, BETA=BB??
4315 DATA PSIS /0.069, 0.036, 0.141, 0.759, 0.759, 0.355, &
4316 0.135, 0.617, 0.263, 0.098, 0.324, 0.468, &
4317 0.355, 0.000, 0.069, 0.036, 0.468, 0.069, 0.069 /
4319 DATA BETA/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, &
4320 6.66, 8.72, 8.17, 10.73, 10.39, 11.55, &
4321 5.25, 0.00, 2.79, 4.26, 11.55, 2.79, 2.79 /
4323 DATA SMCMAX/0.339, 0.421, 0.434, 0.476, 0.476, 0.439, &
4324 0.404, 0.464, 0.465, 0.406, 0.468, 0.468, &
4325 0.439, 1.000, 0.200, 0.421, 0.468, 0.200, 0.339/
4327 DO K=1,NSOIL
4328 DO J=JSTART,JM
4329 DO I=ISTART,IM
4331 !tst
4332 IF (smc(I,K,J) .gt. SMCMAX(ISLTPK(I,J))) then
4333 if (K .eq. 1) then
4334 write(message,*) 'I,J,reducing smc from ' ,I,J,smc(I,K,J), 'to ', SMCMAX(ISLTPK(I,J))
4335 CALL wrf_debug(100,message)
4336 endif
4337 smc(I,K,J)=SMCMAX(ISLTPK(I,J))
4338 ENDIF
4339 !tst
4341 IF ( (sm(I,J) .lt. 0.5) .and. (sice(I,J) .lt. 0.5) ) THEN
4343 IF (ISLTPK(I,J) .gt. 19) THEN
4344 WRITE(message,*) 'FORCING ISLTPK at : ', I,J
4345 CALL wrf_message(message)
4346 ISLTPK(I,J)=9
4347 ELSEIF (ISLTPK(I,J) .le. 0) then
4348 WRITE(message,*) 'FORCING ISLTPK at : ', I,J
4349 CALL wrf_message(message)
4350 ISLTPK(I,J)=1
4351 ENDIF
4354 ! cold start: determine liquid soil water content (sh2o)
4355 ! sh2o <= smc for t < 273.149K (-0.001C)
4357 IF (stc(I,K,J) .LT. 273.149) THEN
4359 ! first guess following explicit solution for Flerchinger Eqn from Koren
4360 ! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
4362 BX = BETA(ISLTPK(I,J))
4363 IF ( BETA(ISLTPK(I,J)) .GT. BLIM ) BX = BLIM
4365 if ( GRAV*(-PSIS(ISLTPK(I,J))) .eq. 0 ) then
4366 write(message,*) 'TROUBLE'
4367 CALL wrf_message(message)
4368 write(message,*) 'I,J: ', i,J
4369 CALL wrf_message(message)
4370 write(message,*) 'grav, isltpk, psis(isltpk): ', grav,isltpk(I,J),&
4371 psis(isltpk(I,J))
4372 CALL wrf_message(message)
4373 endif
4375 if (BX .eq. 0 .or. stc(I,K,J) .eq. 0) then
4376 write(message,*) 'TROUBLE -- I,J,BX, stc: ', I,J,BX,stc(I,K,J)
4377 CALL wrf_message(message)
4378 endif
4379 FK = (((HLICE/(GRAV*(-PSIS(ISLTPK(I,J)))))* &
4380 ((stc(I,K,J)-T0)/stc(I,K,J)))** &
4381 (-1/BX))*SMCMAX(ISLTPK(I,J))
4382 IF (FK .LT. 0.02) FK = 0.02
4383 sh2o(I,K,J) = MIN ( FK, smc(I,K,J) )
4384 ! ----------------------------------------------------------------------
4385 ! now use iterative solution for liquid soil water content using
4386 ! FUNCTION FRH2O (from the Eta "NOAH" land-surface model) with the
4387 ! initial guess for sh2o from above explicit first guess.
4389 sh2o(I,K,J)=FRH2O_init(stc(I,K,J),smc(I,K,J),sh2o(I,K,J), &
4390 SMCMAX(ISLTPK(I,J)),BETA(ISLTPK(I,J)), &
4391 PSIS(ISLTPK(I,J)))
4393 ELSE ! above freezing
4394 sh2o(I,K,J)=smc(I,K,J)
4395 ENDIF
4398 ELSE ! water point
4399 sh2o(I,K,J)=smc(I,K,J)
4401 ENDIF ! test on land/ice/sea
4402 if (sh2o(I,K,J) .gt. SMCMAX(ISLTPK(I,J))) then
4403 write(message,*) 'sh2o > THAN SMCMAX ', I,J,sh2o(I,K,J),SMCMAX(ISLTPK(I,J)),smc(I,K,J)
4404 CALL wrf_message(message)
4405 endif
4407 ENDDO
4408 ENDDO
4409 ENDDO
4411 END SUBROUTINE NMM_SH2O
4413 !-------------------------------------------------------------------
4415 FUNCTION FRH2O_init(TKELV,smc,sh2o,SMCMAX,B,PSIS)
4417 IMPLICIT NONE
4419 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4420 ! PURPOSE: CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT
4421 ! IF TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION
4422 ! TO SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF
4423 ! KOREN ET AL. (1999, JGR, VOL 104(D16), 19569-19585).
4424 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4425 !
4426 ! New version (JUNE 2001): much faster and more accurate newton iteration
4427 ! achieved by first taking log of eqn cited above -- less than 4
4428 ! (typically 1 or 2) iterations achieves convergence. Also, explicit
4429 ! 1-step solution option for special case of parameter Ck=0, which reduces
4430 ! the original implicit equation to a simpler explicit form, known as the
4431 ! ""Flerchinger Eqn". Improved handling of solution in the limit of
4432 ! freezing point temperature T0.
4433 !
4434 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4435 !
4436 ! INPUT:
4437 !
4438 ! TKELV.........Temperature (Kelvin)
4439 ! smc...........Total soil moisture content (volumetric)
4440 ! sh2o..........Liquid soil moisture content (volumetric)
4441 ! SMCMAX........Saturation soil moisture content (from REDPRM)
4442 ! B.............Soil type "B" parameter (from REDPRM)
4443 ! PSIS..........Saturated soil matric potential (from REDPRM)
4444 !
4445 ! OUTPUT:
4446 ! FRH2O.........supercooled liquid water content.
4447 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4449 REAL B
4450 REAL BLIM
4451 REAL BX
4452 REAL CK
4453 REAL DENOM
4454 REAL DF
4455 REAL DH2O
4456 REAL DICE
4457 REAL DSWL
4458 REAL ERROR
4459 REAL FK
4460 REAL FRH2O_init
4461 REAL GS
4462 REAL HLICE
4463 REAL PSIS
4464 REAL sh2o
4465 REAL smc
4466 REAL SMCMAX
4467 REAL SWL
4468 REAL SWLK
4469 REAL TKELV
4470 REAL T0
4472 INTEGER NLOG
4473 INTEGER KCOUNT
4474 PARAMETER (CK=8.0)
4475 ! PARAMETER (CK=0.0)
4476 PARAMETER (BLIM=5.5)
4477 ! PARAMETER (BLIM=7.0)
4478 PARAMETER (ERROR=0.005)
4480 PARAMETER (HLICE=3.335E5)
4481 PARAMETER (GS = 9.81)
4482 PARAMETER (DICE=920.0)
4483 PARAMETER (DH2O=1000.0)
4484 PARAMETER (T0=273.15)
4486 ! ### LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM) ####
4487 ! ### SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT ####
4488 ! ### IS NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES ####
4489 ! ################################################################
4490 !
4491 BX = B
4492 IF ( B .GT. BLIM ) BX = BLIM
4493 ! ------------------------------------------------------------------
4495 ! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
4496 NLOG=0
4497 KCOUNT=0
4499 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4500 ! C IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), sh2o = smc
4501 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4504 IF (TKELV .GT. (T0 - 1.E-3)) THEN
4506 FRH2O_init=smc
4508 ELSE
4510 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4511 IF (CK .NE. 0.0) THEN
4513 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4514 ! CCCCCCCCC OPTION 1: ITERATED SOLUTION FOR NONZERO CK CCCCCCCCCCC
4515 ! CCCCCCCCCCCC IN KOREN ET AL, JGR, 1999, EQN 17 CCCCCCCCCCCCCCCCC
4517 ! INITIAL GUESS FOR SWL (frozen content)
4518 SWL = smc-sh2o
4519 ! KEEP WITHIN BOUNDS.
4520 IF (SWL .GT. (smc-0.02)) SWL=smc-0.02
4521 IF(SWL .LT. 0.) SWL=0.
4522 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4523 ! C START OF ITERATIONS
4524 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4525 DO WHILE (NLOG .LT. 10 .AND. KCOUNT .EQ. 0)
4526 NLOG = NLOG+1
4527 DF = ALOG(( PSIS*GS/HLICE ) * ( ( 1.+CK*SWL )**2. ) * &
4528 ( SMCMAX/(smc-SWL) )**BX) - ALOG(-(TKELV-T0)/TKELV)
4529 DENOM = 2. * CK / ( 1.+CK*SWL ) + BX / ( smc - SWL )
4530 SWLK = SWL - DF/DENOM
4531 ! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
4532 IF (SWLK .GT. (smc-0.02)) SWLK = smc - 0.02
4533 IF(SWLK .LT. 0.) SWLK = 0.
4534 ! MATHEMATICAL SOLUTION BOUNDS APPLIED.
4535 DSWL=ABS(SWLK-SWL)
4536 SWL=SWLK
4537 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4538 ! CC IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
4539 ! CC WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
4540 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4541 IF ( DSWL .LE. ERROR ) THEN
4542 KCOUNT=KCOUNT+1
4543 END IF
4544 END DO
4545 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4546 ! C END OF ITERATIONS
4547 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4548 ! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
4549 FRH2O_init = smc - SWL
4551 ! CCCCCCCCCCCCCCCCCCCCCCCC END OPTION 1 CCCCCCCCCCCCCCCCCCCCCCCCCCC
4553 ENDIF
4555 IF (KCOUNT .EQ. 0) THEN
4556 ! Print*,'Flerchinger used in NEW version. Iterations=',NLOG
4558 ! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4559 ! CCCCC OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0 CCCCCCCC
4560 ! CCCCCCCCCCCCC IN KOREN ET AL., JGR, 1999, EQN 17 CCCCCCCCCCCCCCC
4562 FK=(((HLICE/(GS*(-PSIS)))*((TKELV-T0)/TKELV))**(-1/BX))*SMCMAX
4563 ! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
4564 IF (FK .LT. 0.02) FK = 0.02
4565 FRH2O_init = MIN ( FK, smc )
4567 ! CCCCCCCCCCCCCCCCCCCCCCCCC END OPTION 2 CCCCCCCCCCCCCCCCCCCCCCCCCC
4569 ENDIF
4571 ENDIF
4573 RETURN
4575 END FUNCTION FRH2O_init
4578 !--------------------------------------------------------------------
4580 SUBROUTINE init_module_initialize
4581 END SUBROUTINE init_module_initialize
4583 !---------------------------------------------------------------------
4585 #ifdef HWRF
4586 ! compute earth lat-lons for before interpolations. This is gopal's doing for ocean coupling
4587 !============================================================================================
4589 SUBROUTINE EARTH_LATLON_hwrf ( HLAT,HLON,VLAT,VLON, & !Earth lat,lon at H and V points
4590 DLMD1,DPHD1,WBD1,SBD1, & !input res,west & south boundaries,
4591 CENTRAL_LAT,CENTRAL_LON, & ! central lat,lon, all in degrees
4592 IDS,IDE,JDS,JDE,KDS,KDE, &
4593 IMS,IME,JMS,JME,KMS,KME, &
4594 ITS,ITE,JTS,JTE,KTS,KTE )
4595 !============================================================================
4596 !
4597 IMPLICIT NONE
4598 INTEGER, INTENT(IN ) :: IDS,IDE,JDS,JDE,KDS,KDE
4599 INTEGER, INTENT(IN ) :: IMS,IME,JMS,JME,KMS,KME
4600 INTEGER, INTENT(IN ) :: ITS,ITE,JTS,JTE,KTS,KTE
4601 REAL, INTENT(IN ) :: DLMD1,DPHD1,WBD1,SBD1
4602 REAL, INTENT(IN ) :: CENTRAL_LAT,CENTRAL_LON
4603 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HLAT,HLON,VLAT,VLON
4605 ! local
4607 INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13)
4608 INTEGER :: I,J
4609 REAL(KIND=KNUM) :: WB,SB,DLM,DPH,TPH0,STPH0,CTPH0
4610 REAL(KIND=KNUM) :: TDLM,TDPH,TLMH,TLMV,TLMH0,TLMV0,TPHH,TPHV,DTR
4611 REAL(KIND=KNUM) :: STPH,CTPH,STPV,CTPV,PI_2
4612 REAL(KIND=KNUM) :: SPHH,CLMH,FACTH,SPHV,CLMV,FACTV
4613 REAL(KIND=KNUM), DIMENSION(IMS:IME,JMS:JME) :: GLATH,GLONH,GLATV,GLONV
4614 !-------------------------------------------------------------------------
4616 !
4617 PI_2 = ACOS(0.)
4618 DTR = PI_2/90.
4619 WB = WBD1 * DTR ! WB: western boundary in radians
4620 SB = SBD1 * DTR ! SB: southern boundary in radians
4621 DLM = DLMD1 * DTR ! DLM: dlamda in radians
4622 DPH = DPHD1 * DTR ! DPH: dphi in radians
4623 TDLM = DLM + DLM ! TDLM: 2.0*dlamda
4624 TDPH = DPH + DPH ! TDPH: 2.0*DPH
4626 ! For earth lat lon only
4628 TPH0 = CENTRAL_LAT*DTR ! TPH0: central lat in radians
4629 STPH0 = SIN(TPH0)
4630 CTPH0 = COS(TPH0)
4632 DO J = JTS,MIN(JTE,JDE) !-1) ! H./ This loop takes care of zig-zag
4633 ! ! \.H starting points along j
4634 TLMH0 = WB - TDLM + MOD(J+1,2) * DLM ! ./ TLMH (rotated lats at H points)
4635 TLMV0 = WB - TDLM + MOD(J,2) * DLM ! H (//ly for V points)
4636 TPHH = SB + (J-1)*DPH ! TPHH (rotated lons at H points) are simple trans.
4637 TPHV = TPHH ! TPHV (rotated lons at V points) are simple trans.
4638 STPH = SIN(TPHH)
4639 CTPH = COS(TPHH)
4640 STPV = SIN(TPHV)
4641 CTPV = COS(TPHV)
4643 ! .H
4644 DO I = ITS,MIN(ITE,IDE) !-1) ! /
4645 TLMH = TLMH0 + I*TDLM ! \.H .U .H
4646 ! !H./ ----><----
4647 SPHH = CTPH0 * STPH + STPH0 * CTPH * COS(TLMH) ! DLM + DLM
4648 GLATH(I,J)=ASIN(SPHH) ! GLATH: Earth Lat in radians
4649 CLMH = CTPH*COS(TLMH)/(COS(GLATH(I,J))*CTPH0) &
4650 - TAN(GLATH(I,J))*TAN(TPH0)
4651 IF(CLMH .GT. 1.) CLMH = 1.0
4652 IF(CLMH .LT. -1.) CLMH = -1.0
4653 FACTH = 1.
4654 IF(TLMH .GT. 0.) FACTH = -1.
4655 GLONH(I,J) = -CENTRAL_LON*DTR + FACTH*ACOS(CLMH)
4657 ENDDO
4659 DO I = ITS,MIN(ITE,IDE) !-1)
4660 TLMV = TLMV0 + I*TDLM
4661 SPHV = CTPH0 * STPV + STPH0 * CTPV * COS(TLMV)
4662 GLATV(I,J) = ASIN(SPHV)
4663 CLMV = CTPV*COS(TLMV)/(COS(GLATV(I,J))*CTPH0) &
4664 - TAN(GLATV(I,J))*TAN(TPH0)
4665 IF(CLMV .GT. 1.) CLMV = 1.
4666 IF(CLMV .LT. -1.) CLMV = -1.
4667 FACTV = 1.
4668 IF(TLMV .GT. 0.) FACTV = -1.
4669 GLONV(I,J) = -CENTRAL_LON*DTR + FACTV*ACOS(CLMV)
4671 ENDDO
4673 ENDDO
4675 ! Conversion to degrees (may not be required, eventually)
4677 DO J = JTS, MIN(JTE,JDE) !-1)
4678 DO I = ITS, MIN(ITE,IDE) !-1)
4679 HLAT(I,J) = GLATH(I,J) / DTR
4680 HLON(I,J)= -GLONH(I,J)/DTR
4681 IF(HLON(I,J) .GT. 180.) HLON(I,J) = HLON(I,J) - 360.
4682 IF(HLON(I,J) .LT. -180.) HLON(I,J) = HLON(I,J) + 360.
4683 !
4684 VLAT(I,J) = GLATV(I,J) / DTR
4685 VLON(I,J) = -GLONV(I,J) / DTR
4686 IF(VLON(I,J) .GT. 180.) VLON(I,J) = VLON(I,J) - 360.
4687 IF(VLON(I,J) .LT. -180.) VLON(I,J) = VLON(I,J) + 360.
4689 ENDDO
4690 ENDDO
4692 END SUBROUTINE EARTH_LATLON_hwrf
4694 SUBROUTINE G2T2H_hwrf( SM,HRES_SM, & ! output grid index and weights
4695 HLAT,HLON, & ! target (nest) input lat lon in degrees
4696 DLMD1,DPHD1,WBD1,SBD1, & ! parent res, west and south boundaries
4697 CENTRAL_LAT,CENTRAL_LON, & ! parent central lat,lon, all in degrees
4698 P_IDE,P_JDE,P_IMS,P_IME,P_JMS,P_JME, & ! parent imax and jmax
4699 IDS,IDE,JDS,JDE,KDS,KDE, & ! target (nest) dIMEnsions
4700 IMS,IME,JMS,JME,KMS,KME, &
4701 ITS,ITE,JTS,JTE,KTS,KTE )
4703 !
4704 !*** Tom Black - Initial Version
4705 !*** Gopal - Revised Version for WRF (includes coincident grid points)
4706 !***
4707 !*** GIVEN PARENT CENTRAL LAT-LONS, RESOLUTION AND WESTERN AND SOUTHERN BOUNDARY,
4708 !*** AND THE NESTED GRID LAT-LONS AT H POINTS, THIS ROUTINE FIRST LOCATES THE
4709 !*** INDICES,IIH,JJH, OF THE PARENT DOMAIN'S H POINTS THAT LIES CLOSEST TO THE
4710 !*** h POINTS OF THE NESTED DOMAIN
4711 !
4712 !============================================================================
4713 !
4714 IMPLICIT NONE
4715 INTEGER, INTENT(IN ) :: IDS,IDE,JDS,JDE,KDS,KDE
4716 INTEGER, INTENT(IN ) :: IMS,IME,JMS,JME,KMS,KME
4717 INTEGER, INTENT(IN ) :: ITS,ITE,JTS,JTE,KTS,KTE
4718 INTEGER, INTENT(IN ) :: P_IDE,P_JDE,P_IMS,P_IME,P_JMS,P_JME
4719 REAL, INTENT(IN ) :: DLMD1,DPHD1,WBD1,SBD1
4720 REAL, INTENT(IN ) :: CENTRAL_LAT,CENTRAL_LON
4721 REAL, DIMENSION(P_IMS:P_IME,P_JMS:P_JME), INTENT(IN) :: SM
4722 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(IN) :: HLAT,HLON
4723 REAL, DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HRES_SM
4725 ! local
4727 INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13)
4728 INTEGER :: IMT,JMT,N2R,MK,K,I,J,DSLP0,DSLOPE,N
4729 INTEGER :: NROW,NCOL,KROWS
4730 REAL(KIND=KNUM) :: X,Y,Z,TLAT,TLON
4731 REAL(KIND=KNUM) :: PI_2,D2R,R2D,GLAT,GLON,DPH,DLM,TPH0,TLM0,WB,SB
4732 REAL(KIND=KNUM) :: ROW,COL,SLP0,TLATHC,TLONHC,DENOM,SLOPE
4733 REAL(KIND=KNUM) :: TLAT1,TLAT2,TLON1,TLON2,DLM1,DLM2,DLM3,DLM4,D1,D2
4734 REAL(KIND=KNUM) :: DLA1,DLA2,DLA3,DLA4,S1,R1,DS1,AN1,AN2,AN3 ! Q
4735 REAL(KIND=KNUM) :: DL1,DL2,DL3,DL4,DL1I,DL2I,DL3I,DL4I,SUMDL,TLONO,TLATO
4736 REAL(KIND=KNUM) :: DTEMP
4737 REAL , DIMENSION(IMS:IME,JMS:JME) :: TLATHX,TLONHX
4738 INTEGER, DIMENSION(IMS:IME,JMS:JME) :: KOUTB
4739 REAL SUM,AMAXVAL
4740 REAL, DIMENSION (4, ims:ime, jms:jme ) :: NBWGT
4741 LOGICAL FLIP
4742 REAL, DIMENSION(IMS:IME,JMS:JME) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4
4743 INTEGER, DIMENSION(IMS:IME,JMS:JME) :: IIH,JJH
4744 !-------------------------------------------------------------------------------
4746 IMT=2*P_IDE-2 ! parent i dIMEnsions
4747 JMT=P_JDE/2 ! parent j dIMEnsions
4748 PI_2=ACOS(0.)
4749 D2R=PI_2/90.
4750 R2D=1./D2R
4751 DPH=DPHD1*D2R
4752 DLM=DLMD1*D2R
4753 TPH0= CENTRAL_LAT*D2R
4754 TLM0=-CENTRAL_LON*D2R ! NOTE THE MINUS HERE
4755 WB=WBD1*D2R ! CONVERT NESTED GRID H POINTS FROM GEODETIC
4756 SB=SBD1*D2R
4757 SLP0=DPHD1/DLMD1
4758 DSLP0=NINT(R2D*ATAN(SLP0))
4759 DS1=SQRT(DPH*DPH+DLM*DLM) ! Q
4760 AN1=ASIN(DLM/DS1)
4761 AN2=ASIN(DPH/DS1)
4764 DO J = JTS,MIN(JTE,JDE) !-1)
4765 DO I = ITS,MIN(ITE,IDE) !-1)
4767 !***
4768 !*** LOCATE TARGET h POINTS (HLAT AND HLON) ON THE PARENT DOMAIN AND
4769 !*** DETERMINE THE INDICES IN TERMS OF THE PARENT DOMAIN. FIRST
4770 !*** CONVERT NESTED GRID h POINTS FROM GEODETIC TO TRANSFORMED
4771 !*** COORDINATE ON THE PARENT GRID
4772 !
4774 GLAT=HLAT(I,J)*D2R
4775 GLON= (360. - HLON(I,J))*D2R
4776 X=COS(TPH0)*COS(GLAT)*COS(GLON-TLM0)+SIN(TPH0)*SIN(GLAT)
4777 Y=-COS(GLAT)*SIN(GLON-TLM0)
4778 Z=COS(TPH0)*SIN(GLAT)-SIN(TPH0)*COS(GLAT)*COS(GLON-TLM0)
4779 TLAT=R2D*ATAN(Z/SQRT(X*X+Y*Y))
4780 TLON=R2D*ATAN(Y/X)
4782 !
4783 ROW=TLAT/DPHD1+JMT ! JMT IS THE CENTRAL ROW OF THE PARENT DOMAIN
4784 COL=TLON/DLMD1+P_IDE-1 ! (P_IDE-1) IS THE CENTRAL COLUMN OF THE PARENT DOMAIN
4785 NROW=INT(ROW + 0.001) ! ROUND-OFF IS AVOIDED WITHOUT USING NINT ON PURPOSE
4786 NCOL=INT(COL + 0.001)
4787 TLAT=TLAT*D2R
4788 TLON=TLON*D2R
4790 ! WRITE(60,*)'============================================================'
4791 ! WRITE(60,*)' ','i=',i,'j=',j
4792 !***
4793 !***
4794 !*** FIRST CONSIDER THE SITUATION WHERE THE POINT h IS AT
4795 !***
4796 !*** V H
4797 !***
4798 !***
4799 !*** h
4800 !*** H V
4801 !***
4802 !*** THEN LOCATE THE NEAREST H POINT ON THE PARENT GRID
4803 !***
4804 IF(MOD(NROW,2).EQ.1.AND.MOD(NCOL,2).EQ.1.OR. &
4805 MOD(NROW,2).EQ.0.AND.MOD(NCOL,2).EQ.0)THEN
4806 TLAT1=(NROW-JMT)*DPH
4807 TLAT2=TLAT1+DPH
4808 TLON1=(NCOL-(P_IDE-1))*DLM
4809 TLON2=TLON1+DLM
4810 DLM1=TLON-TLON1
4811 DLM2=TLON-TLON2
4812 ! D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
4813 ! D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
4814 DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
4815 D1=ACOS(DTEMP)
4816 DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
4817 D2=ACOS(DTEMP)
4818 IF(D1.GT.D2)THEN
4819 NROW=NROW+1 ! FIND THE NEAREST H ROW
4820 NCOL=NCOL+1 ! FIND THE NEAREST H COLUMN
4821 ENDIF
4822 ! WRITE(60,*)'NEAREST PARENT IS:','col=',COL,'row=',ROW,'ncol=',NCOL,'nrow=',NROW
4823 ELSE
4824 !***
4825 !*** NOW CONSIDER THE SITUATION WHERE THE POINT h IS AT
4826 !***
4827 !*** H V
4828 !***
4829 !***
4830 !*** h
4831 !*** V H
4832 !***
4833 !*** THEN LOCATE THE NEAREST H POINT ON THE PARENT GRID
4834 !***
4835 !***
4836 TLAT1=(NROW+1-JMT)*DPH
4837 TLAT2=TLAT1-DPH
4838 TLON1=(NCOL-(P_IDE-1))*DLM
4839 TLON2=TLON1+DLM
4840 DLM1=TLON-TLON1
4841 DLM2=TLON-TLON2
4842 ! D1=ACOS(COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
4843 ! D2=ACOS(COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
4844 DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT1)*COS(DLM1)+SIN(TLAT)*SIN(TLAT1))
4845 D1=ACOS(DTEMP)
4846 DTEMP=MIN(1.0_KNUM,COS(TLAT)*COS(TLAT2)*COS(DLM2)+SIN(TLAT)*SIN(TLAT2))
4847 D2=ACOS(DTEMP)
4848 IF(D1.LT.D2)THEN
4849 NROW=NROW+1 ! FIND THE NEAREST H ROW
4850 ELSE
4851 NCOL=NCOL+1 ! FIND THE NEAREST H COLUMN
4852 ENDIF
4853 ! WRITE(60,*)'NEAREST PARENT IS:','col=',COL,'row=',ROW,'ncol=',NCOL,'nrow=',NROW
4854 ENDIF
4856 KROWS=((NROW-1)/2)*IMT
4857 IF(MOD(NROW,2).EQ.1)THEN
4858 K=KROWS+(NCOL+1)/2
4859 ELSE
4860 K=KROWS+P_IDE-1+NCOL/2
4861 ENDIF
4863 !***
4864 !*** WE NOW KNOW THAT THE INNER GRID POINT IN QUESTION IS
4865 !*** NEAREST TO THE CENTER K AS SEEN BELOW. WE MUST FIND
4866 !*** WHICH OF THE FOUR H-BOXES (OF WHICH THIS H POINT IS
4867 !*** A VERTEX) SURROUNDS THE INNER GRID h POINT IN QUESTION.
4868 !***
4869 !**
4870 !*** H
4871 !***
4872 !***
4873 !***
4874 !*** H V H
4875 !***
4876 !***
4877 !*** h
4878 !*** H V H V H
4879 !***
4880 !***
4881 !***
4882 !*** H V H
4883 !***
4884 !***
4885 !***
4886 !*** H
4887 !***
4888 !***
4889 !*** FIND THE SLOPE OF THE LINE CONNECTING h AND THE CENTER H.
4890 !***
4891 N2R=K/IMT
4892 MK=MOD(K,IMT)
4893 !
4894 IF(MK.EQ.0)THEN
4895 TLATHC=SB+(2*N2R-1)*DPH
4896 ELSE
4897 TLATHC=SB+(2*N2R+(MK-1)/(P_IDE-1))*DPH
4898 ENDIF
4899 !
4900 IF(MK.LE.(P_IDE-1))THEN
4901 TLONHC=WB+2*(MK-1)*DLM
4902 ELSE
4903 TLONHC=WB+(2*(MK-(P_IDE-1))-1)*DLM
4904 ENDIF
4906 !
4907 !*** EXECUTE CAUTION IF YOU NEED TO CHANGE THESE CONDITIONS. SINCE WE ARE
4908 !*** DEALING WITH SLOPES TO GENERATE DIAMOND SHAPE H BOXES, WE NEED TO BE
4909 !*** CAREFUL HERE
4910 !
4912 IF(ABS(TLON-TLONHC) .LE. 1.E-4)TLONHC=TLON
4913 IF(ABS(TLAT-TLATHC) .LE. 1.E-4)TLATHC=TLAT
4914 DENOM=(TLON-TLONHC)
4915 !
4916 !***
4917 !***STORE THE LOCATION OF THE WESTERNMOST VERTEX OF THE H-BOX ON
4918 !***THE OUTER GRID THAT SURROUNDS THE h POINT ON THE INNER GRID.
4919 !***
4920 !*** COINCIDENT CONDITIONS
4922 IF(DENOM.EQ.0.0)THEN
4924 IF(TLATHC.EQ.TLAT)THEN
4925 KOUTB(I,J)=K
4926 IIH(I,J) = NCOL
4927 JJH(I,J) = NROW
4928 TLATHX(I,J)=TLATHC
4929 TLONHX(I,J)=TLONHC
4930 HBWGT1(I,J)=1.0
4931 HBWGT2(I,J)=0.0
4932 HBWGT3(I,J)=0.0
4933 HBWGT4(I,J)=0.0
4934 ! WRITE(60,*)'TRIVIAL SOLUTION'
4935 ELSE ! SAME LONGITUDE BUT DIFFERENT LATS
4936 !
4937 IF(TLATHC .GT. TLAT)THEN ! NESTED POINT SOUTH OF PARENT
4938 KOUTB(I,J)=K-(P_IDE-1)
4939 IIH(I,J) = NCOL-1
4940 JJH(I,J) = NROW-1
4941 TLATHX(I,J)=TLATHC-DPH
4942 TLONHX(I,J)=TLONHC-DLM
4943 ! WRITE(60,*)'VANISHING SLOPE, -ve: TLATHC-DPH, TLONHC-DLM'
4944 ELSE ! NESTED POINT NORTH OF PARENT
4945 KOUTB(I,J)=K+(P_IDE-1)-1
4946 IIH(I,J) = NCOL-1
4947 JJH(I,J) = NROW+1
4948 TLATHX(I,J)=TLATHC+DPH
4949 TLONHX(I,J)=TLONHC-DLM
4950 ! WRITE(60,*)'VANISHING SLOPE, +ve: TLATHC+DPH, TLONHC-DLM'
4951 ENDIF
4952 !***
4953 !***
4954 !*** 4
4955 !***
4956 !*** h
4957 !*** 1 2
4958 !***
4959 !*** 3
4960 !*** DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX
4962 TLATO=TLATHX(I,J)
4963 TLONO=TLONHX(I,J)
4964 DLM1=TLON-TLONO
4965 DLA1=TLAT-TLATO ! Q
4966 ! DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
4967 DL1=SQRT(DLM1*DLM1+DLA1*DLA1) ! Q
4968 !
4969 TLATO=TLATHX(I,J)
4970 TLONO=TLONHX(I,J)+2.*DLM
4971 DLM2=TLON-TLONO
4972 DLA2=TLAT-TLATO ! Q
4973 ! DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
4974 DL2=SQRT(DLM2*DLM2+DLA2*DLA2) ! Q
4975 !
4976 TLATO=TLATHX(I,J)-DPH
4977 TLONO=TLONHX(I,J)+DLM
4978 DLM3=TLON-TLONO
4979 DLA3=TLAT-TLATO ! Q
4980 ! DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
4981 DL3=SQRT(DLM3*DLM3+DLA3*DLA3) ! Q
4983 TLATO=TLATHX(I,J)+DPH
4984 TLONO=TLONHX(I,J)+DLM
4985 DLM4=TLON-TLONO
4986 DLA4=TLAT-TLATO ! Q
4987 ! DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
4988 DL4=SQRT(DLM4*DLM4+DLA4*DLA4) ! Q
4991 ! THE BILINEAR WEIGHTS
4992 !***
4993 !***
4994 AN3=ATAN2(DLA1,DLM1) ! Q
4995 R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
4996 S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
4997 R1=R1/DS1
4998 S1=S1/DS1
4999 DL1I=(1.-R1)*(1.-S1)
5000 DL2I=R1*S1
5001 DL3I=R1*(1.-S1)
5002 DL4I=(1.-R1)*S1
5003 !
5004 HBWGT1(I,J)=DL1I
5005 HBWGT2(I,J)=DL2I
5006 HBWGT3(I,J)=DL3I
5007 HBWGT4(I,J)=DL4I
5008 !
5009 ENDIF
5011 ELSE
5012 !
5013 !*** NON-COINCIDENT POINTS
5014 !
5015 SLOPE=(TLAT-TLATHC)/DENOM
5016 DSLOPE=NINT(R2D*ATAN(SLOPE))
5018 IF(DSLOPE.LE.DSLP0.AND.DSLOPE.GE.-DSLP0)THEN
5019 IF(TLON.GT.TLONHC)THEN
5020 ! IF(TLONHC.GE.-WB-DLM)CALL wrf_error_fatal("1H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
5021 KOUTB(I,J)=K
5022 IIH(I,J) = NCOL
5023 JJH(I,J) = NROW
5024 TLATHX(I,J)=TLATHC
5025 TLONHX(I,J)=TLONHC
5026 ! WRITE(60,*)'HERE WE GO1: TLATHC, TLONHC'
5027 ELSE
5028 ! IF(TLONHC.LE.WB+DLM)CALL wrf_error_fatal("2H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
5029 KOUTB(I,J)=K-1
5030 IIH(I,J) = NCOL-2
5031 JJH(I,J) = NROW
5032 TLATHX(I,J)=TLATHC
5033 TLONHX(I,J)=TLONHC -2.*DLM
5034 ! WRITE(60,*)'HERE WE GO2: TLATHC, TLONHC -2.*DLM'
5035 ENDIF
5037 !
5038 ELSEIF(DSLOPE.GT.DSLP0)THEN
5039 IF(TLON.GT.TLONHC)THEN
5040 ! IF(TLATHC.GE.-SB-DPH)CALL wrf_error_fatal("3H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
5041 KOUTB(I,J)=K+(P_IDE-1)-1
5042 IIH(I,J) = NCOL-1
5043 JJH(I,J) = NROW+1
5044 TLATHX(I,J)=TLATHC+DPH
5045 TLONHX(I,J)=TLONHC-DLM
5046 ! WRITE(60,*)'HERE WE GO3: TLATHC+DPH, TLONHC-DLM'
5047 ELSE
5048 ! IF(TLATHC.LE.SB+DPH)CALL wrf_error_fatal("4H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
5049 KOUTB(I,J)=K-(P_IDE-1)
5050 IIH(I,J) = NCOL-1
5051 JJH(I,J) = NROW-1
5052 TLATHX(I,J)=TLATHC-DPH
5053 TLONHX(I,J)=TLONHC-DLM
5054 ! WRITE(60,*)'HERE WE GO4: TLATHC-DPH, TLONHC-DLM'
5055 ENDIF
5057 !
5058 ELSEIF(DSLOPE.LT.-DSLP0)THEN
5059 IF(TLON.GT.TLONHC)THEN
5060 ! IF(TLATHC.LE.SB+DPH)CALL wrf_error_fatal("5H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
5061 KOUTB(I,J)=K-(P_IDE-1)
5062 IIH(I,J) = NCOL-1
5063 JJH(I,J) = NROW-1
5064 TLATHX(I,J)=TLATHC-DPH
5065 TLONHX(I,J)=TLONHC-DLM
5066 ! WRITE(60,*)'HERE WE GO5: TLATHC-DPH, TLONHC-DLM'
5067 ELSE
5068 ! IF(TLATHC.GE.-SB-DPH)CALL wrf_error_fatal("6H:NESTED DOMAIN TOO CLOSE TO THE BOUNDARY OF PARENT")
5069 KOUTB(I,J)=K+(P_IDE-1)-1
5070 IIH(I,J) = NCOL-1
5071 JJH(I,J) = NROW+1
5072 TLATHX(I,J)=TLATHC+DPH
5073 TLONHX(I,J)=TLONHC-DLM
5074 ! WRITE(60,*)'HERE WE GO6: TLATHC+DPH, TLONHC-DLM'
5075 ENDIF
5076 ENDIF
5078 !
5079 !*** NOW WE WILL MOVE AS FOLLOWS:
5080 !***
5081 !***
5082 !*** 4
5083 !***
5084 !***
5085 !***
5086 !*** h
5087 !*** 1 2
5088 !***
5089 !***
5090 !***
5091 !***
5092 !*** 3
5093 !***
5094 !***
5095 !***
5096 !*** DL 1-4 ARE THE ANGULAR DISTANCES FROM h TO EACH VERTEX
5098 TLATO=TLATHX(I,J)
5099 TLONO=TLONHX(I,J)
5100 DLM1=TLON-TLONO
5101 DLA1=TLAT-TLATO ! Q
5102 ! DL1=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM1)+SIN(TLAT)*SIN(TLATO)) ! Q
5103 DL1=SQRT(DLM1*DLM1+DLA1*DLA1) ! Q
5104 !
5105 TLATO=TLATHX(I,J) ! redundant computations
5106 TLONO=TLONHX(I,J)+2.*DLM
5107 DLM2=TLON-TLONO
5108 DLA2=TLAT-TLATO ! Q
5109 ! DL2=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM2)+SIN(TLAT)*SIN(TLATO)) ! Q
5110 DL2=SQRT(DLM2*DLM2+DLA2*DLA2) ! Q
5111 !
5112 TLATO=TLATHX(I,J)-DPH
5113 TLONO=TLONHX(I,J)+DLM
5114 DLM3=TLON-TLONO
5115 DLA3=TLAT-TLATO ! Q
5116 ! DL3=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM3)+SIN(TLAT)*SIN(TLATO)) ! Q
5117 DL3=SQRT(DLM3*DLM3+DLA3*DLA3) ! Q
5118 !
5119 TLATO=TLATHX(I,J)+DPH
5120 TLONO=TLONHX(I,J)+DLM
5121 DLM4=TLON-TLONO
5122 DLA4=TLAT-TLATO ! Q
5123 ! DL4=ACOS(COS(TLAT)*COS(TLATO)*COS(DLM4)+SIN(TLAT)*SIN(TLATO)) ! Q
5124 DL4=SQRT(DLM4*DLM4+DLA4*DLA4) ! Q
5126 ! THE BILINEAR WEIGHTS
5127 !***
5128 AN3=ATAN2(DLA1,DLM1) ! Q
5129 R1=DL1*SIN(AN2-AN3)/SIN(2.*AN1)
5130 S1=DL1*SIN(2.*PI_2-2*AN1-AN2+AN3)/SIN(2.*AN1)
5131 R1=R1/DS1
5132 S1=S1/DS1
5133 DL1I=(1.-R1)*(1.-S1)
5134 DL2I=R1*S1
5135 DL3I=R1*(1.-S1)
5136 DL4I=(1.-R1)*S1
5137 !
5138 HBWGT1(I,J)=DL1I
5139 HBWGT2(I,J)=DL2I
5140 HBWGT3(I,J)=DL3I
5141 HBWGT4(I,J)=DL4I
5142 !
5143 ENDIF
5144 !
5145 !*** FINALLY STORE IIH IN TERMS OF E-GRID INDEX
5146 !
5147 IIH(I,J)=NINT(0.5*IIH(I,J))
5149 ENDDO
5150 ENDDO
5152 !
5153 !*** EXTENSION TO NEAREST NEIGHBOR
5154 !
5155 DO J = JTS,MIN(JTE,JDE) !-1)
5156 DO I = ITS,MIN(ITE,IDE) !-1)
5157 NBWGT(1,I,J)=HBWGT1(I,J)
5158 NBWGT(2,I,J)=HBWGT2(I,J)
5159 NBWGT(3,I,J)=HBWGT3(I,J)
5160 NBWGT(4,I,J)=HBWGT4(I,J)
5161 ENDDO
5162 ENDDO
5164 DO J = JTS,MIN(JTE,JDE) !-1)
5165 DO I = ITS,MIN(ITE,IDE) !-1)
5166 AMAXVAL=0.
5167 DO N=1,4
5168 AMAXVAL=amax1(NBWGT(N,I,J),AMAXVAL)
5169 ENDDO
5170 !
5171 FLIP=.TRUE.
5172 SUM=0.0
5173 DO N=1,4
5174 IF(AMAXVAL .EQ. NBWGT(N,I,J) .AND. FLIP)THEN
5175 NBWGT(N,I,J)=1.0
5176 FLIP=.FALSE.
5177 ELSE
5178 NBWGT(N,I,J)=0.0
5179 ENDIF
5180 SUM=SUM+NBWGT(N,I,J)
5181 IF(SUM .GT. 1.0)CALL wrf_error_fatal ( "horizontal interp error - interp_hnear_nmm" )
5182 ENDDO
5185 IF((NBWGT(1,I,J)+NBWGT(2,I,J)+NBWGT(3,I,J)+NBWGT(4,I,J)) .NE. 1)THEN
5186 WRITE(0,*)'------------------------------------------------------------------------'
5187 WRITE(0,*)'FATAL: SOMETHING IS WRONG WITH THE WEIGHTS IN module_initialize_real.F'
5188 WRITE(0,*)'------------------------------------------------------------------------'
5189 STOP
5190 ENDIF
5192 ! WRITE(66,*)I,J,NBWGT(1,I,J),NBWGT(2,I,J),NBWGT(3,I,J),NBWGT(4,I,J)
5194 ENDDO
5195 ENDDO
5198 DO J=MAX(3,JTS),MIN(JTE,JDE) !-1)
5199 DO I=MAX(3,ITS),MIN(ITE,IDE) !-1)
5200 IF(MOD(JJH(I,J),2) .NE. 0)THEN ! 1,3,5,7
5201 HRES_SM(I,J) = NBWGT(1,I,J)*SM(IIH(I,J),JJH(I,J) ) &
5202 + NBWGT(2,I,J)*SM(IIH(I,J)+1, JJH(I,J) ) &
5203 + NBWGT(3,I,J)*SM(IIH(I,J), JJH(I,J)-1) &
5204 + NBWGT(4,I,J)*SM(IIH(I,J), JJH(I,J)+1)
5205 ! WRITE(68,*)I,J,SM(IIH(I,J),JJH(I,J)),SM(IIH(I,J)+1, JJH(I,J)), &
5206 ! SM(IIH(I,J), JJH(I,J)-1),SM(IIH(I,J), JJH(I,J)+1),HRES_SM(I,J)
5207 ELSE
5208 HRES_SM(I,J) = NBWGT(1,I,J)*SM(IIH(I,J), JJH(I,J) ) &
5209 + NBWGT(2,I,J)*SM(IIH(I,J)+1, JJH(I,J) ) &
5210 + NBWGT(3,I,J)*SM(IIH(I,J)+1, JJH(I,J)-1) &
5211 + NBWGT(4,I,J)*SM(IIH(I,J)+1, JJH(I,J)+1)
5213 ! WRITE(68,*)I,J,SM(IIH(I,J),JJH(I,J)),SM(IIH(I,J)+1, JJH(I,J)), &
5214 ! SM(IIH(I,J)+1, JJH(I,J)-1),SM(IIH(I,J)+1, JJH(I,J)+1),HRES_SM(I,J)
5216 ENDIF
5218 ENDDO
5219 ENDDO
5220 ! Boundary treatment in J direction
5221 DO J=MAX(3,JTS),MIN(JTE,JDE)
5222 HRES_SM(2,J)=HRES_SM(3,J)
5223 HRES_SM(1,J)=HRES_SM(2,J)
5224 END DO
5225 ! Boundary treatment in J direction and 4 corners
5226 DO I=ITS,MIN(ITE,IDE)
5227 HRES_SM(I,2)=HRES_SM(I,3)
5228 HRES_SM(I,1)=HRES_SM(I,2)
5229 END DO
5232 RETURN
5233 END SUBROUTINE G2T2H_hwrf
5234 !========================================================================================
5235 ! end gopal's doing for ocean coupling
5236 !============================================================================================
5237 #endif
5238 END MODULE module_initialize_real