Bugfix in the calculations of immediate melting and supersaturation dynamical tendenc...
[WRF.git] / main / ndown_em.F
blobb5ae20efad6b48d1002328b08e2d260f814ab4ea
1 !WRF:DRIVER_LAYER:MAIN
4 PROGRAM ndown_em
6    USE module_machine
7    USE module_domain, ONLY : domain, head_grid, alloc_and_configure_domain, &
8                              domain_clock_set, domain_clock_get, get_ijk_from_grid
9    USE module_domain_type, ONLY : program_name
10    USE module_initialize_real, ONLY : wrfu_initialize, rebalance_driver
11    USE module_integrate
12    USE module_driver_constants
13    USE module_configure, ONLY : grid_config_rec_type, model_config_rec
14    USE module_io_domain
15    USE module_utility
16    USE module_check_a_mundo
18    USE module_timing
19    USE module_wrf_error
20 #ifdef DM_PARALLEL
21    USE module_dm
22 #endif
24 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
25 !new for bc
26    USE module_bc
27    USE module_big_step_utilities_em
28    USE module_get_file_names
29 #if ( WRF_CHEM == 1 )
30 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
31 ! for chemistry
32    USE module_input_chem_data
33 !  USE module_input_chem_bioemiss
34 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
35 #endif
37    IMPLICIT NONE
38  ! interface
39    INTERFACE
40      ! mediation-supplied
41      SUBROUTINE med_read_wrf_chem_bioemiss ( grid , config_flags)
42        USE module_domain
43        TYPE (domain) grid
44        TYPE (grid_config_rec_type) config_flags
45      END SUBROUTINE med_read_wrf_chem_bioemiss
47      SUBROUTINE init_domain_constants_em_ptr ( parent , nest )
48        USE module_domain
49        USE module_configure
50        TYPE(domain), POINTER  :: parent , nest
51      END SUBROUTINE init_domain_constants_em_ptr
53      SUBROUTINE vertical_interp (nested_grid,znw_c,znu_c,cf1_c,cf2_c,cf3_c,cfn_c,cfn1_c,k_dim_c,c3h,c4h,c3f,c4f)
54          USE module_domain
55          USE module_configure
56          TYPE(domain), POINTER ::  nested_grid
57          INTEGER , INTENT (IN) :: k_dim_c
58          REAL , INTENT (IN) :: cf1_c,cf2_c,cf3_c,cfn_c,cfn1_c
59          REAL , DIMENSION(k_dim_c) , INTENT (IN) ::  znw_c,znu_c
60          REAL , DIMENSION(k_dim_c) , INTENT (IN) ::  c3h,c4h,c3f,c4f
61       END SUBROUTINE vertical_interp
64    END INTERFACE
65    INTEGER :: ids , ide , jds , jde , kds , kde
66    INTEGER :: ims , ime , jms , jme , kms , kme
67    INTEGER :: ips , ipe , jps , jpe , kps , kpe
68    INTEGER :: its , ite , jts , jte , kts , kte
69    INTEGER :: nids, nide, njds, njde, nkds, nkde,    &
70               nims, nime, njms, njme, nkms, nkme,    &
71               nips, nipe, njps, njpe, nkps, nkpe
72    INTEGER :: spec_bdy_width
73    INTEGER :: i , j , k , nvchem, nvmoist, nvscalar
74    INTEGER :: time_loop_max , time_loop
75    INTEGER :: total_time_sec , file_counter
76    INTEGER :: julyr , julday , iswater , map_proj
77    INTEGER :: icnt
79    REAL    :: dt , new_bdy_frq
80    REAL    :: gmt , cen_lat , cen_lon , dx , dy , truelat1 , truelat2 , moad_cen_lat , stand_lon
82    REAL , DIMENSION(:,:,:) , ALLOCATABLE :: ubdy3dtemp1 , vbdy3dtemp1 , tbdy3dtemp1 , pbdy3dtemp1 , qbdy3dtemp1
83    REAL , DIMENSION(:,:,:) , ALLOCATABLE :: mbdy2dtemp1
84    REAL , DIMENSION(:,:,:) , ALLOCATABLE :: ubdy3dtemp2 , vbdy3dtemp2 , tbdy3dtemp2 , pbdy3dtemp2 , qbdy3dtemp2
85    REAL , DIMENSION(:,:,:) , ALLOCATABLE :: mbdy2dtemp2
86    REAL , DIMENSION(:,:,:) , ALLOCATABLE :: cbdy3dtemp1 , cbdy3dtemp2
87    REAL , DIMENSION(:,:,:) , ALLOCATABLE :: sbdy3dtemp1 , sbdy3dtemp2
88    REAL , DIMENSION(:,:,:,:) , ALLOCATABLE :: cbdy3dtemp0
89    REAL , DIMENSION(:,:,:,:) , ALLOCATABLE :: qbdy3dtemp1_coupled, qbdy3dtemp2_coupled
90    REAL , DIMENSION(:,:,:,:) , ALLOCATABLE :: sbdy3dtemp1_coupled, sbdy3dtemp2_coupled
92    CHARACTER(LEN=19) :: start_date_char , current_date_char , end_date_char
93    CHARACTER(LEN=19) :: stopTimeStr
95 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
97    INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
99    REAL    :: time
100    INTEGER :: rc
102    INTEGER :: loop , levels_to_process
103    INTEGER , PARAMETER :: max_sanity_file_loop = 100
105    TYPE (domain) , POINTER :: keep_grid, grid_ptr, null_domain, parent_grid , nested_grid
106    TYPE (domain)           :: dummy
107    TYPE (grid_config_rec_type)              :: config_flags
108    INTEGER                 :: number_at_same_level
109    INTEGER                 :: time_step_begin_restart
111    INTEGER :: max_dom , domain_id , fid , fido, fidb , oid , idum1 , idum2 , ierr
112    INTEGER :: status_next_var
113    INTEGER :: debug_level
114    LOGICAL :: input_from_file , need_new_file
115    CHARACTER (LEN=19) :: date_string
117 #ifdef DM_PARALLEL
118    INTEGER                 :: nbytes
119    INTEGER, PARAMETER      :: configbuflen = 4* CONFIG_BUF_LEN
120    INTEGER                 :: configbuf( configbuflen )
121    LOGICAL , EXTERNAL      :: wrf_dm_on_monitor
122 #endif
124    INTEGER                 :: idsi
125    CHARACTER (LEN=256)     :: inpname , outname , bdyname
126    CHARACTER (LEN=256)     :: si_inpname
127 character *19 :: temp19
128 character *24 :: temp24 , temp24b
129 character(len=24) :: start_date_hold
131    CHARACTER (LEN=256)     :: message
132 integer :: ii
134 #include "version_decl"
135 #include "commit_decl"
137 !!!!!!!!!!!!!!!!!!!!!    mousta   
138     integer ::  n_ref_m,k_dim_c,k_dim_n
139 real :: cf1_c,cf2_c,cf3_c,cfn_c,cfn1_c
140    REAL , DIMENSION(:) , ALLOCATABLE :: znw_c,znu_c
141    REAL , DIMENSION(:) , ALLOCATABLE :: c3h,c4h,c3f,c4f
142 !!!!!!!!!!!!!!!!!!!!!!!!!!11
144    !  Interface block for routine that passes pointers and needs to know that they
145    !  are receiving pointers.
147    INTERFACE
149       SUBROUTINE med_interp_domain ( parent_grid , nested_grid )
150          USE module_domain
151          USE module_configure
152          TYPE(domain), POINTER :: parent_grid , nested_grid
153       END SUBROUTINE med_interp_domain
155       SUBROUTINE Setup_Timekeeping( parent_grid )
156          USE module_domain
157          TYPE(domain), POINTER :: parent_grid
158       END SUBROUTINE Setup_Timekeeping
160       SUBROUTINE vert_cor(parent_grid,znw_c,k_dim_c)
161          USE module_domain
162          TYPE(domain), POINTER :: parent_grid
163          integer , intent(in) :: k_dim_c
164          real , dimension(k_dim_c), INTENT(IN) :: znw_c
165       END SUBROUTINE vert_cor
166    END INTERFACE
168    !  Define the name of this program (program_name defined in module_domain)
170    program_name = "NDOWN_EM " // TRIM(release_version) // " PREPROCESSOR"
172 #ifdef DM_PARALLEL
173    CALL disable_quilting
174 #else
175    CALL wrf_error_fatal ( 'NDOWN :  HAVE TO BUILD FOR NESTING' )
176 #endif
178    !  Initialize the modules used by the WRF system.  Many of the CALLs made from the
179    !  init_modules routine are NO-OPs.  Typical initializations are: the size of a
180    !  REAL, setting the file handles to a pre-use value, defining moisture and
181    !  chemistry indices, etc.
183    CALL init_modules(1)   ! Phase 1 returns after MPI_INIT() (if it is called)
184 #ifdef NO_LEAP_CALENDAR
185    CALL WRFU_Initialize( defaultCalKind=WRFU_CAL_NOLEAP, rc=rc )
186 #else
187    CALL WRFU_Initialize( defaultCalKind=WRFU_CAL_GREGORIAN, rc=rc )
188 #endif
189    CALL init_modules(2)   ! Phase 2 resumes after MPI_INIT() (if it is called)
191    !  Get the NAMELIST data.  This is handled in the initial_config routine.  All of the
192    !  NAMELIST input variables are assigned to the model_config_rec structure.  Below,
193    !  note for parallel processing, only the monitor processor handles the raw Fortran
194    !  I/O, and then broadcasts the info to each of the other nodes.
196 #ifdef DM_PARALLEL
197    IF ( wrf_dm_on_monitor() ) THEN
198      CALL initial_config
199    ENDIF
200    CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
201    CALL wrf_dm_bcast_bytes( configbuf, nbytes )
202    CALL set_config_as_buffer( configbuf, configbuflen )
203    CALL wrf_dm_initialize
204 #else
205    CALL initial_config
206 #endif
208    CALL setup_physics_suite
209    CALL check_nml_consistency
210    CALL set_physics_rconfigs
212    !  If we are running ndown, and that is WHERE we are now, make sure that we account
213    !  for folks forgetting to say that the aux_input2 stream is in place.
215    IF ( model_config_rec%io_form_auxinput2 .EQ. 0 ) THEN
216       CALL wrf_error_fatal( 'ndown: Please set io_form_auxinput2 in the time_control namelist record (probably to 2).')
217    END IF
219 !!!!!!!!!!!!!!! mousta
220    n_ref_m = model_config_rec % vert_refine_fact
221    k_dim_c = model_config_rec % e_vert(1)
222    k_dim_n = k_dim_c
223    if (n_ref_m .ge. 2) k_dim_n = (k_dim_c - 1) *  n_ref_m + 1
224    model_config_rec % e_vert(1) = k_dim_n
225    model_config_rec % e_vert(2) = k_dim_n
226    ALLOCATE(znw_c(k_dim_c))
227    ALLOCATE(znu_c(k_dim_c))
228    ALLOCATE(c3h(k_dim_c))
229    ALLOCATE(c4h(k_dim_c))
230    ALLOCATE(c3f(k_dim_c))
231    ALLOCATE(c4f(k_dim_c))
232    WRITE ( message , FMT = '(A,3I5)' ) 'KDIM_C', k_dim_c , model_config_rec % e_vert(1) , model_config_rec % e_vert(2)
233    CALL       wrf_debug (  99,message )
234 !!!!!!!!!!!!!!! mousta
236    !  And here is an instance of using the information in the NAMELIST.  
238    CALL nl_get_debug_level ( 1, debug_level )
239    CALL set_wrf_debug_level ( debug_level )
241    !  Allocated and configure the mother domain.  Since we are in the nesting down
242    !  mode, we know a) we got a nest, and b) we only got 1 nest.
244    NULLIFY( null_domain )
246    CALL       wrf_message ( program_name )
247    CALL       wrf_message ( commit_version )
248    CALL       wrf_debug ( 100 , 'ndown_em: calling alloc_and_configure_domain coarse ' )
249    CALL alloc_and_configure_domain ( domain_id  = 1 ,                  &
250                                      grid       = head_grid ,          &
251                                      parent     = null_domain ,        &
252                                      kid        = -1                   )
254    parent_grid => head_grid
256    !  Set up time initializations.
258    CALL Setup_Timekeeping ( parent_grid )
260    CALL domain_clock_set( head_grid, &
261                           time_step_seconds=model_config_rec%interval_seconds )
262    CALL       wrf_debug ( 100 , 'ndown_em: calling model_to_grid_config_rec ' )
263    CALL model_to_grid_config_rec ( parent_grid%id , model_config_rec , config_flags )
264    CALL       wrf_debug ( 100 , 'ndown_em: calling set_scalar_indices_from_config ' )
265    CALL set_scalar_indices_from_config ( parent_grid%id , idum1, idum2 )
267    !  Initialize the I/O for WRF.
269    CALL       wrf_debug ( 100 , 'ndown_em: calling init_wrfio' )
270    CALL init_wrfio
272    !  Some of the configuration values may have been modified from the initial READ
273    !  of the NAMELIST, so we re-broadcast the configuration records.
275 #ifdef DM_PARALLEL
276    CALL get_config_as_buffer( configbuf, configbuflen, nbytes )
277    CALL wrf_dm_bcast_bytes( configbuf, nbytes )
278    CALL set_config_as_buffer( configbuf, configbuflen )
279 #endif
281    !  We need to current and starting dates for the output files.  The times need to be incremented
282    !  so that the lateral BC files are not overwritten.
284 #ifdef PLANET
285    WRITE ( start_date_char , FMT = '(I4.4,"-",I5.5,"_",I2.2,":",I2.2,":",I2.2)' ) &
286            model_config_rec%start_year  (parent_grid%id) , &
287            model_config_rec%start_day   (parent_grid%id) , &
288            model_config_rec%start_hour  (parent_grid%id) , &
289            model_config_rec%start_minute(parent_grid%id) , &
290            model_config_rec%start_second(parent_grid%id)
292    WRITE (   end_date_char , FMT = '(I4.4,"-",I5.5,"_",I2.2,":",I2.2,":",I2.2)' ) &
293            model_config_rec%  end_year  (parent_grid%id) , &
294            model_config_rec%  end_day   (parent_grid%id) , &
295            model_config_rec%  end_hour  (parent_grid%id) , &
296            model_config_rec%  end_minute(parent_grid%id) , &
297            model_config_rec%  end_second(parent_grid%id)
298 #else
299    WRITE ( start_date_char , FMT = '(I4.4,"-",I2.2,"-",I2.2,"_",I2.2,":",I2.2,":",I2.2)' ) &
300            model_config_rec%start_year  (parent_grid%id) , &
301            model_config_rec%start_month (parent_grid%id) , &
302            model_config_rec%start_day   (parent_grid%id) , &
303            model_config_rec%start_hour  (parent_grid%id) , &
304            model_config_rec%start_minute(parent_grid%id) , &
305            model_config_rec%start_second(parent_grid%id)
307    WRITE (   end_date_char , FMT = '(I4.4,"-",I2.2,"-",I2.2,"_",I2.2,":",I2.2,":",I2.2)' ) &
308            model_config_rec%  end_year  (parent_grid%id) , &
309            model_config_rec%  end_month (parent_grid%id) , &
310            model_config_rec%  end_day   (parent_grid%id) , &
311            model_config_rec%  end_hour  (parent_grid%id) , &
312            model_config_rec%  end_minute(parent_grid%id) , &
313            model_config_rec%  end_second(parent_grid%id)
314 #endif
316    !  Override stop time with value computed above.
317    CALL domain_clock_set( parent_grid, stop_timestr=end_date_char )
319    CALL geth_idts ( end_date_char , start_date_char , total_time_sec )
321    new_bdy_frq = model_config_rec%interval_seconds
322    time_loop_max = total_time_sec / model_config_rec%interval_seconds + 1
324    start_date        = start_date_char // '.0000'
325    current_date      = start_date_char // '.0000'
326    start_date_hold   = start_date_char // '.0000'
327    current_date_char = start_date_char
329    !  Get a list of available file names to try.  This fills up the eligible_file_name
330    !  array with number_of_eligible_files entries.  This routine issues a nonstandard
331    !  call (system).
333    file_counter = 1
334    need_new_file = .FALSE.
335    CALL unix_ls ( 'wrfout' , parent_grid%id )
337    !  Open the input data (wrfout_d01_xxxxxx) for reading.
338    
339    CALL wrf_debug          ( 100 , 'ndown_em main: calling open_r_dataset for ' // TRIM(eligible_file_name(file_counter)) )
340    CALL open_r_dataset     ( fid, TRIM(eligible_file_name(file_counter)) , head_grid , config_flags , "DATASET=AUXINPUT1", ierr )
341    IF ( ierr .NE. 0 ) THEN
342       WRITE( wrf_err_message , FMT='(A,A,A,I8)' ) 'program ndown: error opening ',TRIM(eligible_file_name(file_counter)), &
343                                                   ' for reading ierr=',ierr
344       CALL WRF_ERROR_FATAL ( wrf_err_message )
345    ENDIF
347    !  We know how many time periods to process, so we begin.
349    big_time_loop_thingy : DO time_loop = 1 , time_loop_max
351       !  Which date are we currently soliciting?
353       CALL geth_newdate ( date_string , start_date_char , ( time_loop - 1 ) * NINT ( new_bdy_frq) )
354       WRITE ( message , FMT = '(A,I5," ",A,A)' ) '-------->>>  Processing data: loop=',time_loop,'  date/time = ',date_string
355       CALL       wrf_debug (  99,message )
356       current_date_char = date_string
357       current_date      = date_string // '.0000'
358       start_date        = date_string // '.0000'
359       WRITE ( message , FMT = '(A,I5," ",A,A)' ) 'loopmax = ', time_loop_max, '   ending date = ',end_date_char
360       CALL       wrf_debug (  99,message )
361       CALL domain_clock_set( parent_grid, &
362                              current_timestr=current_date(1:19) )
364       !  Which times are in this file, and more importantly, are any of them the
365       !  ones that we want?  We need to loop over times in each files, loop
366       !  over files.
368       get_the_right_time : DO
369       
370          CALL wrf_get_next_time ( fid , date_string , status_next_var )
371          WRITE ( message , FMT = '(A,A,A,A,A,I5)' ) 'file date/time = ',date_string,'     desired date = ',&
372          current_date_char,'     status = ', status_next_var
373          CALL       wrf_debug (  99,message )
375          IF      (  status_next_var .NE. 0 ) THEN
376             CALL wrf_debug          ( 100 , 'ndown_em main: calling close_dataset  for ' // TRIM(eligible_file_name(file_counter)) )
377             CALL close_dataset      ( fid , config_flags , "DATASET=INPUT" )
378             file_counter = file_counter + 1
379             IF ( file_counter .GT. number_of_eligible_files ) THEN
380                WRITE( wrf_err_message , FMT='(A)' ) 'program ndown: opening too many files'
381                CALL wrf_message ( TRIM(wrf_err_message) )
382                WRITE( wrf_err_message , FMT='(A)' ) 'Usually, this is caused by trying to run ndown past the last time available d01 model output'
383                CALL wrf_message ( TRIM(wrf_err_message) )
384                WRITE( wrf_err_message , FMT='(A)' ) 'The CG model output is used to supply lateral boundary conditions'
385                CALL wrf_message ( TRIM(wrf_err_message) )
386                WRITE( wrf_err_message , FMT='(A)' ) 'The ndown program uses the start and end times, the WRF model for d01 likely used the run times option.'
387                CALL wrf_message ( TRIM(wrf_err_message) )
388                WRITE( wrf_err_message , FMT='(A)' ) 'Check the namelist.input time_control section for inconsistencies.'
389                CALL WRF_ERROR_FATAL ( wrf_err_message )
390             END IF
391             CALL wrf_debug      ( 100 , 'ndown_em main: calling open_r_dataset for ' // TRIM(eligible_file_name(file_counter)) )
392             CALL open_r_dataset ( fid, TRIM(eligible_file_name(file_counter)) , head_grid , config_flags , "DATASET=INPUT", ierr )
393             IF ( ierr .NE. 0 ) THEN
394                WRITE( wrf_err_message , FMT='(A,A,A,I8)' ) 'program ndown: error opening ',TRIM(eligible_file_name(file_counter)), &
395                                                            ' for reading ierr=',ierr
396                CALL WRF_ERROR_FATAL ( wrf_err_message )
397             ENDIF
398             CYCLE get_the_right_time
399          ELSE IF ( TRIM(date_string) .LT. TRIM(current_date_char) ) THEN
400             CYCLE get_the_right_time
401          ELSE IF ( TRIM(date_string) .EQ. TRIM(current_date_char) ) THEN
402             EXIT get_the_right_time
403          ELSE IF ( TRIM(date_string) .GT. TRIM(current_date_char) ) THEN
404             WRITE( wrf_err_message , FMT='(A,A,A,A,A)' ) 'Found ',TRIM(date_string),' before I found ',TRIM(current_date_char),'.'
405             CALL WRF_ERROR_FATAL ( wrf_err_message )
406          END IF
407       END DO get_the_right_time
409       CALL wrf_debug          ( 100 , 'wrf: calling input_history' )
410       CALL wrf_get_previous_time ( fid , date_string , status_next_var )
411       WRITE ( message , * ) 'CFB' ,cf1_c,cf2_c,cf3_c,cfn_c,cfn1_c,znw_c(1),znu_c(1)
412       CALL       wrf_debug (  99,message )
413       CALL input_history ( fid , head_grid , config_flags, ierr)
414 !!!!!!!!!!!!!1   mousta
415       cf1_c = head_grid%cf1
416       cf2_c = head_grid%cf2
417       cf3_c = head_grid%cf3
419       cfn_c = head_grid%cfn
420       cfn1_c = head_grid%cfn1
422       do k = 1,k_dim_c
423       znw_c(k) = head_grid%znw(k)
424       znu_c(k) = head_grid%znu(k)
425       enddo
426       do k = 1,k_dim_c
427       c3h(k) = head_grid%c3h(k)
428       c4h(k) = head_grid%c4h(k)
429       c3f(k) = head_grid%c3f(k)
430       c4f(k) = head_grid%c4f(k)
431       enddo
432       call vert_cor(head_grid,znw_c,k_dim_c)
433       WRITE ( message , * ) 'CFA' ,cf1_c,cf2_c,cf3_c,cfn_c,cfn1_c,znw_c(1),znu_c(1)
434       CALL       wrf_debug (  99,message )
435       WRITE ( message , * ) 'CFV' ,head_grid%cf1,head_grid%cf2,head_grid%cf3,head_grid%cfn,head_grid%cfn1,&
436       head_grid%znw(1),head_grid%znu(1) , head_grid%e_vert , parent_grid%cf1 , parent_grid%znw(1) , parent_grid%znu(1)
437       CALL       wrf_debug (  99,message )
438 !!!!!!!!!!!!!1   mousta
439       CALL wrf_debug          ( 100 , 'wrf: back from input_history' )
441       !  Get the coarse grid info for later transfer to the fine grid domain.
443       CALL wrf_get_dom_ti_integer ( fid , 'MAP_PROJ' , map_proj , 1 , icnt , ierr )
444       CALL wrf_get_dom_ti_real    ( fid , 'DX'  , dx  , 1 , icnt , ierr )
445       CALL wrf_get_dom_ti_real    ( fid , 'DY'  , dy  , 1 , icnt , ierr )
446       CALL wrf_get_dom_ti_real    ( fid , 'CEN_LAT' , cen_lat , 1 , icnt , ierr )
447       CALL wrf_get_dom_ti_real    ( fid , 'CEN_LON' , cen_lon , 1 , icnt , ierr )
448       CALL wrf_get_dom_ti_real    ( fid , 'TRUELAT1' , truelat1 , 1 , icnt , ierr )
449       CALL wrf_get_dom_ti_real    ( fid , 'TRUELAT2' , truelat2 , 1 , icnt , ierr )
450       CALL wrf_get_dom_ti_real    ( fid , 'MOAD_CEN_LAT' , moad_cen_lat , 1 , icnt , ierr )
451       CALL wrf_get_dom_ti_real    ( fid , 'STAND_LON' , stand_lon , 1 , icnt , ierr )
452 !     CALL wrf_get_dom_ti_real    ( fid , 'GMT' , gmt , 1 , icnt , ierr )
453 !     CALL wrf_get_dom_ti_integer ( fid , 'JULYR' , julyr , 1 , icnt , ierr )
454 !     CALL wrf_get_dom_ti_integer ( fid , 'JULDAY' , julday , 1 , icnt , ierr )
455       CALL wrf_get_dom_ti_integer ( fid , 'ISWATER' , iswater , 1 , icnt , ierr )
457       !  First time in, do this: allocate sapce for the fine grid, get the config flags, open the
458       !  wrfinput and wrfbdy files.  This COULD be done outside the time loop, I think, so check it
459       !  out and move it up if you can.
461       IF ( time_loop .EQ. 1 ) THEN
463          CALL       wrf_message ( program_name )
464          CALL       wrf_debug ( 100 , 'wrf: calling alloc_and_configure_domain fine ' )
465          CALL alloc_and_configure_domain ( domain_id  = 2 ,                  &
466                                            grid       = nested_grid ,        &
467                                            parent     = parent_grid ,        &
468                                            kid        = 1                   )
469    
470          CALL       wrf_debug ( 100 , 'wrf: calling model_to_grid_config_rec ' )
471          CALL model_to_grid_config_rec ( nested_grid%id , model_config_rec , config_flags )
472          CALL       wrf_debug ( 100 , 'wrf: calling set_scalar_indices_from_config ' )
473          CALL set_scalar_indices_from_config ( nested_grid%id , idum1, idum2 )
475          !  Set up time initializations for the fine grid.
477          CALL Setup_Timekeeping ( nested_grid )
478          ! Strictly speaking, nest stop time should come from model_config_rec...  
479          CALL domain_clock_get( parent_grid, stop_timestr=stopTimeStr )
480          CALL domain_clock_set( nested_grid,                        &
481                                 current_timestr=current_date(1:19), &
482                                 stop_timestr=stopTimeStr ,          &
483                                 time_step_seconds=                  &
484                                   model_config_rec%interval_seconds )
486          !  Generate an output file from this program, which will be an input file to WRF.
488          CALL nl_set_bdyfrq ( nested_grid%id , new_bdy_frq )
489          config_flags%bdyfrq = new_bdy_frq
491 #if ( WRF_CHEM == 1 )
492 nested_grid%chem_opt    = parent_grid%chem_opt
493 nested_grid%chem_in_opt = parent_grid%chem_in_opt
494 #endif
496          !  Initialize constants and 1d arrays in fine grid from the parent.
498          CALL init_domain_constants_em_ptr ( parent_grid , nested_grid )
500 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
501    
502          CALL wrf_debug          ( 100 , 'ndown_em main: calling open_w_dataset for wrfinput' )
503          CALL construct_filename1( outname , 'wrfinput' , nested_grid%id , 2 )
504          CALL open_w_dataset     ( fido, TRIM(outname) , nested_grid , config_flags , output_input , "DATASET=INPUT", ierr )
505          IF ( ierr .NE. 0 ) THEN
506             WRITE( wrf_err_message , FMT='(A,A,A,I8)' ) 'program ndown: error opening ',TRIM(outname),' for reading ierr=',ierr
507             CALL WRF_ERROR_FATAL ( wrf_err_message )
508          ENDIF
510          !  Various sizes that we need to be concerned about.
512          ids = nested_grid%sd31
513          ide = nested_grid%ed31
514          kds = nested_grid%sd32
515          kde = nested_grid%ed32
516          jds = nested_grid%sd33
517          jde = nested_grid%ed33
519          ims = nested_grid%sm31
520          ime = nested_grid%em31
521          kms = nested_grid%sm32
522          kme = nested_grid%em32
523          jms = nested_grid%sm33
524          jme = nested_grid%em33
526          ips = nested_grid%sp31
527          ipe = nested_grid%ep31
528          kps = nested_grid%sp32
529          kpe = nested_grid%ep32
530          jps = nested_grid%sp33
531          jpe = nested_grid%ep33
534 !        print *, ids , ide , jds , jde , kds , kde
535 !        print *, ims , ime , jms , jme , kms , kme
536 !        print *, ips , ipe , jps , jpe , kps , kpe
538          spec_bdy_width = model_config_rec%spec_bdy_width
540          !  This is the space needed to save the current 3d data for use in computing
541          !  the lateral boundary tendencies.
543          ALLOCATE ( ubdy3dtemp1(ims:ime,kms:kme,jms:jme) )
544          ALLOCATE ( vbdy3dtemp1(ims:ime,kms:kme,jms:jme) )
545          ALLOCATE ( tbdy3dtemp1(ims:ime,kms:kme,jms:jme) )
546          ALLOCATE ( pbdy3dtemp1(ims:ime,kms:kme,jms:jme) )
547          ALLOCATE ( qbdy3dtemp1(ims:ime,kms:kme,jms:jme) )
548          ALLOCATE ( mbdy2dtemp1(ims:ime,1:1,    jms:jme) )
549          ALLOCATE ( ubdy3dtemp2(ims:ime,kms:kme,jms:jme) )
550          ALLOCATE ( vbdy3dtemp2(ims:ime,kms:kme,jms:jme) )
551          ALLOCATE ( tbdy3dtemp2(ims:ime,kms:kme,jms:jme) )
552          ALLOCATE ( pbdy3dtemp2(ims:ime,kms:kme,jms:jme) )
553          ALLOCATE ( qbdy3dtemp2(ims:ime,kms:kme,jms:jme) )
554          ALLOCATE ( qbdy3dtemp1_coupled(ims:ime,kms:kme,jms:jme,1:num_moist) )
555          ALLOCATE ( qbdy3dtemp2_coupled(ims:ime,kms:kme,jms:jme,1:num_moist) )
556          ALLOCATE ( mbdy2dtemp2(ims:ime,1:1,    jms:jme) )
557          ALLOCATE ( cbdy3dtemp0(ims:ime,kms:kme,jms:jme,1:num_chem) )
558          ALLOCATE ( cbdy3dtemp1(ims:ime,kms:kme,jms:jme) )
559          ALLOCATE ( cbdy3dtemp2(ims:ime,kms:kme,jms:jme) )
560          ALLOCATE ( sbdy3dtemp1(ims:ime,kms:kme,jms:jme) )
561          ALLOCATE ( sbdy3dtemp2(ims:ime,kms:kme,jms:jme) )
562          ALLOCATE ( sbdy3dtemp1_coupled(ims:ime,kms:kme,jms:jme,1:num_scalar) )
563          ALLOCATE ( sbdy3dtemp2_coupled(ims:ime,kms:kme,jms:jme,1:num_scalar) )
566       END IF
568       CALL domain_clock_set( nested_grid,                        &
569                              current_timestr=current_date(1:19), &
570                              time_step_seconds=                  &
571                                model_config_rec%interval_seconds )
573       !  Do the horizontal interpolation.
575       nested_grid%imask_nostag = 1
576       nested_grid%imask_xstag = 1
577       nested_grid%imask_ystag = 1
578       nested_grid%imask_xystag = 1
580       CALL med_interp_domain ( head_grid , nested_grid )
581       WRITE ( message , * ) 'MOUSTA_L', nested_grid%mu_2(ips,jps) , nested_grid%u_2(ips,kde-2,jps)
582       CALL       wrf_debug (  99,message )
583       CALL vertical_interp (nested_grid,znw_c,znu_c,cf1_c,cf2_c,cf3_c,cfn_c,cfn1_c,k_dim_c,c3h,c4h,c3f,c4f)
584       WRITE ( message , * ) 'MOUSTA_V', nested_grid%mu_2(ips,jps) , nested_grid%u_2(ips,kde-2,jps)
585       CALL       wrf_debug (  99,message )
586       nested_grid%ht_int = nested_grid%ht
588 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
590       IF ( time_loop .EQ. 1 ) THEN
592          !  Dimension info for fine grid.
594          CALL get_ijk_from_grid (  nested_grid ,                   &
595                                    nids, nide, njds, njde, nkds, nkde,    &
596                                    nims, nime, njms, njme, nkms, nkme,    &
597                                    nips, nipe, njps, njpe, nkps, nkpe    )
599          !  Store horizontally interpolated terrain in temp location
601          CALL  copy_3d_field ( nested_grid%ht_fine , nested_grid%ht , &
602                                nids , nide , njds , njde , 1   , 1   , &
603                                nims , nime , njms , njme , 1   , 1   , &
604                                nips , nipe , njps , njpe , 1   , 1   )
606          !  Open the fine grid SI static file.
607    
608          CALL construct_filename1( si_inpname , 'wrfndi' , nested_grid%id , 2 )
609          CALL       wrf_debug ( 100 , 'med_sidata_input: calling open_r_dataset for ' // TRIM(si_inpname) )
610          CALL open_r_dataset ( idsi, TRIM(si_inpname) , nested_grid , config_flags , "DATASET=INPUT", ierr )
611          IF ( ierr .NE. 0 ) THEN
612             CALL wrf_error_fatal( 'real: error opening FG input for reading: ' // TRIM (si_inpname) )
613          END IF
615          !  Input data.
616    
617          CALL       wrf_debug ( 100 , 'ndown_em: calling input_auxinput2' )
618          CALL input_auxinput2 ( idsi , nested_grid , config_flags , ierr )
619          nested_grid%ht_input = nested_grid%ht
620    
621          !  Close this fine grid static input file.
622    
623          CALL       wrf_debug ( 100 , 'ndown_em: closing fine grid static input' )
624          CALL close_dataset ( idsi , config_flags , "DATASET=INPUT" )
626          !  Blend parent and nest field of terrain.
628          CALL blend_terrain ( nested_grid%ht_fine , nested_grid%ht , &
629                                nids , nide , njds , njde , 1   , 1   , &
630                                nims , nime , njms , njme , 1   , 1   , &
631                                nips , nipe , njps , njpe , 1   , 1   )
633          nested_grid%ht_input = nested_grid%ht
634          nested_grid%ht_int   = nested_grid%ht_fine
636          !  We need a fine grid landuse in the interpolation.  So we need to generate
637          !  that field now.
639          IF      ( ( nested_grid%ivgtyp(ips,jps) .GT. 0 ) .AND. &
640                    ( nested_grid%isltyp(ips,jps) .GT. 0 ) ) THEN
641             DO j = jps, MIN(jde-1,jpe)
642                DO i = ips, MIN(ide-1,ipe)
643                   nested_grid% vegcat(i,j) = nested_grid%ivgtyp(i,j)
644                   nested_grid%soilcat(i,j) = nested_grid%isltyp(i,j)
645                END DO
646             END DO
648          ELSE IF ( ( nested_grid% vegcat(ips,jps) .GT. 0.5 ) .AND. &
649                    ( nested_grid%soilcat(ips,jps) .GT. 0.5 ) ) THEN
650             DO j = jps, MIN(jde-1,jpe)
651                DO i = ips, MIN(ide-1,ipe)
652                   nested_grid%ivgtyp(i,j) = NINT(nested_grid% vegcat(i,j))
653                   nested_grid%isltyp(i,j) = NINT(nested_grid%soilcat(i,j))
654                END DO
655             END DO
657          ELSE
658             num_veg_cat      = SIZE ( nested_grid%landusef , DIM=2 )
659             num_soil_top_cat = SIZE ( nested_grid%soilctop , DIM=2 )
660             num_soil_bot_cat = SIZE ( nested_grid%soilcbot , DIM=2 )
661    
662             CALL land_percentages (  nested_grid%xland , &
663                                      nested_grid%landusef , nested_grid%soilctop , nested_grid%soilcbot , &
664                                      nested_grid%isltyp , nested_grid%ivgtyp , &
665                                      num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
666                                      ids , ide , jds , jde , kds , kde , &
667                                      ims , ime , jms , jme , kms , kme , &
668                                      ips , ipe , jps , jpe , kps , kpe , &
669                                      model_config_rec%iswater(nested_grid%id) )
671           END IF
673           DO j = jps, MIN(jde-1,jpe)
674             DO i = ips, MIN(ide-1,ipe)
675                nested_grid%lu_index(i,j) = nested_grid%ivgtyp(i,j)
676             END DO
677          END DO
679 #ifndef PLANET
680          CALL check_consistency ( nested_grid%ivgtyp , nested_grid%isltyp , nested_grid%landmask , &
681                                   ids , ide , jds , jde , kds , kde , &
682                                   ims , ime , jms , jme , kms , kme , &
683                                   ips , ipe , jps , jpe , kps , kpe , &
684                                   model_config_rec%iswater(nested_grid%id) )
686          CALL check_consistency2( nested_grid%ivgtyp , nested_grid%isltyp , nested_grid%landmask , &
687                                   nested_grid%tmn , nested_grid%tsk , nested_grid%sst , nested_grid%xland , &
688                                   nested_grid%tslb , nested_grid%smois , nested_grid%sh2o , &
689                                   config_flags%num_soil_layers , nested_grid%id , &
690                                   ids , ide , jds , jde , kds , kde , &
691                                   ims , ime , jms , jme , kms , kme , &
692                                   ips , ipe , jps , jpe , kps , kpe , &
693                                   model_config_rec%iswater(nested_grid%id) )
694 #endif
696       END IF
698 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
699    
700       !  We have 2 terrain elevations.  One is from input and the other is from the
701       !  the horizontal interpolation.
703       nested_grid%ht_fine = nested_grid%ht_input
704       nested_grid%ht      = nested_grid%ht_int
706       !  We have both the interpolated fields and the higher-resolution static fields.  From these
707       !  the rebalancing is now done.  Note also that the field nested_grid%ht is now from the
708       !  fine grid input file (after this call is completed).
710       CALL rebalance_driver ( nested_grid )
712       !  Different things happen during the different time loops:
713       !      first loop - write wrfinput file, close data set, copy files to holder arrays
714       !      middle loops - diff 3d/2d arrays, compute and output bc
715       !      last loop - diff 3d/2d arrays, compute and output bc, write wrfbdy file, close wrfbdy file
717       IF ( time_loop .EQ. 1 ) THEN
719          !  Set the time info.
721 !        print *,'current_date = ',current_date
722          CALL domain_clock_set( nested_grid, &
723                                 current_timestr=current_date(1:19) )
724 #if ( WRF_CHEM == 1 )
726 !         Put in chemistry data
728          IF( nested_grid%chem_opt .NE. 0 ) then
729 !           IF( nested_grid%chem_in_opt .EQ. 0 ) then
730              ! Read the chemistry data from a previous wrf forecast (wrfout file)
731               ! Generate chemistry data from a idealized vertical profile
732 !             message = 'STARTING WITH BACKGROUND CHEMISTRY '
733               CALL  wrf_message ( message )
735 !             CALL input_chem_profile ( nested_grid )
737               if(nested_grid%biomass_burn_opt == BIOMASSB) THEN
738                 message = 'READING BIOMASS BURNING EMISSIONS DATA '
739                 CALL  wrf_message ( message )
740                 CALL med_read_wrf_chem_emissopt3 ( nested_grid , config_flags)
741               end if
743               if(nested_grid%dust_opt == 1 .or. nested_grid%dmsemis_opt == 1         &
744                  .or. nested_grid%chem_opt == 300  .or. nested_grid%chem_opt == 301) then
745                  message = 'READING GOCART BG AND/OR DUST and DMS REF FIELDS'
746                  CALL  wrf_message ( message )
747                  CALL med_read_wrf_chem_gocart_bg ( nested_grid , config_flags)
748               end if
750               if( nested_grid%bio_emiss_opt .eq. 2 )then
751                  message = 'READING BEIS3.11 EMISSIONS DATA'
752                  CALL  wrf_message ( message )
753                  CALL med_read_wrf_chem_bioemiss ( nested_grid , config_flags)
754               else if( nested_grid%bio_emiss_opt == 3 ) THEN
755                  message = 'READING MEGAN 2 EMISSIONS DATA'
756                  CALL  wrf_message ( message )
757                  CALL med_read_wrf_chem_bioemiss ( nested_grid , config_flags)
758               endif
759 !           ELSE
760 !             message = 'RUNNING WITHOUT CHEMISTRY INITIALIZATION'
761 !             CALL  wrf_message ( message )
762 !           ENDIF
763          ENDIF
764 #endif
766          !  Output the first time period of the data.
767    
768          CALL output_input ( fido , nested_grid , config_flags , ierr )
770          CALL wrf_put_dom_ti_integer ( fido , 'MAP_PROJ' , map_proj , 1 , ierr )
771 !        CALL wrf_put_dom_ti_real    ( fido , 'DX'  , dx  , 1 , ierr )
772 !        CALL wrf_put_dom_ti_real    ( fido , 'DY'  , dy  , 1 , ierr )
773          CALL wrf_put_dom_ti_real    ( fido , 'CEN_LAT' , cen_lat , 1 , ierr )
774          CALL wrf_put_dom_ti_real    ( fido , 'CEN_LON' , cen_lon , 1 , ierr )
775          CALL wrf_put_dom_ti_real    ( fido , 'TRUELAT1' , truelat1 , 1 , ierr )
776          CALL wrf_put_dom_ti_real    ( fido , 'TRUELAT2' , truelat2 , 1 , ierr )
777          CALL wrf_put_dom_ti_real    ( fido , 'MOAD_CEN_LAT' , moad_cen_lat , 1 , ierr )
778          CALL wrf_put_dom_ti_real    ( fido , 'STAND_LON' , stand_lon , 1 , ierr )
779          CALL wrf_put_dom_ti_integer ( fido , 'ISWATER' , iswater , 1 , ierr )
781          !  These change if the initial time for the nest is not the same as the
782          !  first time period in the WRF output file.
783          !  Now that we know the starting date, we need to set the GMT, JULYR, and JULDAY
784          !  values for the global attributes.  This call is based on the setting of the
785          !  current_date string.
787          CALL geth_julgmt ( julyr , julday , gmt)
788          CALL nl_set_julyr  ( nested_grid%id , julyr  )
789          CALL nl_set_julday ( nested_grid%id , julday )
790          CALL nl_set_gmt    ( nested_grid%id , gmt    )
791          CALL wrf_put_dom_ti_real    ( fido , 'GMT' , gmt , 1 , ierr )
792          CALL wrf_put_dom_ti_integer ( fido , 'JULYR' , julyr , 1 , ierr )
793          CALL wrf_put_dom_ti_integer ( fido , 'JULDAY' , julday , 1 , ierr )
794 !print *,'current_date =',current_date
795 !print *,'julyr=',julyr
796 !print *,'julday=',julday
797 !print *,'gmt=',gmt
798          
799          !  Close the input (wrfout_d01_000000, for example) file.  That's right, the 
800          !  input is an output file.  Who'd've thunk.
801    
802          CALL close_dataset      ( fido , config_flags , "DATASET=INPUT" )
804          !  We need to save the 3d/2d data to compute a difference during the next loop.  Couple the
805          !  3d fields with total mu (mub + mu_2) and the stagger-specific map scale factor.
807          ! u, theta, h, scalars coupled with my, v coupled with mx
808          CALL couple ( nested_grid%mu_2 , nested_grid%mub , ubdy3dtemp1 , nested_grid%u_2                 , &
809                        'u' , nested_grid%msfuy , &
810                        nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
811                        ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
812          CALL couple ( nested_grid%mu_2 , nested_grid%mub , vbdy3dtemp1 , nested_grid%v_2                 , &
813                        'v' , nested_grid%msfvx , &
814                        nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
815                        ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
816          CALL couple ( nested_grid%mu_2 , nested_grid%mub , tbdy3dtemp1 , nested_grid%t_2                 , &
817                        't' , nested_grid%msfty , &
818                        nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
819                        ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
820          CALL couple ( nested_grid%mu_2 , nested_grid%mub , pbdy3dtemp1 , nested_grid%ph_2                , &
821                        'h' , nested_grid%msfty , &
822                        nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
823                        ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
824          DO nvmoist=PARAM_FIRST_SCALAR, num_moist
825            CALL couple ( nested_grid%mu_2 , nested_grid%mub , qbdy3dtemp1 , nested_grid%moist(ims:ime,kms:kme,jms:jme,nvmoist)    , &
826                          't' , nested_grid%msfty , &
827                          nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
828                          ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
829            qbdy3dtemp1_coupled(:,:,:,nvmoist) =  qbdy3dtemp1
830          END DO
831          DO nvscalar=PARAM_FIRST_SCALAR, num_scalar
832            CALL couple ( nested_grid%mu_2 , nested_grid%mub , sbdy3dtemp1 , nested_grid%scalar(ims:ime,kms:kme,jms:jme,nvscalar)    , &
833                          't' , nested_grid%msfty , &
834                          nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
835                          ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
836            sbdy3dtemp1_coupled(:,:,:,nvscalar) =  sbdy3dtemp1
837          END DO
839           DO j = jps , jpe
840              DO i = ips , ipe
841                 mbdy2dtemp1(i,1,j) = nested_grid%mu_2(i,j)
842              END DO
843           END DO
845          !  There are 2 components to the lateral boundaries.  First, there is the starting
846          !  point of this time period - just the outer few rows and columns.
848          CALL stuff_bdy     ( ubdy3dtemp1 , nested_grid%u_bxs, nested_grid%u_bxe,                        &
849                                             nested_grid%u_bys, nested_grid%u_bye,                        &
850                                                                      'U' ,               spec_bdy_width      , &
851                                                                            ids , ide , jds , jde , kds , kde , &
852                                                                            ims , ime , jms , jme , kms , kme , &
853                                                                            ips , ipe , jps , jpe , kps , kpe )
854          CALL stuff_bdy     ( vbdy3dtemp1 , nested_grid%v_bxs, nested_grid%v_bxe,                        &
855                                             nested_grid%v_bys, nested_grid%v_bye,                        &
856                                                                      'V' ,               spec_bdy_width      , &
857                                                                            ids , ide , jds , jde , kds , kde , &
858                                                                            ims , ime , jms , jme , kms , kme , &
859                                                                            ips , ipe , jps , jpe , kps , kpe )
860          CALL stuff_bdy     ( tbdy3dtemp1 , nested_grid%t_bxs, nested_grid%t_bxe,                        &
861                                             nested_grid%t_bys, nested_grid%t_bye,                        &
862                                                                      'T' ,               spec_bdy_width      , &
863                                                                            ids , ide , jds , jde , kds , kde , &
864                                                                            ims , ime , jms , jme , kms , kme , &
865                                                                            ips , ipe , jps , jpe , kps , kpe )
866          CALL stuff_bdy     ( pbdy3dtemp1 , nested_grid%ph_bxs, nested_grid%ph_bxe,                      &
867                                             nested_grid%ph_bys, nested_grid%ph_bye,                      &
868                                                                      'W' ,               spec_bdy_width      , &
869                                                                            ids , ide , jds , jde , kds , kde , &
870                                                                            ims , ime , jms , jme , kms , kme , &
871                                                                            ips , ipe , jps , jpe , kps , kpe )
872          DO nvmoist=PARAM_FIRST_SCALAR, num_moist
873            qbdy3dtemp1 = qbdy3dtemp1_coupled(:,:,:,nvmoist)
874            CALL stuff_bdy     ( qbdy3dtemp1 , nested_grid%moist_bxs(jms:jme,kms:kme,1:spec_bdy_width,nvmoist), &
875                                               nested_grid%moist_bxe(jms:jme,kms:kme,1:spec_bdy_width,nvmoist), &
876                                               nested_grid%moist_bys(ims:ime,kms:kme,1:spec_bdy_width,nvmoist), &
877                                               nested_grid%moist_bye(ims:ime,kms:kme,1:spec_bdy_width,nvmoist), &
878                                                                     'T' ,               spec_bdy_width      , &
879                                                                            ids , ide , jds , jde , kds , kde , &
880                                                                            ims , ime , jms , jme , kms , kme , &
881                                                                            ips , ipe , jps , jpe , kps , kpe )
882          END DO
884          DO nvscalar=PARAM_FIRST_SCALAR, num_scalar
885            sbdy3dtemp1 = sbdy3dtemp1_coupled(:,:,:,nvscalar)
886            CALL stuff_bdy     ( sbdy3dtemp1 , nested_grid%scalar_bxs(jms:jme,kms:kme,1:spec_bdy_width,nvscalar), &
887                                               nested_grid%scalar_bxe(jms:jme,kms:kme,1:spec_bdy_width,nvscalar), &
888                                               nested_grid%scalar_bys(ims:ime,kms:kme,1:spec_bdy_width,nvscalar), &
889                                               nested_grid%scalar_bye(ims:ime,kms:kme,1:spec_bdy_width,nvscalar), &
890                                                                      'T' ,                 spec_bdy_width      , &
891                                                                              ids , ide , jds , jde , kds , kde , &
892                                                                              ims , ime , jms , jme , kms , kme , &
893                                                                              ips , ipe , jps , jpe , kps , kpe )
894          END DO
896          CALL stuff_bdy     ( mbdy2dtemp1 , nested_grid%mu_bxs, nested_grid%mu_bxe,                      &
897                                             nested_grid%mu_bys, nested_grid%mu_bye,                      &
898                                                                      'M' ,               spec_bdy_width      , &
899                                                                            ids , ide , jds , jde , 1 , 1 , &
900                                                                            ims , ime , jms , jme , 1 , 1 , &
901                                                                            ips , ipe , jps , jpe , 1 , 1 )
902 #if ( WRF_CHEM == 1 )
903          do nvchem=1,num_chem
904 !        if(nvchem.eq.p_o3)then
905 !          write(0,*)'fill ch_b',cbdy3dtemp1(5,1,5),nvchem
906 !        endif
907          cbdy3dtemp1(ips:ipe,kps:kpe,jps:jpe)=nested_grid%chem(ips:ipe,kps:kpe,jps:jpe,nvchem)
908 !        if(nvchem.eq.p_o3)then
909 !          write(0,*)'fill ch_b',cbdy3dtemp1(5,1,5)
910 !        endif
911          CALL stuff_bdy     ( cbdy3dtemp1 , nested_grid%chem_bxs(jms:jme,kms:kme,1:spec_bdy_width,nvchem),                                &
912                                             nested_grid%chem_bxe(jms:jme,kms:kme,1:spec_bdy_width,nvchem),                                &
913                                             nested_grid%chem_bys(ims:ime,kms:kme,1:spec_bdy_width,nvchem),                                &
914                                             nested_grid%chem_bye(ims:ime,kms:kme,1:spec_bdy_width,nvchem),                                &
915                                                                      'T' ,               spec_bdy_width      , &
916                                                                            ids , ide , jds , jde , kds , kde , &
917                                                                            ims , ime , jms , jme , kms , kme , &
918                                                                            ips , ipe , jps , jpe , kps , kpe )
919            cbdy3dtemp0(ips:ipe,kps:kpe,jps:jpe,nvchem)=cbdy3dtemp1(ips:ipe,kps:kpe,jps:jpe)
920 !        if(nvchem.eq.p_o3)then
921 !          write(0,*)'filled ch_b',time_loop,cbdy3dtemp1(ips,kps,jps),cbdy3dtemp0(ips,kps,jps,nvchem)
922 !        endif
923          enddo
924 #endif
925       ELSE IF ( ( time_loop .GT. 1 ) .AND. ( time_loop .LT. time_loop_max ) ) THEN
927          ! u, theta, h, scalars coupled with my, v coupled with mx
928          CALL couple ( nested_grid%mu_2 , nested_grid%mub , ubdy3dtemp2 , nested_grid%u_2                 , &
929                        'u' , nested_grid%msfuy , &
930                        nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
931                        ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
932          CALL couple ( nested_grid%mu_2 , nested_grid%mub , vbdy3dtemp2 , nested_grid%v_2                 , &
933                        'v' , nested_grid%msfvx , &
934                        nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
935                        ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
936          CALL couple ( nested_grid%mu_2 , nested_grid%mub , tbdy3dtemp2 , nested_grid%t_2                 , &
937                        't' , nested_grid%msfty , &
938                        nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
939                        ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
940          CALL couple ( nested_grid%mu_2 , nested_grid%mub , pbdy3dtemp2 , nested_grid%ph_2                , &
941                        'h' , nested_grid%msfty , &
942                        nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
943                        ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
944          DO nvmoist=PARAM_FIRST_SCALAR, num_moist
945            CALL couple ( nested_grid%mu_2 , nested_grid%mub , qbdy3dtemp2 , nested_grid%moist(ims:ime,kms:kme,jms:jme,nvmoist)    , &
946                          't' , nested_grid%msfty , &
947                          nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
948                          ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
949            qbdy3dtemp2_coupled(:,:,:,nvmoist) =  qbdy3dtemp2
950          END DO
951          DO nvscalar=PARAM_FIRST_SCALAR, num_scalar
952            CALL couple ( nested_grid%mu_2 , nested_grid%mub , sbdy3dtemp2 , nested_grid%scalar(ims:ime,kms:kme,jms:jme,nvscalar)    , &
953                          't' , nested_grid%msfty , &
954                          nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
955                          ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
956            sbdy3dtemp2_coupled(:,:,:,nvscalar) =  sbdy3dtemp2
957          END DO
959           DO j = jps , jpe
960              DO i = ips , ipe
961                 mbdy2dtemp2(i,1,j) = nested_grid%mu_2(i,j)
962              END DO
963           END DO
965          !  During all of the loops after the first loop, we first compute the boundary
966          !  tendencies with the current data values and the previously save information
967          !  stored in the *bdy3dtemp1 arrays.
969          CALL stuff_bdytend ( ubdy3dtemp2 , ubdy3dtemp1 , new_bdy_frq ,                               &
970                                             nested_grid%u_btxs, nested_grid%u_btxe   ,          &
971                                             nested_grid%u_btys, nested_grid%u_btye   ,          &
972                                                                   'U'  , &
973                                                                                 spec_bdy_width      , &
974                                                                   ids , ide , jds , jde , kds , kde , &
975                                                                   ims , ime , jms , jme , kms , kme , &
976                                                                   ips , ipe , jps , jpe , kps , kpe )
977          CALL stuff_bdytend ( vbdy3dtemp2 , vbdy3dtemp1 , new_bdy_frq ,                               &
978                                             nested_grid%v_btxs, nested_grid%v_btxe   ,          &
979                                             nested_grid%v_btys, nested_grid%v_btye   ,          &
980                                                                   'V'  , &
981                                                                                 spec_bdy_width      , &
982                                                                   ids , ide , jds , jde , kds , kde , &
983                                                                   ims , ime , jms , jme , kms , kme , &
984                                                                   ips , ipe , jps , jpe , kps , kpe )
985          CALL stuff_bdytend ( tbdy3dtemp2 , tbdy3dtemp1 , new_bdy_frq ,                               &
986                                             nested_grid%t_btxs, nested_grid%t_btxe   ,          &
987                                             nested_grid%t_btys, nested_grid%t_btye   ,          &
988                                                                   'T'  , &
989                                                                                 spec_bdy_width      , &
990                                                                   ids , ide , jds , jde , kds , kde , &
991                                                                   ims , ime , jms , jme , kms , kme , &
992                                                                   ips , ipe , jps , jpe , kps , kpe )
993          CALL stuff_bdytend ( pbdy3dtemp2 , pbdy3dtemp1 , new_bdy_frq ,                               &
994                                             nested_grid%ph_btxs, nested_grid%ph_btxe   ,        &
995                                             nested_grid%ph_btys, nested_grid%ph_btye   ,        &
996                                                                   'W' , &
997                                                                                 spec_bdy_width      , &
998                                                                   ids , ide , jds , jde , kds , kde , &
999                                                                   ims , ime , jms , jme , kms , kme , &
1000                                                                   ips , ipe , jps , jpe , kps , kpe )
1001          DO nvmoist=PARAM_FIRST_SCALAR, num_moist
1002            qbdy3dtemp1 = qbdy3dtemp1_coupled(:,:,:,nvmoist)
1003            qbdy3dtemp2 = qbdy3dtemp2_coupled(:,:,:,nvmoist)
1004            CALL stuff_bdytend ( qbdy3dtemp2 , qbdy3dtemp1 , new_bdy_frq ,                               &
1005                                               nested_grid%moist_btxs(jms:jme,kms:kme,1:spec_bdy_width,nvmoist), &
1006                                               nested_grid%moist_btxe(jms:jme,kms:kme,1:spec_bdy_width,nvmoist), &
1007                                               nested_grid%moist_btys(ims:ime,kms:kme,1:spec_bdy_width,nvmoist), &
1008                                               nested_grid%moist_btye(ims:ime,kms:kme,1:spec_bdy_width,nvmoist), &
1009                                                                   'T' , &
1010                                                                                 spec_bdy_width      , &
1011                                                                   ids , ide , jds , jde , kds , kde , &
1012                                                                   ims , ime , jms , jme , kms , kme , &
1013                                                                   ips , ipe , jps , jpe , kps , kpe )
1014          END DO
1016          DO nvscalar=PARAM_FIRST_SCALAR, num_scalar
1017            sbdy3dtemp1 = sbdy3dtemp1_coupled(:,:,:,nvscalar)
1018            sbdy3dtemp2 = sbdy3dtemp2_coupled(:,:,:,nvscalar)
1019            CALL stuff_bdytend ( sbdy3dtemp2 , sbdy3dtemp1 , new_bdy_frq ,  &
1020                                               nested_grid%scalar_btxs(jms:jme,kms:kme,1:spec_bdy_width,nvscalar), &
1021                                               nested_grid%scalar_btxe(jms:jme,kms:kme,1:spec_bdy_width,nvscalar), &
1022                                               nested_grid%scalar_btys(ims:ime,kms:kme,1:spec_bdy_width,nvscalar), &
1023                                               nested_grid%scalar_btye(ims:ime,kms:kme,1:spec_bdy_width,nvscalar), &
1024                                                                  'T' , &
1025                                                                                 spec_bdy_width,       &
1026                                                                   ids , ide , jds , jde , kds , kde , &
1027                                                                   ims , ime , jms , jme , kms , kme , &
1028                                                                   ips , ipe , jps , jpe , kps , kpe )
1029          END DO
1031          CALL stuff_bdytend ( mbdy2dtemp2 , mbdy2dtemp1 , new_bdy_frq ,                               &
1032                                             nested_grid%mu_btxs, nested_grid%mu_btxe   ,        &
1033                                             nested_grid%mu_btys, nested_grid%mu_btye   ,        &
1034                                                                   'M' , &
1035                                                                                 spec_bdy_width      , &
1036                                                                   ids , ide , jds , jde , 1 , 1 , &
1037                                                                   ims , ime , jms , jme , 1 , 1 , &
1038                                                                   ips , ipe , jps , jpe , 1 , 1 )
1039 #if ( WRF_CHEM == 1 )
1040          do nvchem=1,num_chem
1041          cbdy3dtemp1(ips:ipe,kps:kpe,jps:jpe)=cbdy3dtemp0(ips:ipe,kps:kpe,jps:jpe,nvchem)
1042          cbdy3dtemp2(ips:ipe,kps:kpe,jps:jpe)=nested_grid%chem(ips:ipe,kps:kpe,jps:jpe,nvchem)
1043 !        if(nvchem.eq.p_o3)then
1044 !          write(0,*)'fill 1ch_b2',time_loop,cbdy3dtemp1(ips,kps,jps),cbdy3dtemp0(ips,kps,jps,nvchem),cbdy3dtemp2(ips,kps,jps)
1045 !        endif
1046          CALL stuff_bdytend ( cbdy3dtemp2 , cbdy3dtemp1 , new_bdy_frq ,  &
1047                                             nested_grid%chem_btxs(jms:jme,kms:kme,1:spec_bdy_width,nvchem), &
1048                                             nested_grid%chem_btxe(jms:jme,kms:kme,1:spec_bdy_width,nvchem), &
1049                                             nested_grid%chem_btys(ims:ime,kms:kme,1:spec_bdy_width,nvchem), &
1050                                             nested_grid%chem_btye(ims:ime,kms:kme,1:spec_bdy_width,nvchem), &
1051                                                                  'T' , &
1052                                                                                 spec_bdy_width      , &
1053                                                                   ids , ide , jds , jde , kds , kde , &
1054                                                                   ims , ime , jms , jme , kms , kme , &
1055                                                                   ips , ipe , jps , jpe , kps , kpe )
1056          cbdy3dtemp0(ips:ipe,kps:kpe,jps:jpe,nvchem)=cbdy3dtemp2(ips:ipe,kps:kpe,jps:jpe)
1057 !        if(nvchem.eq.p_o3)then
1058 !          write(0,*)'fill 2ch_b2',cbdy3dtemp1(ips,kps,jps),cbdy3dtemp0(ips,kps,jps,nvchem),cbdy3dtemp2(ips,kps,jps)
1059 !        endif
1060          enddo
1061 #endif
1062          IF ( time_loop .EQ. 2 ) THEN
1063    
1064             !  Generate an output file from this program, which will be an input file to WRF.
1066             CALL wrf_debug          ( 100 , 'ndown_em main: calling open_w_dataset for wrfbdy' )
1067             CALL construct_filename1( bdyname , 'wrfbdy' , nested_grid%id , 2 )
1068             CALL open_w_dataset     ( fidb, TRIM(bdyname) , nested_grid , config_flags , output_boundary , &
1069                                       "DATASET=BOUNDARY", ierr )
1070             IF ( ierr .NE. 0 ) THEN
1071                WRITE( wrf_err_message , FMT='(A,A,A,I8)' ) 'program ndown: error opening ',TRIM(bdyname),' for reading ierr=',ierr
1072                CALL WRF_ERROR_FATAL ( wrf_err_message )
1073             ENDIF
1075          END IF
1077          !  Both pieces of the boundary data are now available to be written.
1078          
1079       CALL domain_clock_set( nested_grid, &
1080                              current_timestr=current_date(1:19) )
1081       temp24= current_date
1082       temp24b=start_date_hold
1083       start_date = start_date_hold
1084       CALL geth_newdate ( temp19 , temp24b(1:19) , (time_loop-2) * model_config_rec%interval_seconds )
1085       current_date = temp19 //  '.0000'
1086       CALL geth_julgmt ( julyr , julday , gmt)
1087       CALL nl_set_julyr  ( nested_grid%id , julyr  )
1088       CALL nl_set_julday ( nested_grid%id , julday )
1089       CALL nl_set_gmt    ( nested_grid%id , gmt    )
1090       CALL wrf_put_dom_ti_real    ( fidb , 'GMT' , gmt , 1 , ierr )
1091       CALL wrf_put_dom_ti_integer ( fidb , 'JULYR' , julyr , 1 , ierr )
1092       CALL wrf_put_dom_ti_integer ( fidb , 'JULDAY' , julday , 1 , ierr )
1093       CALL domain_clock_set( nested_grid, &
1094                              current_timestr=current_date(1:19) )
1095 print *,'bdy time = ',time_loop-1,'  bdy date = ',current_date,' ',start_date
1096       CALL output_boundary ( fidb , nested_grid , config_flags , ierr )
1097       current_date = temp24
1098       start_date = temp24b
1099       CALL domain_clock_set( nested_grid, &
1100                              current_timestr=current_date(1:19) )
1102          IF ( time_loop .EQ. 2 ) THEN
1103             CALL wrf_put_dom_ti_real    ( fidb , 'BDYFRQ' , new_bdy_frq , 1 , ierr )
1104          END IF
1106          !  We need to save the 3d data to compute a difference during the next loop.  Couple the
1107          !  3d fields with total mu (mub + mu_2) and the stagger-specific map scale factor.
1108          !  We load up the boundary data again for use in the next loop.
1110           DO j = jps , jpe
1111              DO k = kps , kpe
1112                 DO i = ips , ipe
1113                    ubdy3dtemp1(i,k,j) = ubdy3dtemp2(i,k,j)
1114                    vbdy3dtemp1(i,k,j) = vbdy3dtemp2(i,k,j)
1115                    tbdy3dtemp1(i,k,j) = tbdy3dtemp2(i,k,j)
1116                    pbdy3dtemp1(i,k,j) = pbdy3dtemp2(i,k,j)
1117                    qbdy3dtemp1_coupled(i,k,j,:) = qbdy3dtemp2_coupled(i,k,j,:)
1118                    sbdy3dtemp1_coupled(i,k,j,:) = sbdy3dtemp2_coupled(i,k,j,:)
1119                 END DO
1120              END DO
1121           END DO
1123           DO j = jps , jpe
1124              DO i = ips , ipe
1125                 mbdy2dtemp1(i,1,j) = mbdy2dtemp2(i,1,j)
1126              END DO
1127           END DO
1129          !  There are 2 components to the lateral boundaries.  First, there is the starting
1130          !  point of this time period - just the outer few rows and columns.
1132          CALL stuff_bdy     ( ubdy3dtemp1 , &
1133                               nested_grid%u_bxs, nested_grid%u_bxe     ,                   &
1134                               nested_grid%u_bys, nested_grid%u_bye     ,                   &
1135                                                        'U' ,               spec_bdy_width      , &
1136                                                                            ids , ide , jds , jde , kds , kde , &
1137                                                                            ims , ime , jms , jme , kms , kme , &
1138                                                                            ips , ipe , jps , jpe , kps , kpe )
1139          CALL stuff_bdy     ( vbdy3dtemp1 , &
1140                               nested_grid%v_bxs, nested_grid%v_bxe     ,                   &
1141                               nested_grid%v_bys, nested_grid%v_bye     ,                   &
1142                                                        'V' ,               spec_bdy_width      , &
1143                                                                            ids , ide , jds , jde , kds , kde , &
1144                                                                            ims , ime , jms , jme , kms , kme , &
1145                                                                            ips , ipe , jps , jpe , kps , kpe )
1146          CALL stuff_bdy     ( tbdy3dtemp1 , &
1147                               nested_grid%t_bxs, nested_grid%t_bxe     ,                   &
1148                               nested_grid%t_bys, nested_grid%t_bye     ,                   &
1149                                                        'T' ,               spec_bdy_width      , &
1150                                                                            ids , ide , jds , jde , kds , kde , &
1151                                                                            ims , ime , jms , jme , kms , kme , &
1152                                                                            ips , ipe , jps , jpe , kps , kpe )
1153          CALL stuff_bdy     ( pbdy3dtemp1 , &
1154                               nested_grid%ph_bxs, nested_grid%ph_bxe     ,                   &
1155                               nested_grid%ph_bys, nested_grid%ph_bye     ,                   &
1156                                                        'W' ,               spec_bdy_width      , &
1157                                                                            ids , ide , jds , jde , kds , kde , &
1158                                                                            ims , ime , jms , jme , kms , kme , &
1159                                                                            ips , ipe , jps , jpe , kps , kpe )
1160          DO nvmoist=PARAM_FIRST_SCALAR, num_moist
1161            qbdy3dtemp1 = qbdy3dtemp1_coupled(:,:,:,nvmoist)
1162            CALL stuff_bdy     ( qbdy3dtemp1 , &
1163                                 nested_grid%moist_bxs(jms:jme,kms:kme,1:spec_bdy_width,nvmoist), &
1164                                 nested_grid%moist_bxe(jms:jme,kms:kme,1:spec_bdy_width,nvmoist), &
1165                                 nested_grid%moist_bys(ims:ime,kms:kme,1:spec_bdy_width,nvmoist), &
1166                                 nested_grid%moist_bye(ims:ime,kms:kme,1:spec_bdy_width,nvmoist), &
1167                                                        'T' ,               spec_bdy_width      , &
1168                                                                            ids , ide , jds , jde , kds , kde , &
1169                                                                            ims , ime , jms , jme , kms , kme , &
1170                                                                            ips , ipe , jps , jpe , kps , kpe )
1171          END DO
1173          DO nvscalar=PARAM_FIRST_SCALAR, num_scalar
1174            sbdy3dtemp1 = sbdy3dtemp1_coupled(:,:,:,nvscalar)
1175            CALL stuff_bdy     ( sbdy3dtemp1 , &
1176                                 nested_grid%scalar_bxs(jms:jme,kms:kme,1:spec_bdy_width,nvscalar), &
1177                                 nested_grid%scalar_bxe(jms:jme,kms:kme,1:spec_bdy_width,nvscalar), &
1178                                 nested_grid%scalar_bys(ims:ime,kms:kme,1:spec_bdy_width,nvscalar), &
1179                                 nested_grid%scalar_bye(ims:ime,kms:kme,1:spec_bdy_width,nvscalar), &
1180                                                       'T' ,                  spec_bdy_width,   &
1181                                                                            ids , ide , jds , jde , kds , kde , &
1182                                                                            ims , ime , jms , jme , kms , kme , &
1183                                                                            ips , ipe , jps , jpe , kps , kpe )
1184          END DO
1187 #if ( WRF_CHEM == 1 )
1188          do nvchem=1,num_chem
1189          cbdy3dtemp1(ips:ipe,kps:kpe,jps:jpe)=cbdy3dtemp0(ips:ipe,kps:kpe,jps:jpe,nvchem)
1190 !        if(nvchem.eq.p_o3)then
1191 !          write(0,*)'fill 2ch_b3',cbdy3dtemp1(ips,kps,jps),cbdy3dtemp0(ips,kps,jps,nvchem)
1192 !        endif
1193          CALL stuff_bdy     ( cbdy3dtemp1 , &
1194                               nested_grid%chem_bxs(jms:jme,kms:kme,1:spec_bdy_width,nvchem), &
1195                               nested_grid%chem_bxe(jms:jme,kms:kme,1:spec_bdy_width,nvchem),     &
1196                               nested_grid%chem_bys(ims:ime,kms:kme,1:spec_bdy_width,nvchem), &
1197                               nested_grid%chem_bye(ims:ime,kms:kme,1:spec_bdy_width,nvchem),     &
1198                                                                     'T' ,               spec_bdy_width      , &
1199                                                                            ids , ide , jds , jde , kds , kde , &
1200                                                                            ims , ime , jms , jme , kms , kme , &
1201                                                                            ips , ipe , jps , jpe , kps , kpe )
1202 !          cbdy3dtemp0(ips:ipe,kps:kpe,jps:jpe,nvchem)=cbdy3dtemp1(ips:ipe,kps:kpe,jps:jpe)
1203 !        if(nvchem.eq.p_o3)then
1204 !          write(0,*)'fill 2ch_b3',cbdy3dtemp1(ips,kps,jps),cbdy3dtemp0(ips,kps,jps,nvchem)
1205 !        endif
1206          enddo
1207 #endif
1208          CALL stuff_bdy     ( mbdy2dtemp1 , &
1209                               nested_grid%mu_bxs, nested_grid%mu_bxe    ,  &
1210                               nested_grid%mu_bys, nested_grid%mu_bye    ,  &
1211                                                                      'M' ,               spec_bdy_width      , &
1212                                                                            ids , ide , jds , jde , 1 , 1 , &
1213                                                                            ims , ime , jms , jme , 1 , 1 , &
1214                                                                            ips , ipe , jps , jpe , 1 , 1 )
1216       ELSE IF ( time_loop .EQ. time_loop_max ) THEN
1218          ! u, theta, h, scalars coupled with my, v coupled with mx
1219          CALL couple ( nested_grid%mu_2 , nested_grid%mub , ubdy3dtemp2 , nested_grid%u_2                 , &
1220                        'u' , nested_grid%msfuy , &
1221                        nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
1222                        ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
1223          CALL couple ( nested_grid%mu_2 , nested_grid%mub , vbdy3dtemp2 , nested_grid%v_2                 , &
1224                        'v' , nested_grid%msfvx , &
1225                        nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
1226                        ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
1227          CALL couple ( nested_grid%mu_2 , nested_grid%mub , tbdy3dtemp2 , nested_grid%t_2                 , &
1228                        't' , nested_grid%msfty , &
1229                        nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
1230                        ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
1231          CALL couple ( nested_grid%mu_2 , nested_grid%mub , pbdy3dtemp2 , nested_grid%ph_2                , &
1232                        'h' , nested_grid%msfty , &
1233                        nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
1234                        ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
1235          DO nvmoist=PARAM_FIRST_SCALAR, num_moist
1236            CALL couple ( nested_grid%mu_2 , nested_grid%mub , qbdy3dtemp2 , nested_grid%moist(ims:ime,kms:kme,jms:jme,nvmoist)    , &
1237                          't' , nested_grid%msfty , &
1238                          nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
1239                          ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
1240            qbdy3dtemp2_coupled(:,:,:,nvmoist) = qbdy3dtemp2
1241          END DO
1242          DO nvscalar=PARAM_FIRST_SCALAR, num_scalar
1243            CALL couple ( nested_grid%mu_2 , nested_grid%mub , sbdy3dtemp2 , nested_grid%scalar(ims:ime,kms:kme,jms:jme,nvscalar)    , &
1244                          't' , nested_grid%msfty , &
1245                          nested_grid%c1h, nested_grid%c2h, nested_grid%c1f, nested_grid%c2f, &
1246                          ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, ips, ipe, jps, jpe, kps, kpe )
1247            sbdy3dtemp2_coupled(:,:,:,nvscalar) = sbdy3dtemp2
1248          END DO
1249          mbdy2dtemp2(:,1,:) = nested_grid%mu_2(:,:)
1251          !  During all of the loops after the first loop, we first compute the boundary
1252          !  tendencies with the current data values and the previously save information
1253          !  stored in the *bdy3dtemp1 arrays.
1254 #if ( WRF_CHEM == 1 )
1255          do nvchem=1,num_chem
1256          cbdy3dtemp1(ips:ipe,kps:kpe,jps:jpe)=cbdy3dtemp0(ips:ipe,kps:kpe,jps:jpe,nvchem)
1257          cbdy3dtemp2(ips:ipe,kps:kpe,jps:jpe)=nested_grid%chem(ips:ipe,kps:kpe,jps:jpe,nvchem)
1258 !        if(nvchem.eq.p_o3)then
1259 !          write(0,*)'fill 1ch_b4',cbdy3dtemp1(ips,kps,jps),cbdy3dtemp0(ips,kps,jps,nvchem),cbdy3dtemp2(ips,kps,jps)
1260 !        endif
1261          CALL stuff_bdytend ( cbdy3dtemp2 , cbdy3dtemp1 , new_bdy_frq ,  &
1262                               nested_grid%chem_btxs(jms:jme,kms:kme,1:spec_bdy_width,nvchem),  &
1263                               nested_grid%chem_btxe(jms:jme,kms:kme,1:spec_bdy_width,nvchem), &
1264                               nested_grid%chem_btys(ims:ime,kms:kme,1:spec_bdy_width,nvchem),  &
1265                               nested_grid%chem_btye(ims:ime,kms:kme,1:spec_bdy_width,nvchem), &
1266                                                                   'T' , &
1267                                                                                 spec_bdy_width      , &
1268                                                                   ids , ide , jds , jde , kds , kde , &
1269                                                                   ims , ime , jms , jme , kms , kme , &
1270                                                                   ips , ipe , jps , jpe , kps , kpe )
1271          cbdy3dtemp0(ips:ipe,kps:kpe,jps:jpe,nvchem)=cbdy3dtemp2(ips:ipe,kps:kpe,jps:jpe)
1272 !        if(nvchem.eq.p_o3)then
1273 !          write(0,*)'fill 2ch_b4',cbdy3dtemp1(ips,kps,jps),cbdy3dtemp0(ips,kps,jps,nvchem),cbdy3dtemp2(ips,kps,jps)
1274 !        endif
1275          enddo
1276 #endif
1278          CALL stuff_bdytend ( ubdy3dtemp2 , ubdy3dtemp1 , new_bdy_frq , &
1279                               nested_grid%u_btxs  , nested_grid%u_btxe , &
1280                               nested_grid%u_btys  , nested_grid%u_btye , &
1281                                                              'U'  , &
1282                                                                                 spec_bdy_width      , &
1283                                                                   ids , ide , jds , jde , kds , kde , &
1284                                                                   ims , ime , jms , jme , kms , kme , &
1285                                                                   ips , ipe , jps , jpe , kps , kpe )
1286          CALL stuff_bdytend ( vbdy3dtemp2 , vbdy3dtemp1 , new_bdy_frq , &
1287                               nested_grid%v_btxs  , nested_grid%v_btxe , &
1288                               nested_grid%v_btys  , nested_grid%v_btye , &
1289                                                              'V'  , &
1290                                                                                 spec_bdy_width      , &
1291                                                                   ids , ide , jds , jde , kds , kde , &
1292                                                                   ims , ime , jms , jme , kms , kme , &
1293                                                                   ips , ipe , jps , jpe , kps , kpe )
1294          CALL stuff_bdytend ( tbdy3dtemp2 , tbdy3dtemp1 , new_bdy_frq , &
1295                               nested_grid%t_btxs  , nested_grid%t_btxe , &
1296                               nested_grid%t_btys  , nested_grid%t_btye , &
1297                                                              'T'  , &
1298                                                                                 spec_bdy_width      , &
1299                                                                   ids , ide , jds , jde , kds , kde , &
1300                                                                   ims , ime , jms , jme , kms , kme , &
1301                                                                   ips , ipe , jps , jpe , kps , kpe )
1302          CALL stuff_bdytend ( pbdy3dtemp2 , pbdy3dtemp1 , new_bdy_frq , &
1303                               nested_grid%ph_btxs  , nested_grid%ph_btxe , &
1304                               nested_grid%ph_btys  , nested_grid%ph_btye , &
1305                                                              'W' , &
1306                                                                                 spec_bdy_width      , &
1307                                                                   ids , ide , jds , jde , kds , kde , &
1308                                                                   ims , ime , jms , jme , kms , kme , &
1309                                                                   ips , ipe , jps , jpe , kps , kpe )
1310          DO nvmoist=PARAM_FIRST_SCALAR, num_moist
1311            qbdy3dtemp1 = qbdy3dtemp1_coupled(:,:,:,nvmoist)
1312            qbdy3dtemp2 = qbdy3dtemp2_coupled(:,:,:,nvmoist)
1313            CALL stuff_bdytend ( qbdy3dtemp2 , qbdy3dtemp1 , new_bdy_frq , &
1314                                 nested_grid%moist_btxs(jms:jme,kms:kme,1:spec_bdy_width,nvmoist) , &
1315                                 nested_grid%moist_btxe(jms:jme,kms:kme,1:spec_bdy_width,nvmoist) , &
1316                                 nested_grid%moist_btys(ims:ime,kms:kme,1:spec_bdy_width,nvmoist) , &
1317                                 nested_grid%moist_btye(ims:ime,kms:kme,1:spec_bdy_width,nvmoist) , &
1318                                                              'T' , &
1319                                                                                 spec_bdy_width      , &
1320                                                                   ids , ide , jds , jde , kds , kde , &
1321                                                                   ims , ime , jms , jme , kms , kme , &
1322                                                                   ips , ipe , jps , jpe , kps , kpe )
1323          END DO
1325          DO nvscalar=PARAM_FIRST_SCALAR, num_scalar
1326            sbdy3dtemp1 = sbdy3dtemp1_coupled(:,:,:,nvscalar)
1327            sbdy3dtemp2 = sbdy3dtemp2_coupled(:,:,:,nvscalar)
1328            CALL stuff_bdytend ( sbdy3dtemp2 , sbdy3dtemp1 , new_bdy_frq ,  &
1329                               nested_grid%scalar_btxs(jms:jme,kms:kme,1:spec_bdy_width,nvscalar), &
1330                               nested_grid%scalar_btxe(jms:jme,kms:kme,1:spec_bdy_width,nvscalar), &
1331                               nested_grid%scalar_btys(ims:ime,kms:kme,1:spec_bdy_width,nvscalar), &
1332                               nested_grid%scalar_btye(ims:ime,kms:kme,1:spec_bdy_width,nvscalar), &
1333                                                                   'T' , &
1334                                                                                 spec_bdy_width ,  &
1335                                                                   ids , ide , jds , jde , kds , kde , &
1336                                                                   ims , ime , jms , jme , kms , kme , &
1337                                                                   ips , ipe , jps , jpe , kps , kpe )
1338          END DO  
1340          CALL stuff_bdytend ( mbdy2dtemp2 , mbdy2dtemp1 , new_bdy_frq , &
1341                               nested_grid%mu_btxs  , nested_grid%mu_btxe , &
1342                               nested_grid%mu_btys  , nested_grid%mu_btye , &
1343                                                              'M' , &
1344                                                                                 spec_bdy_width      , &
1345                                                                   ids , ide , jds , jde , 1 , 1 , &
1346                                                                   ims , ime , jms , jme , 1 , 1 , &
1347                                                                   ips , ipe , jps , jpe , 1 , 1 )
1349          IF ( time_loop .EQ. 2 ) THEN
1350    
1351             !  Generate an output file from this program, which will be an input file to WRF.
1353             CALL wrf_debug          ( 100 , 'ndown_em main: calling open_w_dataset for wrfbdy' )
1354             CALL construct_filename1( bdyname , 'wrfbdy' , nested_grid%id , 2 )
1355             CALL open_w_dataset     ( fidb, TRIM(bdyname) , nested_grid , config_flags , output_boundary , &
1356                                       "DATASET=BOUNDARY", ierr )
1357             IF ( ierr .NE. 0 ) THEN
1358                WRITE( wrf_err_message , FMT='(A,A,A,I8)' ) 'program ndown: error opening ',TRIM(bdyname),' for reading ierr=',ierr
1359                CALL WRF_ERROR_FATAL ( wrf_err_message )
1360             ENDIF
1362          END IF
1364          !  Both pieces of the boundary data are now available to be written.
1366       CALL domain_clock_set( nested_grid, &
1367                              current_timestr=current_date(1:19) )
1368       temp24= current_date
1369       temp24b=start_date_hold
1370       start_date = start_date_hold
1371       CALL geth_newdate ( temp19 , temp24b(1:19) , (time_loop-2) * model_config_rec%interval_seconds )
1372       current_date = temp19 //  '.0000'
1373       CALL geth_julgmt ( julyr , julday , gmt)
1374       CALL nl_set_julyr  ( nested_grid%id , julyr  )
1375       CALL nl_set_julday ( nested_grid%id , julday )
1376       CALL nl_set_gmt    ( nested_grid%id , gmt    )
1377       CALL wrf_put_dom_ti_real    ( fidb , 'GMT' , gmt , 1 , ierr )
1378       CALL wrf_put_dom_ti_integer ( fidb , 'JULYR' , julyr , 1 , ierr )
1379       CALL wrf_put_dom_ti_integer ( fidb , 'JULDAY' , julday , 1 , ierr )
1380       CALL domain_clock_set( nested_grid, &
1381                              current_timestr=current_date(1:19) )
1382       CALL output_boundary ( fidb , nested_grid , config_flags , ierr )
1383       current_date = temp24
1384       start_date = temp24b
1385       CALL domain_clock_set( nested_grid, &
1386                              current_timestr=current_date(1:19) )
1388          IF ( time_loop .EQ. 2 ) THEN
1389             CALL wrf_put_dom_ti_real    ( fidb , 'BDYFRQ' , new_bdy_frq , 1 , ierr )
1390          END IF
1392          !  Since this is the last time through here, we need to close the boundary file.
1394          CALL model_to_grid_config_rec ( nested_grid%id , model_config_rec , config_flags )
1395          CALL close_dataset ( fidb , config_flags , "DATASET=BOUNDARY" )
1398       END IF
1400       !  Process which time now?
1402    END DO big_time_loop_thingy
1403    DEALLOCATE(znw_c)
1404    DEALLOCATE(znu_c)
1405    CALL model_to_grid_config_rec ( parent_grid%id , model_config_rec , config_flags )
1406    CALL med_shutdown_io ( parent_grid , config_flags )
1408    CALL wrf_debug ( 0 , 'ndown_em: SUCCESS COMPLETE NDOWN_EM INIT' )
1410    CALL wrf_shutdown
1412    CALL WRFU_Finalize( rc=rc )
1414 END PROGRAM ndown_em
1416 SUBROUTINE land_percentages ( xland , &
1417                               landuse_frac , soil_top_cat , soil_bot_cat , &
1418                               isltyp , ivgtyp , &
1419                               num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1420                               ids , ide , jds , jde , kds , kde , &
1421                               ims , ime , jms , jme , kms , kme , &
1422                               its , ite , jts , jte , kts , kte , &
1423                               iswater )
1424    USE module_soil_pre
1426    IMPLICIT NONE
1428    INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
1429                            ims , ime , jms , jme , kms , kme , &
1430                            its , ite , jts , jte , kts , kte , &
1431                            iswater
1433    INTEGER , INTENT(IN) :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
1434    REAL , DIMENSION(ims:ime,1:num_veg_cat,jms:jme) , INTENT(INOUT):: landuse_frac
1435    REAL , DIMENSION(ims:ime,1:num_soil_top_cat,jms:jme) , INTENT(IN):: soil_top_cat
1436    REAL , DIMENSION(ims:ime,1:num_soil_bot_cat,jms:jme) , INTENT(IN):: soil_bot_cat
1437    INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: isltyp , ivgtyp
1438    REAL , DIMENSION(ims:ime,jms:jme) , INTENT(OUT) :: xland
1440    CALL process_percent_cat_new ( xland , &
1441                                   landuse_frac , soil_top_cat , soil_bot_cat , &
1442                                   isltyp , ivgtyp , &
1443                                   num_veg_cat , num_soil_top_cat , num_soil_bot_cat , &
1444                                   ids , ide , jds , jde , kds , kde , &
1445                                   ims , ime , jms , jme , kms , kme , &
1446                                   its , ite , jts , jte , kts , kte , &
1447                                   iswater )
1449 END SUBROUTINE land_percentages
1451 SUBROUTINE check_consistency ( ivgtyp , isltyp , landmask , &
1452                                   ids , ide , jds , jde , kds , kde , &
1453                                   ims , ime , jms , jme , kms , kme , &
1454                                   its , ite , jts , jte , kts , kte , &
1455                                   iswater )
1457    IMPLICIT NONE
1459    INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
1460                            ims , ime , jms , jme , kms , kme , &
1461                            its , ite , jts , jte , kts , kte , &
1462                            iswater
1463    INTEGER , DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: isltyp , ivgtyp
1464    REAL    , DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: landmask
1466    LOGICAL :: oops
1467    INTEGER :: oops_count , i , j
1469    oops = .FALSE.
1470    oops_count = 0
1472    DO j = jts, MIN(jde-1,jte)
1473       DO i = its, MIN(ide-1,ite)
1474          IF ( ( ( landmask(i,j) .LT. 0.5 ) .AND. ( ivgtyp(i,j) .NE. iswater ) ) .OR. &
1475               ( ( landmask(i,j) .GT. 0.5 ) .AND. ( ivgtyp(i,j) .EQ. iswater ) ) ) THEN
1476             print *,'mismatch in landmask and veg type'
1477             print *,'i,j=',i,j, '  landmask =',NINT(landmask(i,j)),'  ivgtyp=',ivgtyp(i,j)
1478             oops = .TRUE.
1479             oops_count = oops_count + 1
1480 landmask(i,j) = 0
1481 ivgtyp(i,j)=16
1482 isltyp(i,j)=14
1483          END IF
1484       END DO
1485    END DO
1487    IF ( oops ) THEN
1488       CALL wrf_debug( 0, 'mismatch in check_consistency, turned to water points, be careful' )
1489    END IF
1491 END SUBROUTINE check_consistency
1493 SUBROUTINE check_consistency2( ivgtyp , isltyp , landmask , &
1494                                tmn , tsk , sst , xland , &
1495                                tslb , smois , sh2o , &
1496                                num_soil_layers , id , &
1497                                ids , ide , jds , jde , kds , kde , &
1498                                ims , ime , jms , jme , kms , kme , &
1499                                its , ite , jts , jte , kts , kte , &
1500                                iswater )
1502    USE module_configure
1503    USE module_optional_input
1505    INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , &
1506                            ims , ime , jms , jme , kms , kme , &
1507                            its , ite , jts , jte , kts , kte
1508    INTEGER , INTENT(IN) :: num_soil_layers , id
1510    INTEGER , DIMENSION(ims:ime,jms:jme) :: ivgtyp , isltyp
1511    REAL    , DIMENSION(ims:ime,jms:jme) :: landmask , tmn , tsk , sst , xland
1512    REAL    , DIMENSION(ims:ime,num_soil_layers,jms:jme) :: tslb , smois , sh2o
1514    INTEGER :: oops1 , oops2
1515    INTEGER :: i , j , k
1517       fix_tsk_tmn : SELECT CASE ( model_config_rec%sf_surface_physics(id) )
1519          CASE ( SLABSCHEME , LSMSCHEME , RUCLSMSCHEME )
1520             DO j = jts, MIN(jde-1,jte)
1521                DO i = its, MIN(ide-1,ite)
1522                   IF ( ( landmask(i,j) .LT. 0.5 ) .AND. ( flag_sst .EQ. 1 ) ) THEN
1523                      tmn(i,j) = sst(i,j)
1524                      tsk(i,j) = sst(i,j)
1525                   ELSE IF ( landmask(i,j) .LT. 0.5 ) THEN
1526                      tmn(i,j) = tsk(i,j)
1527                   END IF
1528                END DO
1529             END DO
1530       END SELECT fix_tsk_tmn
1532       !  Is the TSK reasonable?
1534       DO j = jts, MIN(jde-1,jte)
1535          DO i = its, MIN(ide-1,ite)
1536             IF ( tsk(i,j) .LT. 170 .or. tsk(i,j) .GT. 400. ) THEN
1537                print *,'error in the TSK'
1538                print *,'i,j=',i,j
1539                print *,'landmask=',landmask(i,j)
1540                print *,'tsk, sst, tmn=',tsk(i,j),sst(i,j),tmn(i,j)
1541                if(tmn(i,j).gt.170. .and. tmn(i,j).lt.400.)then
1542                   tsk(i,j)=tmn(i,j)
1543                else if(sst(i,j).gt.170. .and. sst(i,j).lt.400.)then
1544                   tsk(i,j)=sst(i,j)
1545                else
1546                   CALL wrf_error_fatal ( 'TSK unreasonable' )
1547                end if
1548             END IF
1549          END DO
1550       END DO
1552       !  Is the TMN reasonable?
1554       DO j = jts, MIN(jde-1,jte)
1555          DO i = its, MIN(ide-1,ite)
1556             IF ( ( ( tmn(i,j) .LT. 170. ) .OR. ( tmn(i,j) .GT. 400. ) ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
1557                   print *,'error in the TMN'
1558                   print *,'i,j=',i,j
1559                   print *,'landmask=',landmask(i,j)
1560                   print *,'tsk, sst, tmn=',tsk(i,j),sst(i,j),tmn(i,j)
1561                if(tsk(i,j).gt.170. .and. tsk(i,j).lt.400.)then
1562                   tmn(i,j)=tsk(i,j)
1563                else if(sst(i,j).gt.170. .and. sst(i,j).lt.400.)then
1564                   tmn(i,j)=sst(i,j)
1565                else
1566                   CALL wrf_error_fatal ( 'TMN unreasonable' )
1567                endif
1568             END IF
1569          END DO
1570       END DO
1572       !  Is the TSLB reasonable?
1574       DO j = jts, MIN(jde-1,jte)
1575          DO i = its, MIN(ide-1,ite)
1576             IF ( ( ( tslb(i,1,j) .LT. 170. ) .OR. ( tslb(i,1,j) .GT. 400. ) ) .AND. ( landmask(i,j) .GT. 0.5 ) ) THEN
1577                   print *,'error in the TSLB'
1578                   print *,'i,j=',i,j
1579                   print *,'landmask=',landmask(i,j)
1580                   print *,'tsk, sst, tmn=',tsk(i,j),sst(i,j),tmn(i,j)
1581                   print *,'tslb = ',tslb(i,:,j)
1582                   print *,'old smois = ',smois(i,:,j)
1583                   DO l = 1 , num_soil_layers
1584                      sh2o(i,l,j) = 0.0
1585                   END DO
1586                   DO l = 1 , num_soil_layers
1587                      smois(i,l,j) = 0.3
1588                   END DO
1589                   if(tsk(i,j).gt.170. .and. tsk(i,j).lt.400.)then
1590                      DO l = 1 , num_soil_layers
1591                         tslb(i,l,j)=tsk(i,j)
1592                      END DO
1593                   else if(sst(i,j).gt.170. .and. sst(i,j).lt.400.)then
1594                      DO l = 1 , num_soil_layers
1595                         tslb(i,l,j)=sst(i,j)
1596                      END DO
1597                   else if(tmn(i,j).gt.170. .and. tmn(i,j).lt.400.)then
1598                      DO l = 1 , num_soil_layers
1599                         tslb(i,l,j)=tmn(i,j)
1600                      END DO
1601                   else
1602                      CALL wrf_error_fatal ( 'TSLB unreasonable' )
1603                   endif
1604             END IF
1605          END DO
1606       END DO
1608       !  Let us make sure (again) that the landmask and the veg/soil categories match.
1610 oops1=0
1611 oops2=0
1612       DO j = jts, MIN(jde-1,jte)
1613          DO i = its, MIN(ide-1,ite)
1614             IF ( ( ( landmask(i,j) .LT. 0.5 ) .AND. ( ivgtyp(i,j) .NE. iswater .OR. isltyp(i,j) .NE. 14 ) ) .OR. &
1615                  ( ( landmask(i,j) .GT. 0.5 ) .AND. ( ivgtyp(i,j) .EQ. iswater .OR. isltyp(i,j) .EQ. 14 ) ) ) THEN
1616                IF ( tslb(i,1,j) .GT. 1. ) THEN
1617 oops1=oops1+1
1618                   ivgtyp(i,j) = 5
1619                   isltyp(i,j) = 8
1620                   landmask(i,j) = 1
1621                   xland(i,j) = 1
1622                ELSE IF ( sst(i,j) .GT. 1. ) THEN
1623 oops2=oops2+1
1624                   ivgtyp(i,j) = iswater
1625                   isltyp(i,j) = 14
1626                   landmask(i,j) = 0
1627                   xland(i,j) = 2
1628                ELSE
1629                   print *,'the landmask and soil/veg cats do not match'
1630                   print *,'i,j=',i,j
1631                   print *,'landmask=',landmask(i,j)
1632                   print *,'ivgtyp=',ivgtyp(i,j)
1633                   print *,'isltyp=',isltyp(i,j)
1634                   print *,'iswater=', iswater
1635                   print *,'tslb=',tslb(i,:,j)
1636                   print *,'sst=',sst(i,j)
1637                   CALL wrf_error_fatal ( 'mismatch_landmask_ivgtyp' )
1638                END IF
1639             END IF
1640          END DO
1641       END DO
1642 if (oops1.gt.0) then
1643 print *,'points artificially set to land : ',oops1
1644 endif
1645 if(oops2.gt.0) then
1646 print *,'points artificially set to water: ',oops2
1647 endif
1649 END SUBROUTINE check_consistency2
1650 !!!!!!!!!!!!!!!!!!!!!!!!!!!11
1651       SUBROUTINE vert_cor(parent,znw_c,k_dim_c)
1652          USE module_domain
1653          USE module_model_constants
1654    IMPLICIT NONE
1655          TYPE(domain), POINTER :: parent
1656          integer , intent(in) :: k_dim_c
1657          real , dimension(k_dim_c), INTENT(IN) :: znw_c
1659        integer :: kde_c , kde_n ,n_refine,ii,kkk,k
1660        INTEGER :: ids, ide, jds, jde, kds, kde, &
1661                   ims, ime, jms, jme, kms, kme, &
1662                   ips, ipe, jps, jpe, kps, kpe
1663        real :: dznw_m,cof1,cof2
1665         kde_c = k_dim_c
1666         kde_n = parent%e_vert
1667         n_refine = parent % vert_refine_fact
1670          kkk = 0
1671          do k = 1 , kde_c-1
1672          dznw_m = znw_c(k+1) - znw_c(k)
1673          do ii = 1,n_refine
1674          kkk = kkk + 1
1675          parent%znw(kkk) = znw_c(k) + float(ii-1)/float(n_refine)*dznw_m
1676          enddo
1677          enddo
1678          parent%znw(kde_n) = znw_c(kde_c)
1679          parent%znw(1) = znw_c(1)
1681       DO k=1, kde_n-1
1682          parent%dnw(k) = parent%znw(k+1) - parent%znw(k)
1683          parent%rdnw(k) = 1./parent%dnw(k)
1684          parent%znu(k) = 0.5*(parent%znw(k+1)+parent%znw(k))
1685       END DO
1687       DO k=2, kde_n-1
1688          parent%dn(k) = 0.5*(parent%dnw(k)+parent%dnw(k-1))
1689          parent%rdn(k) = 1./parent%dn(k)
1690          parent%fnp(k) = .5* parent%dnw(k  )/parent%dn(k)
1691          parent%fnm(k) = .5* parent%dnw(k-1)/parent%dn(k)
1692       END DO
1694   cof1 = (2.*parent%dn(2)+parent%dn(3))/(parent%dn(2)+parent%dn(3))*parent%dnw(1)/parent%dn(2)
1695   cof2 =     parent%dn(2)        /(parent%dn(2)+parent%dn(3))*parent%dnw(1)/parent%dn(3)
1697       parent%cf1  = parent%fnp(2) + cof1
1698       parent%cf2  = parent%fnm(2) - cof1 - cof2
1699       parent%cf3  = cof2
1701       parent%cfn  = (.5*parent%dnw(kde_n-1)+parent%dn(kde_n-1))/parent%dn(kde_n-1)
1702       parent%cfn1 = -.5*parent%dnw(kde_n-1)/parent%dn(kde_n-1)
1704       ids = parent%sd31
1705       ide = parent%ed31
1706       kds = parent%sd32
1707       jds = parent%sd33
1708       jde = parent%ed33
1710       ims = parent%sm31
1711       ime = parent%em31
1712       kms = parent%sm32
1713       jms = parent%sm33
1714       jme = parent%em33
1716       ips = parent%sp31
1717       ipe = parent%ep31
1718       kps = parent%sp32
1719       jps = parent%sp33
1720       jpe = parent%ep33
1722       CALL compute_vcoord_1d_coeffs ( parent%ht, parent%etac, parent%znw, &
1723                                       parent%hybrid_opt, &
1724                                       r_d, g, p1000mb, &
1725                                       parent%p_top, parent%p00, parent%t00, parent%tlp, &
1726                                       ids, ide, jds, jde, kds, kde_n, &
1727                                       ims, ime, jms, jme, kms, kde_n, &
1728                                       ips, ipe, jps, jpe, kps, kde_n, &
1729                                       parent%znu, &
1730                                       parent%c1f, parent%c2f, parent%c3f, parent%c4f, &
1731                                       parent%c1h, parent%c2h, parent%c3h, parent%c4h )
1733       END SUBROUTINE vert_cor
1736 SUBROUTINE init_domain_constants_em_ptr ( parent , nest )
1737    USE module_domain
1738    USE module_configure
1739    IMPLICIT NONE
1740    TYPE(domain), POINTER  :: parent , nest
1741    INTERFACE
1742    SUBROUTINE init_domain_constants_em ( parent , nest )
1743       USE module_domain
1744       USE module_configure
1745       TYPE(domain)  :: parent , nest
1746    END SUBROUTINE init_domain_constants_em
1747    END INTERFACE
1748    CALL init_domain_constants_em ( parent , nest )
1749 END SUBROUTINE init_domain_constants_em_ptr
1752      SUBROUTINE vertical_interp (parent_grid,znw_c,znu_c,cf1_c,cf2_c,cf3_c,cfn_c,cfn1_c,k_dim_c,c3h,c4h,c3f,c4f)
1753          USE module_domain
1754          USE module_configure
1755    IMPLICIT NONE
1756          TYPE(domain), POINTER ::  parent_grid
1757          INTEGER , INTENT (IN) :: k_dim_c
1758          REAL , INTENT (IN) :: cf1_c,cf2_c,cf3_c,cfn_c,cfn1_c
1759          REAL , DIMENSION(k_dim_c) , INTENT (IN) ::  znw_c,znu_c
1760          REAL , DIMENSION(k_dim_c) , INTENT (IN) ::  c3h,c4h,c3f,c4f
1762        integer :: kde_c , kde_n ,n_refine,ii,kkk
1763        integer :: i , j, k , itrace
1764        real :: p_top_m  , p_surf_m , mu_m , hsca_m , pre_c ,pre_n
1766        real, allocatable, dimension(:) ::  alt_w_c , alt_u_c ,pro_w_c , pro_u_c
1767        real, allocatable, dimension(:) ::  alt_w_n , alt_u_n ,pro_w_n , pro_u_n
1769    INTEGER :: nids, nide, njds, njde, nkds, nkde, &
1770               nims, nime, njms, njme, nkms, nkme, &
1771               nips, nipe, njps, njpe, nkps, nkpe
1773       
1774          hsca_m = 6.7
1775          n_refine = model_config_rec % vert_refine_fact
1776          kde_c = k_dim_c
1777          kde_n = parent_grid%e_vert
1779          CALL get_ijk_from_grid ( parent_grid , &
1780                                    nids, nide, njds, njde, nkds, nkde, &
1781                                    nims, nime, njms, njme, nkms, nkme, &
1782                                    nips, nipe, njps, njpe, nkps, nkpe )
1784 !     print * , 'MOUSTA_VER  ',parent_grid%e_vert , kde_c , kde_n
1785 !        print *, nids , nide , njds , njde , nkds , nkde
1786 !        print *, nims , nime , njms , njme , nkms , nkme
1787 !        print *, nips , nipe , njps , njpe , nkps , nkpe
1791          allocate (alt_w_c(kde_c))
1792          allocate (alt_u_c(kde_c+1))
1793          allocate (pro_w_c(kde_c))
1794          allocate (pro_u_c(kde_c+1))
1796          allocate (alt_w_n(kde_n))
1797          allocate (alt_u_n(kde_n+1))
1798          allocate (pro_w_n(kde_n))
1799          allocate (pro_u_n(kde_n+1))
1801 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11111
1802          p_top_m = parent_grid%p_top
1803          p_surf_m = 1.e5
1804          mu_m = p_surf_m - p_top_m
1805 !        print * , 'p_top_m', p_top_m
1806 !    parent
1807          do  k = 1,kde_c
1808          pre_c = mu_m * c3f(k) + c4f(k) + p_top_m
1809          alt_w_c(k) =  -hsca_m * alog(pre_c/p_surf_m)
1810          enddo
1811          do  k = 1,kde_c-1
1812          pre_c = mu_m * c3h(k) + c4h(k) + p_top_m
1813          alt_u_c(k+1) =  -hsca_m * alog(pre_c/p_surf_m)
1814          enddo
1815          alt_u_c(1) =  alt_w_c(1)
1816          alt_u_c(kde_c+1) =  alt_w_c(kde_c)
1817 !    nest
1818          do  k = 1,kde_n
1819          pre_n = mu_m * parent_grid%c3f(k) + parent_grid%c4f(k) + p_top_m
1820          alt_w_n(k) =  -hsca_m * alog(pre_n/p_surf_m)
1821          enddo
1822          do  k = 1,kde_n-1
1823          pre_n = mu_m * parent_grid%c3h(k) + parent_grid%c4h(k) + p_top_m
1824          alt_u_n(k+1) =  -hsca_m * alog(pre_n/p_surf_m)
1825          enddo
1826          alt_u_n(1) =  alt_w_n(1)
1827          alt_u_n(kde_n+1) =  alt_w_n(kde_n)
1828 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1830 IF ( SIZE( parent_grid%u_2, 1 ) * SIZE( parent_grid%u_2, 3 ) .GT. 1 ) THEN
1831 do j = njms , njme
1832 do i = nims , nime
1833     
1834     do k = 1,kde_c-1
1835     pro_u_c(k+1) = parent_grid%u_2(i,k,j)
1836     enddo
1837     pro_u_c(1      ) = cf1_c*parent_grid%u_2(i,1,j) &
1838                      + cf2_c*parent_grid%u_2(i,2,j) &
1839                      + cf3_c*parent_grid%u_2(i,3,j)
1841     pro_u_c(kde_c+1) = cfn_c *parent_grid%u_2(i,kde_c-1,j) &
1842                      + cfn1_c*parent_grid%u_2(i,kde_c-2,j)
1844     call inter(pro_u_c,alt_u_c,kde_c+1,pro_u_n,alt_u_n,kde_n+1)
1846     do k = 1,kde_n-1
1847     parent_grid%u_2(i,k,j) = pro_u_n(k+1)
1848     enddo
1850 enddo
1851 enddo
1852 ENDIF
1854 IF ( SIZE( parent_grid%v_2, 1 ) * SIZE( parent_grid%v_2, 3 ) .GT. 1 ) THEN
1855 do j = njms , njme
1856 do i = nims , nime
1858     do k = 1,kde_c-1
1859     pro_u_c(k+1) = parent_grid%v_2(i,k,j)
1860     enddo
1861     pro_u_c(1      ) = cf1_c*parent_grid%v_2(i,1,j) &
1862                      + cf2_c*parent_grid%v_2(i,2,j) &
1863                      + cf3_c*parent_grid%v_2(i,3,j)
1865     pro_u_c(kde_c+1) = cfn_c *parent_grid%v_2(i,kde_c-1,j) &
1866                      + cfn1_c*parent_grid%v_2(i,kde_c-2,j)
1868     call inter(pro_u_c,alt_u_c,kde_c+1,pro_u_n,alt_u_n,kde_n+1)
1870     do k = 1,kde_n-1
1871     parent_grid%v_2(i,k,j) = pro_u_n(k+1)
1872     enddo
1874 enddo
1875 enddo
1876 ENDIF
1878 IF ( SIZE( parent_grid%w_2, 1 ) * SIZE( parent_grid%w_2, 3 ) .GT. 1 ) THEN
1879 do j = njms , njme
1880 do i = nims , nime
1882     do k = 1,kde_c
1883     pro_w_c(k) = parent_grid%w_2(i,k,j)
1884     enddo
1885     call inter(pro_w_c,alt_w_c,kde_c,pro_w_n,alt_w_n,kde_n)
1887     do k = 1,kde_n
1888     parent_grid%w_2(i,k,j) = pro_w_n(k)
1889     enddo
1890 enddo
1891 enddo
1892 ENDIF
1894 IF ( SIZE( parent_grid%t_2, 1 ) * SIZE( parent_grid%t_2, 3 ) .GT. 1 ) THEN
1895 do j = njms , njme
1896 do i = nims , nime
1898     do k = 1,kde_c-1
1899     pro_u_c(k+1) = parent_grid%t_2(i,k,j)
1900     enddo
1901     pro_u_c(1      ) = cf1_c*parent_grid%t_2(i,1,j) &
1902                      + cf2_c*parent_grid%t_2(i,2,j) &
1903                      + cf3_c*parent_grid%t_2(i,3,j)
1905     pro_u_c(kde_c+1) = cfn_c *parent_grid%t_2(i,kde_c-1,j) &
1906                      + cfn1_c*parent_grid%t_2(i,kde_c-2,j)
1908     call inter(pro_u_c,alt_u_c,kde_c+1,pro_u_n,alt_u_n,kde_n+1)
1910     do k = 1,kde_n-1
1911     parent_grid%t_2(i,k,j) = pro_u_n(k+1)
1912     enddo
1914 enddo
1915 enddo
1916 ENDIF
1918 DO itrace = PARAM_FIRST_SCALAR, num_moist
1919 IF ( SIZE( parent_grid%moist, 1 ) * SIZE( parent_grid%moist, 3 ) .GT. 1 ) THEN
1920 do j = njms , njme
1921 do i = nims , nime
1923     do k = 1,kde_c-1
1924     pro_u_c(k+1) = parent_grid%moist(i,k,j,itrace)
1925     enddo
1926     pro_u_c(1      ) = cf1_c*parent_grid%moist(i,1,j,itrace) &
1927                      + cf2_c*parent_grid%moist(i,2,j,itrace) &
1928                      + cf3_c*parent_grid%moist(i,3,j,itrace)
1930     pro_u_c(kde_c+1) = cfn_c *parent_grid%moist(i,kde_c-1,j,itrace) &
1931                      + cfn1_c*parent_grid%moist(i,kde_c-2,j,itrace)
1933     call inter(pro_u_c,alt_u_c,kde_c+1,pro_u_n,alt_u_n,kde_n+1)
1935     do k = 1,kde_n-1
1936     parent_grid%moist(i,k,j,itrace) = pro_u_n(k+1)
1937     enddo
1939 enddo
1940 enddo
1941 ENDIF
1942 ENDDO
1944 DO itrace = PARAM_FIRST_SCALAR, num_dfi_moist
1945 IF ( SIZE( parent_grid%dfi_moist, 1 ) * SIZE( parent_grid%dfi_moist, 3 ) .GT. 1 ) THEN
1946 do j = njms , njme
1947 do i = nims , nime
1949     do k = 1,kde_c-1
1950     pro_u_c(k+1) = parent_grid%dfi_moist(i,k,j,itrace)
1951     enddo
1952     pro_u_c(1      ) = cf1_c*parent_grid%dfi_moist(i,1,j,itrace) &
1953                      + cf2_c*parent_grid%dfi_moist(i,2,j,itrace) &
1954                      + cf3_c*parent_grid%dfi_moist(i,3,j,itrace)
1956     pro_u_c(kde_c+1) = cfn_c *parent_grid%dfi_moist(i,kde_c-1,j,itrace) &
1957                      + cfn1_c*parent_grid%dfi_moist(i,kde_c-2,j,itrace)
1959     call inter(pro_u_c,alt_u_c,kde_c+1,pro_u_n,alt_u_n,kde_n+1)
1961     do k = 1,kde_n-1
1962     parent_grid%dfi_moist(i,k,j,itrace) = pro_u_n(k+1)
1963     enddo
1965 enddo
1966 enddo
1967 ENDIF
1968 ENDDO
1970 DO itrace = PARAM_FIRST_SCALAR, num_scalar
1971 IF ( SIZE( parent_grid%scalar, 1 ) * SIZE( parent_grid%scalar, 3 ) .GT. 1 ) THEN
1972 do j = njms , njme
1973 do i = nims , nime
1975     do k = 1,kde_c-1
1976     pro_u_c(k+1) = parent_grid%scalar(i,k,j,itrace)
1977     enddo
1978     pro_u_c(1      ) = cf1_c*parent_grid%scalar(i,1,j,itrace) &
1979                      + cf2_c*parent_grid%scalar(i,2,j,itrace) &
1980                      + cf3_c*parent_grid%scalar(i,3,j,itrace)
1982     pro_u_c(kde_c+1) = cfn_c *parent_grid%scalar(i,kde_c-1,j,itrace) &
1983                      + cfn1_c*parent_grid%scalar(i,kde_c-2,j,itrace)
1985     call inter(pro_u_c,alt_u_c,kde_c+1,pro_u_n,alt_u_n,kde_n+1)
1987     do k = 1,kde_n-1
1988     parent_grid%scalar(i,k,j,itrace) = pro_u_n(k+1)
1989     enddo
1991 enddo
1992 enddo
1993 ENDIF
1994 ENDDO
1996 DO itrace = PARAM_FIRST_SCALAR, num_dfi_scalar
1997 IF ( SIZE( parent_grid%dfi_scalar, 1 ) * SIZE( parent_grid%dfi_scalar, 3 ) .GT. 1 ) THEN
1998 do j = njms , njme
1999 do i = nims , nime
2001     do k = 1,kde_c-1
2002     pro_u_c(k+1) = parent_grid%dfi_scalar(i,k,j,itrace)
2003     enddo
2004     pro_u_c(1      ) = cf1_c*parent_grid%dfi_scalar(i,1,j,itrace) &
2005                      + cf2_c*parent_grid%dfi_scalar(i,2,j,itrace) &
2006                      + cf3_c*parent_grid%dfi_scalar(i,3,j,itrace)
2008     pro_u_c(kde_c+1) = cfn_c *parent_grid%dfi_scalar(i,kde_c-1,j,itrace) &
2009                      + cfn1_c*parent_grid%dfi_scalar(i,kde_c-2,j,itrace)
2011     call inter(pro_u_c,alt_u_c,kde_c+1,pro_u_n,alt_u_n,kde_n+1)
2013     do k = 1,kde_n-1
2014     parent_grid%dfi_scalar(i,k,j,itrace) = pro_u_n(k+1)
2015     enddo
2017 enddo
2018 enddo
2019 ENDIF
2020 ENDDO
2024 IF ( SIZE( parent_grid%f_ice_phy, 1 ) * SIZE( parent_grid%f_ice_phy, 3 ) .GT. 1 ) THEN
2025 do j = njms , njme
2026 do i = nims , nime
2028     do k = 1,kde_c-1
2029     pro_u_c(k+1) = parent_grid%f_ice_phy(i,k,j)
2030     enddo
2031     pro_u_c(1      ) = cf1_c*parent_grid%f_ice_phy(i,1,j) &
2032                      + cf2_c*parent_grid%f_ice_phy(i,2,j) &
2033                      + cf3_c*parent_grid%f_ice_phy(i,3,j)
2035     pro_u_c(kde_c+1) = cfn_c *parent_grid%f_ice_phy(i,kde_c-1,j) &
2036                      + cfn1_c*parent_grid%f_ice_phy(i,kde_c-2,j)
2038     call inter(pro_u_c,alt_u_c,kde_c+1,pro_u_n,alt_u_n,kde_n+1)
2040     do k = 1,kde_n-1
2041     parent_grid%f_ice_phy(i,k,j) = pro_u_n(k+1)
2042     enddo
2044 enddo
2045 enddo
2046 ENDIF
2048 IF ( SIZE( parent_grid%f_rain_phy, 1 ) * SIZE( parent_grid%f_rain_phy, 3 ) .GT. 1 ) THEN
2049 do j = njms , njme
2050 do i = nims , nime
2052     do k = 1,kde_c-1
2053     pro_u_c(k+1) = parent_grid%f_rain_phy(i,k,j)
2054     enddo
2055     pro_u_c(1      ) = cf1_c*parent_grid%f_rain_phy(i,1,j) &
2056                      + cf2_c*parent_grid%f_rain_phy(i,2,j) &
2057                      + cf3_c*parent_grid%f_rain_phy(i,3,j)
2059     pro_u_c(kde_c+1) = cfn_c *parent_grid%f_rain_phy(i,kde_c-1,j) &
2060                      + cfn1_c*parent_grid%f_rain_phy(i,kde_c-2,j)
2062     call inter(pro_u_c,alt_u_c,kde_c+1,pro_u_n,alt_u_n,kde_n+1)
2064     do k = 1,kde_n-1
2065     parent_grid%f_rain_phy(i,k,j) = pro_u_n(k+1)
2066     enddo
2068 enddo
2069 enddo
2070 ENDIF
2073 IF ( SIZE( parent_grid%f_rimef_phy, 1 ) * SIZE( parent_grid%f_rimef_phy, 3 ) .GT. 1 ) THEN
2074 do j = njms , njme
2075 do i = nims , nime
2077     do k = 1,kde_c-1
2078     pro_u_c(k+1) = parent_grid%f_rimef_phy(i,k,j)
2079     enddo
2080     pro_u_c(1      ) = cf1_c*parent_grid%f_rimef_phy(i,1,j) &
2081                      + cf2_c*parent_grid%f_rimef_phy(i,2,j) &
2082                      + cf3_c*parent_grid%f_rimef_phy(i,3,j)
2084     pro_u_c(kde_c+1) = cfn_c *parent_grid%f_rimef_phy(i,kde_c-1,j) &
2085                      + cfn1_c*parent_grid%f_rimef_phy(i,kde_c-2,j)
2087     call inter(pro_u_c,alt_u_c,kde_c+1,pro_u_n,alt_u_n,kde_n+1)
2089     do k = 1,kde_n-1
2090     parent_grid%f_rimef_phy(i,k,j) = pro_u_n(k+1)
2091     enddo
2093 enddo
2094 enddo
2095 ENDIF
2097 IF ( SIZE( parent_grid%h_diabatic, 1 ) * SIZE( parent_grid%h_diabatic, 3 ) .GT. 1 ) THEN
2098 do j = njms , njme
2099 do i = nims , nime
2101     do k = 1,kde_c-1
2102     pro_u_c(k+1) = parent_grid%h_diabatic(i,k,j)
2103     enddo
2104     pro_u_c(1      ) = cf1_c*parent_grid%h_diabatic(i,1,j) &
2105                      + cf2_c*parent_grid%h_diabatic(i,2,j) &
2106                      + cf3_c*parent_grid%h_diabatic(i,3,j)
2108     pro_u_c(kde_c+1) = cfn_c *parent_grid%h_diabatic(i,kde_c-1,j) &
2109                      + cfn1_c*parent_grid%h_diabatic(i,kde_c-2,j)
2111     call inter(pro_u_c,alt_u_c,kde_c+1,pro_u_n,alt_u_n,kde_n+1)
2113     do k = 1,kde_n-1
2114     parent_grid%h_diabatic(i,k,j) = pro_u_n(k+1)
2115     enddo
2117 enddo
2118 enddo
2119 ENDIF
2121 IF ( SIZE( parent_grid%rthraten, 1 ) * SIZE( parent_grid%rthraten, 3 ) .GT. 1 ) THEN
2122 do j = njms , njme
2123 do i = nims , nime
2125     do k = 1,kde_c-1
2126     pro_u_c(k+1) =  parent_grid%rthraten(i,k,j)
2127     enddo
2128     pro_u_c(1      ) = cf1_c*parent_grid%rthraten(i,1,j) &
2129                      + cf2_c*parent_grid%rthraten(i,2,j) &
2130                      + cf3_c*parent_grid%rthraten(i,3,j)
2132     pro_u_c(kde_c+1) = cfn_c *parent_grid%rthraten(i,kde_c-1,j) &
2133                      + cfn1_c*parent_grid%rthraten(i,kde_c-2,j)
2135     call inter(pro_u_c,alt_u_c,kde_c+1,pro_u_n,alt_u_n,kde_n+1)
2137     do k = 1,kde_n-1
2138     parent_grid%rthraten(i,k,j) = pro_u_n(k+1)
2139     enddo
2141 enddo
2142 enddo
2143 ENDIF
2153 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2154 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2155          deallocate (alt_w_c)
2156          deallocate (alt_u_c)
2157          deallocate (pro_w_c)
2158          deallocate (pro_u_c)
2160          deallocate (alt_w_n)
2161          deallocate (alt_u_n)
2162          deallocate (pro_w_n)
2163          deallocate (pro_u_n)
2166       END SUBROUTINE vertical_interp
2168 !!!!!!!!!!!!!!!!!!!!!!!!11
2169     SUBROUTINE inter(pro_c,alt_c,kde_c,pro_n,alt_n,kde_n)
2171   IMPLICIT NONE
2172    INTEGER , INTENT(IN) :: kde_c,kde_n
2173    REAL , DIMENSION(kde_c) , INTENT(IN ) :: pro_c,alt_c
2174    REAL , DIMENSION(kde_n) , INTENT(IN ) :: alt_n
2175    REAL , DIMENSION(kde_n) , INTENT(OUT) :: pro_n
2177       real ,dimension(kde_c) :: a,b,c,d
2178       real :: p
2179       integer :: i,j
2182       call coeff_mon(alt_c,pro_c,a,b,c,d,kde_c)
2184        do i = 1,kde_n-1
2186           do j=1,kde_c-1
2188                if ( (alt_n(i) .ge. alt_c(j)).and.(alt_n(i) .lt. alt_c(j+1)) ) then
2189                 p =  alt_n(i)-alt_c(j)
2190                 pro_n(i) = p*( p*(a(j)*p+b(j))+c(j)) + d(j)
2191                 goto 20
2192                endif
2193            enddo
2194 20             continue
2195        enddo
2197        pro_n(kde_n) = pro_c(kde_c)
2200    END SUBROUTINE inter
2202 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
2204      subroutine  coeff_mon(x,y,a,b,c,d,n)
2206       implicit none
2208       integer :: n
2209       real ,dimension(n) :: x,y,a,b,c,d
2210       real ,dimension(n) :: h,s,p,yp
2212        integer :: i
2215       do i=1,n-1
2216       h(i) = (x(i+1)-x(i))
2217       s(i) = (y(i+1)-y(i)) / h(i)
2218       enddo
2220       do i=2,n-1
2221       p(i) = (s(i-1)*h(i)+s(i)*h(i-1)) / (h(i-1)+h(i))
2222       enddo
2224       p(1) = s(1)
2225       p(n) = s(n-1)
2227       do i=1,n
2228       yp(i) = p(i)
2229       enddo
2230 !!!!!!!!!!!!!!!!!!!!!
2232       do i=2,n-1
2233       yp(i) = (sign(1.,s(i-1))+sign(1.,s(i)))* min( abs(s(i-1)),abs(s(i)),0.5*abs(p(i)))
2234       enddo
2236       do i = 1,n-1
2237       a(i) = (yp(i)+yp(i+1)-2.*s(i))/(h(i)*h(i))
2238       b(i) = (3.*s(i)-2.*yp(i)-yp(i+1))/h(i)
2239       c(i) = yp(i)
2240       d(i) = y(i)
2241       enddo
2243       end  subroutine  coeff_mon