source: lmdz_wrf/trunk/WRFV3/dyn_em/start_em.F @ 1531

Last change on this file since 1531 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 56.2 KB
Line 
1!-------------------------------------------------------------------
2
3   SUBROUTINE start_domain_em ( grid, allowed_to_read &
4! Actual arguments generated from Registry
5# include "dummy_new_args.inc"
6!
7)
8
9   USE module_domain, ONLY : domain, wrfu_timeinterval, get_ijk_from_grid, &
10        domain_setgmtetc
11   USE module_state_description
12   USE module_model_constants
13   USE module_bc, ONLY : boundary_condition_check, set_physical_bc2d
14   USE module_bc_em
15   USE module_configure, ONLY : grid_config_rec_type
16   USE module_tiles, ONLY : set_tiles
17#ifdef DM_PARALLEL
18   USE module_dm, ONLY : wrf_dm_min_real, wrf_dm_maxval, ntasks_x, ntasks_y, &
19        local_communicator_periodic, local_communicator, mytask, ntasks
20#else
21   USE module_dm, ONLY : wrf_dm_min_real
22#endif
23   USE module_comm_dm
24
25   USE module_physics_init
26   USE module_fr_sfire_driver_wrf, ONLY : sfire_driver_em_init
27   USE module_stoch, ONLY : SETUP_STOCH,rand_seed, update_stoch
28#ifdef WRF_CHEM
29   USE module_aerosols_sorgam, ONLY: sum_pm_sorgam
30   USE module_gocart_aerosols, ONLY: sum_pm_gocart
31   USE module_mosaic_driver, ONLY: sum_pm_mosaic
32   USE module_input_tracer, ONLY: initialize_tracer
33#endif
34
35!!debug
36!USE module_compute_geop
37
38   USE module_model_constants
39   USE module_avgflx_em, ONLY : zero_avgflx
40
41   IMPLICIT NONE
42   !  Input data.
43   TYPE (domain)          :: grid
44
45   LOGICAL , INTENT(IN)   :: allowed_to_read
46
47   !  Definitions of dummy arguments to this routine (generated from Registry).
48# include "dummy_new_decl.inc"
49
50   !  Structure that contains run-time configuration (namelist) data for domain
51   TYPE (grid_config_rec_type)              :: config_flags
52
53   !  Local data
54   INTEGER                             ::                       &
55                                  ids, ide, jds, jde, kds, kde, &
56                                  ims, ime, jms, jme, kms, kme, &
57                                  ips, ipe, jps, jpe, kps, kpe, &
58                                  its, ite, jts, jte, kts, kte, &
59                                  ij,i,j,k,ii,jj,kk,loop,error,l
60
61   INTEGER ::              imsx, imex, jmsx, jmex, kmsx, kmex,    &
62                           ipsx, ipex, jpsx, jpex, kpsx, kpex,    &
63                           imsy, imey, jmsy, jmey, kmsy, kmey,    &
64                           ipsy, ipey, jpsy, jpey, kpsy, kpey
65
66   INTEGER     :: i_m
67
68   REAL        :: p00, t00, a, tiso, p_surf, pd_surf, temp
69#ifdef WRF_CHEM
70   REAL        RGASUNIV ! universal gas constant [ J/mol-K ]
71   PARAMETER ( RGASUNIV = 8.314510 )
72      REAL,DIMENSION(grid%sm31:grid%em31,grid%sm32:grid%em32,grid%sm33:grid%em33) :: &
73                      z_at_w,convfac
74   REAL :: tempfac
75#endif
76
77   REAL :: qvf1, qvf2, qvf
78   REAL :: MPDT
79   REAL :: spongeweight
80   LOGICAL :: first_trip_for_this_domain, start_of_simulation, fill_w_flag
81   LOGICAL, EXTERNAL :: wrf_dm_on_monitor
82
83#ifndef WRF_CHEM
84      REAL,ALLOCATABLE,DIMENSION(:,:,:) :: cldfra_old
85#endif
86
87   REAL :: lat1 , lat2 , lat3 , lat4
88   REAL :: lon1 , lon2 , lon3 , lon4
89   INTEGER :: num_points_lat_lon , iloc , jloc
90   CHARACTER (LEN=132) :: message
91   TYPE(WRFU_TimeInterval) :: stepTime
92   REAL, DIMENSION(:,:), ALLOCATABLE :: clat_glob
93   logical :: f_flux  ! flag for computing averaged fluxes in cu_gd
94
95   INTEGER :: idex, jdex
96   INTEGER :: im1,ip1,jm1,jp1
97   REAL :: hx,hy,pi
98
99   CALL get_ijk_from_grid ( grid ,                              &
100                           ids, ide, jds, jde, kds, kde,        &
101                           ims, ime, jms, jme, kms, kme,        &
102                           ips, ipe, jps, jpe, kps, kpe,        &
103                           imsx, imex, jmsx, jmex, kmsx, kmex,  &
104                           ipsx, ipex, jpsx, jpex, kpsx, kpex,  &
105                           imsy, imey, jmsy, jmey, kmsy, kmey,  &
106                           ipsy, ipey, jpsy, jpey, kpsy, kpey )
107
108   kts = kps ; kte = kpe     ! note that tile is entire patch
109   its = ips ; ite = ipe     ! note that tile is entire patch
110   jts = jps ; jte = jpe    ! note that tile is entire patch
111#ifndef WRF_CHEM
112         ALLOCATE(CLDFRA_OLD(IMS:IME,KMS:KME,JMS:JME),STAT=I)  ; CLDFRA_OLD = 0.
113#endif
114   CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
115
116   IF ( ( MOD (ide-ids,config_flags%parent_grid_ratio) .NE. 0 ) .OR. &
117        ( MOD (jde-jds,config_flags%parent_grid_ratio) .NE. 0 ) ) THEN
118      WRITE(message, FMT='("Nested dimensions are illegal for domain ",I2,":  Both &
119         &MOD(",I4,"-",I1,",",I2,") and MOD(",I4,"-",I1,",",I2,") must = 0" )') &
120         grid%id,ide,ids,config_flags%parent_grid_ratio,jde,jds,config_flags%parent_grid_ratio
121      CALL wrf_error_fatal ( message )
122   END IF
123
124!   --- SETUP AND INITIALIZE STOCHASTIC KINETIC ENERGY BACKSCATTER SCHEME ---
125
126   IF ( (grid%id.eq.1) .AND. ( .NOT. grid%did_stoch ) ) THEN
127   grid%did_stoch = .TRUE.
128   stoch_force_select: SELECT CASE(config_flags%stoch_force_opt)
129
130   CASE (STOCH_BACKSCATTER)
131
132     IF ( wrf_dm_on_monitor () ) THEN
133       CALL rand_seed ( config_flags, grid%SEED1, grid%SEED2, grid%NENS )
134     ENDIF
135#ifdef DM_PARALLEL
136     CALL wrf_dm_bcast_bytes ( grid%SEED1, IWORDSIZE )
137     CALL wrf_dm_bcast_bytes ( grid%SEED2, IWORDSIZE )
138#endif
139
140     call SETUP_STOCH(grid%VERTSTRUCC,grid%VERTSTRUCS,                   &
141                      grid%SPT_AMP,grid%SPSTREAM_AMP,                    &
142                      grid%stoch_vertstruc_opt,                          &
143                      grid%SEED1,grid%SEED2,grid%time_step,              &
144                      grid%DX,grid%DY,                                   &
145                      grid%TOT_BACKSCAT_PSI,grid%TOT_BACKSCAT_T,         &
146                      ids, ide, jds, jde, kds, kde,                      &
147                      ims, ime, jms, jme, kms, kme,                      &
148                      its, ite, jts, jte, kts, kte                       )
149   
150    END SELECT stoch_force_select
151    END IF
152!   --- END SETUP STOCHASTIC KINETIC ENERGY BACKSCATTER SCHEME ----------
153
154   IF ( config_flags%polar ) THEN
155!write(0,*)__FILE__,__LINE__,' clat ',ips,ipe,jps,jpe
156!do j = jps,jpe
157!write(0,*)__FILE__,__LINE__,' clat ',ids,j,grid%clat(ips,j)
158!enddo
159
160#ifdef DM_PARALLEL
161! WARNING:  this might present scaling issues on very large numbers of processors
162     ALLOCATE( clat_glob(ids:ide,jds:jde) )
163
164     CALL wrf_patch_to_global_real ( grid%clat, clat_glob, grid%domdesc, 'xy', 'xy', &
165                                     ids, ide, jds, jde, 1, 1, &
166                                     ims, ime, jms, jme, 1, 1, &
167                                     its, ite, jts, jte, 1, 1 )
168
169     CALL wrf_dm_bcast_real ( clat_glob , (ide-ids+1)*(jde-jds+1) )
170
171     grid%clat_xxx(ipsx:ipex,jpsx:jpex) = clat_glob(ipsx:ipex,jpsx:jpex)
172
173     DEALLOCATE( clat_glob )
174#endif
175   ENDIF
176
177! here we check to see if the boundary conditions are set properly
178
179   CALL boundary_condition_check( config_flags, bdyzone, error, grid%id )
180
181!kludge - need to stop CG from resetting precip and phys tendencies to zero
182!         when we are in here due to a nest being spawned, we want to still
183!         recompute the base state, but that is about it
184   !  This is temporary and will need to be changed when grid%itimestep is removed.
185
186   IF ( grid%itimestep .EQ. 0 ) THEN
187      first_trip_for_this_domain = .TRUE.
188   ELSE
189      first_trip_for_this_domain = .FALSE.
190   END IF
191
192   IF ( .not. ( config_flags%restart .or. grid%moved ) ) THEN
193       grid%itimestep=0
194   ENDIF
195
196   IF ( config_flags%restart .or. grid%moved ) THEN
197       first_trip_for_this_domain = .TRUE.
198   ENDIF
199
200! wig: Add a combined exponential+linear weight on the mother boundaries
201!      following code changes by Ruby Leung. For the nested grid, there
202!      appears to be some problems when a sponge is used. The points where
203!      processors meet have problematic values.
204
205   CALL lbc_fcx_gcx ( grid%fcx , grid%gcx , grid%spec_bdy_width , &
206                      grid%spec_zone , grid%relax_zone , grid%dt , config_flags%spec_exp , &
207                      config_flags%specified , config_flags%nested )
208
209   IF ( config_flags%nested ) THEN
210     grid%dtbc = 0.
211   ENDIF
212
213   IF ( ( grid%id .NE. 1 ) .AND. ( .NOT. config_flags%input_from_file ) ) THEN
214
215      !  Every time a domain starts or every time a domain moves, this routine is called.  We want
216      !  the center (middle) lat/lon of the grid for the metacode.  The lat/lon values are
217      !  defined at mass points.  Depending on the even/odd points in the SN and WE directions,
218      !  we end up with the middle point as either 1 point or an average of either 2 or 4 points.
219      !  Add to this, the need to make sure that we are on the correct patch to retrieve the
220      !  value of the lat/lon, AND that the lat/lons (for an average) may not all be on the same
221      !  patch.  Once we find the correct value for lat lon, we need to keep it around on all patches,
222      !  which is where the wrf_dm_min_real calls come in.
223      !  If this is the most coarse domain, we do not go in here.  Also, if there is an input file
224      !  (which has the right values for the middle lat/lon) we do not go in this IF test.
225
226      IF      ( ( MOD(ide,2) .EQ. 0 ) .AND. ( MOD(jde,2) .EQ. 0 ) ) THEN
227         num_points_lat_lon = 1
228         iloc = ide/2
229         jloc = jde/2
230         IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
231              ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
232            lat1 = grid%xlat (iloc,jloc)
233            lon1 = grid%xlong(iloc,jloc)
234         ELSE
235            lat1 = 99999.
236            lon1 = 99999.
237         END IF
238         lat1 = wrf_dm_min_real ( lat1 )
239         lon1 = wrf_dm_min_real ( lon1 )
240         CALL nl_set_cen_lat ( grid%id , lat1 )
241         CALL nl_set_cen_lon ( grid%id , lon1 )
242      ELSE IF ( ( MOD(ide,2) .NE. 0 ) .AND. ( MOD(jde,2) .EQ. 0 ) ) THEN
243         num_points_lat_lon = 2
244         iloc = (ide-1)/2
245         jloc =  jde   /2
246         IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
247              ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
248            lat1 = grid%xlat (iloc,jloc)
249            lon1 = grid%xlong(iloc,jloc)
250         ELSE
251            lat1 = 99999.
252            lon1 = 99999.
253         END IF
254         lat1 = wrf_dm_min_real ( lat1 )
255         lon1 = wrf_dm_min_real ( lon1 )
256
257         iloc = (ide+1)/2
258         jloc =  jde   /2
259         IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
260              ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
261            lat2 = grid%xlat (iloc,jloc)
262            lon2 = grid%xlong(iloc,jloc)
263         ELSE
264            lat2 = 99999.
265            lon2 = 99999.
266         END IF
267         lat2 = wrf_dm_min_real ( lat2 )
268         lon2 = wrf_dm_min_real ( lon2 )
269
270         CALL nl_set_cen_lat ( grid%id , ( lat1 + lat2 ) * 0.50 )
271         CALL nl_set_cen_lon ( grid%id , ( lon1 + lon2 ) * 0.50 )
272      ELSE IF ( ( MOD(ide,2) .EQ. 0 ) .AND. ( MOD(jde,2) .NE. 0 ) ) THEN
273         num_points_lat_lon = 2
274         iloc =  ide   /2
275         jloc = (jde-1)/2
276         IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
277              ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
278            lat1 = grid%xlat (iloc,jloc)
279            lon1 = grid%xlong(iloc,jloc)
280         ELSE
281            lat1 = 99999.
282            lon1 = 99999.
283         END IF
284         lat1 = wrf_dm_min_real ( lat1 )
285         lon1 = wrf_dm_min_real ( lon1 )
286
287         iloc =  ide   /2
288         jloc = (jde+1)/2
289         IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
290              ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
291            lat2 = grid%xlat (iloc,jloc)
292            lon2 = grid%xlong(iloc,jloc)
293         ELSE
294            lat2 = 99999.
295            lon2 = 99999.
296         END IF
297         lat2 = wrf_dm_min_real ( lat2 )
298         lon2 = wrf_dm_min_real ( lon2 )
299
300         CALL nl_set_cen_lat ( grid%id , ( lat1 + lat2 ) * 0.50 )
301         CALL nl_set_cen_lon ( grid%id , ( lon1 + lon2 ) * 0.50 )
302      ELSE IF ( ( MOD(ide,2) .NE. 0 ) .AND. ( MOD(jde,2) .NE. 0 ) ) THEN
303         num_points_lat_lon = 4
304         iloc = (ide-1)/2
305         jloc = (jde-1)/2
306         IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
307              ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
308            lat1 = grid%xlat (iloc,jloc)
309            lon1 = grid%xlong(iloc,jloc)
310         ELSE
311            lat1 = 99999.
312            lon1 = 99999.
313         END IF
314         lat1 = wrf_dm_min_real ( lat1 )
315         lon1 = wrf_dm_min_real ( lon1 )
316
317         iloc = (ide+1)/2
318         jloc = (jde-1)/2
319         IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
320              ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
321            lat2 = grid%xlat (iloc,jloc)
322            lon2 = grid%xlong(iloc,jloc)
323         ELSE
324            lat2 = 99999.
325            lon2 = 99999.
326         END IF
327         lat2 = wrf_dm_min_real ( lat2 )
328         lon2 = wrf_dm_min_real ( lon2 )
329
330         iloc = (ide-1)/2
331         jloc = (jde+1)/2
332         IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
333              ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
334            lat3 = grid%xlat (iloc,jloc)
335            lon3 = grid%xlong(iloc,jloc)
336         ELSE
337            lat3 = 99999.
338            lon3 = 99999.
339         END IF
340         lat3 = wrf_dm_min_real ( lat3 )
341         lon3 = wrf_dm_min_real ( lon3 )
342
343         iloc = (ide+1)/2
344         jloc = (jde+1)/2
345         IF ( ( ips .LE. iloc ) .AND. ( ipe .GE. iloc ) .AND. &
346              ( jps .LE. jloc ) .AND. ( jpe .GE. jloc ) ) THEN
347            lat4 = grid%xlat (iloc,jloc)
348            lon4 = grid%xlong(iloc,jloc)
349         ELSE
350            lat4 = 99999.
351            lon4 = 99999.
352         END IF
353         lat4 = wrf_dm_min_real ( lat4 )
354         lon4 = wrf_dm_min_real ( lon4 )
355
356         CALL nl_set_cen_lat ( grid%id , ( lat1 + lat2 + lat3 + lat4 ) * 0.25 )
357         CALL nl_set_cen_lon ( grid%id , ( lon1 + lon2 + lon3 + lon4 ) * 0.25 )
358      END IF
359   END IF
360
361   IF ( .NOT. config_flags%restart .AND. &
362        (( config_flags%input_from_hires ) .OR. ( config_flags%input_from_file ))) THEN
363
364      IF ( config_flags%map_proj .EQ. 0 ) THEN
365         CALL wrf_error_fatal ( 'start_domain: Idealized case cannot have a separate nested input file' )
366      END IF
367
368      IF ( config_flags%use_baseparam_fr_nml ) then
369      CALL nl_get_base_pres  ( 1 , p00 )
370      CALL nl_get_base_temp  ( 1 , t00 )
371      CALL nl_get_base_lapse ( 1 , a   )
372      CALL nl_get_iso_temp   ( 1 , tiso )
373
374      ELSE
375      ! get these constants from model data
376
377      t00  = grid%t00
378      p00  = grid%p00
379      a    = grid%tlp
380      tiso = grid%tiso
381
382      IF (t00 .LT. 100. .or. p00 .LT. 10000.) THEN
383      WRITE(wrf_err_message,*)&
384      'start_em: did not find base state parameters in wrfinput. Add use_baseparam_fr_nml = .t. in &dynamics and rerun'
385      CALL wrf_error_fatal(TRIM(wrf_err_message))
386      ENDIF
387
388      ENDIF
389
390      !  Base state potential temperature and inverse density (alpha = 1/rho) from
391      !  the half eta levels and the base-profile surface pressure.  Compute 1/rho
392      !  from equation of state.  The potential temperature is a perturbation from t0.
393
394      DO j = jts, MIN(jte,jde-1)
395         DO i = its, MIN(ite,ide-1)
396
397            !  Base state pressure is a function of eta level and terrain, only, plus
398            !  the hand full of constants: p00 (sea level pressure, Pa), t00 (sea level
399            !  temperature, K), and A (temperature difference, from 1000 mb to 300 mb, K).
400
401            p_surf = p00 * EXP ( -t00/a + ( (t00/a)**2 - 2.*g*grid%ht(i,j)/a/r_d ) **0.5 )
402
403            DO k = 1, kte-1
404               grid%pb(i,k,j) = grid%znu(k)*(p_surf - grid%p_top) + grid%p_top
405               temp = MAX ( tiso, t00 + A*LOG(grid%pb(i,k,j)/p00) )
406               grid%t_init(i,k,j) = temp*(p00/grid%pb(i,k,j))**(r_d/cp) - t0
407!              grid%t_init(i,k,j) = (t00 + A*LOG(grid%pb(i,k,j)/p00))*(p00/grid%pb(i,k,j))**(r_d/cp) - t0
408               grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
409            END DO
410
411            !  Base state mu is defined as base state surface pressure minus grid%p_top
412
413            grid%mub(i,j) = p_surf - grid%p_top
414
415            !  Integrate base geopotential, starting at terrain elevation.  This assures that
416            !  the base state is in exact hydrostatic balance with respect to the model equations.
417            !  This field is on full levels.
418
419            grid%phb(i,1,j) = grid%ht(i,j) * g
420            DO k  = 2,kte
421               grid%phb(i,k,j) = grid%phb(i,k-1,j) - grid%dnw(k-1)*grid%mub(i,j)*grid%alb(i,k-1,j)
422            END DO
423         END DO
424      END DO
425
426   ENDIF
427
428   IF(.not.config_flags%restart)THEN
429
430!  if this is for a nested domain, the defined/interpolated fields are the _2
431
432     IF ( first_trip_for_this_domain ) THEN
433
434! data that is expected to be zero must be explicitly initialized as such
435!    grid%h_diabatic = 0.
436
437       DO j = jts,min(jte,jde-1)
438       DO k = kts,kte-1
439       DO i = its, min(ite,ide-1)
440           IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
441             grid%t_1(i,k,j)=grid%t_2(i,k,j)
442           ENDIF
443       ENDDO
444       ENDDO
445       ENDDO
446
447       DO j = jts,min(jte,jde-1)
448       DO k = kts,kte
449       DO i = its, min(ite,ide-1)
450          grid%ph_1(i,k,j)=grid%ph_2(i,k,j)
451       ENDDO
452       ENDDO
453       ENDDO
454
455       DO j = jts,min(jte,jde-1)
456       DO i = its, min(ite,ide-1)
457           IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
458             grid%mu_1(i,j)=grid%mu_2(i,j)
459           ENDIF
460       ENDDO
461       ENDDO
462     END IF
463
464!  reconstitute base-state fields
465
466    IF(config_flags%max_dom .EQ. 1)THEN
467! with single domain, grid%t_init from wrfinput is OK to use
468     DO j = jts,min(jte,jde-1)
469     DO k = kts,kte-1
470     DO i = its, min(ite,ide-1)
471       IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
472         grid%pb(i,k,j) = grid%znu(k)*grid%mub(i,j)+grid%p_top
473         grid%alb(i,k,j) = (r_d/p1000mb)*(grid%t_init(i,k,j)+t0)*(grid%pb(i,k,j)/p1000mb)**cvpm
474       ENDIF
475     ENDDO
476     ENDDO
477     ENDDO
478    ELSE
479! with nests, grid%t_init generally needs recomputations (since it is not interpolated)
480     DO j = jts,min(jte,jde-1)
481     DO k = kts,kte-1
482     DO i = its, min(ite,ide-1)
483       IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
484         grid%pb(i,k,j) = grid%znu(k)*grid%mub(i,j)+grid%p_top
485         grid%alb(i,k,j) = -grid%rdnw(k)*(grid%phb(i,k+1,j)-grid%phb(i,k,j))/grid%mub(i,j)
486         grid%t_init(i,k,j) = grid%alb(i,k,j)*(p1000mb/r_d)/((grid%pb(i,k,j)/p1000mb)**cvpm) - t0
487       ENDIF
488     ENDDO
489     ENDDO
490     ENDDO
491    ENDIF
492
493! Use equations from calc_p_rho_phi to derive p and al from ph
494
495    DO j=jts,min(jte,jde-1)
496    DO k=kts,kte-1
497    DO i=its,min(ite,ide-1)
498       qvf = 1.+rvovrd*moist(i,k,j,P_QV)
499       grid%al(i,k,j)=-1./(grid%mub(i,j)+grid%mu_1(i,j))*(grid%alb(i,k,j)*grid%mu_1(i,j)  &
500                  +grid%rdnw(k)*(grid%ph_1(i,k+1,j)-grid%ph_1(i,k,j)))
501       grid%p(i,k,j)=p1000mb*( (r_d*(t0+grid%t_1(i,k,j))*qvf)/                     &
502                  (p1000mb*(grid%al(i,k,j)+grid%alb(i,k,j))) )**cpovcv  &
503                  -grid%pb(i,k,j)
504       grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(i,k,j)
505       grid%alt(i,k,j) = grid%al(i,k,j) + grid%alb(i,k,j)
506    ENDDO
507    ENDDO
508    ENDDO
509
510#if 0
511     DO j = jts,min(jte,jde-1)
512
513       k = kte-1
514       DO i = its, min(ite,ide-1)
515         IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
516           qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k,j,P_QV))
517           qvf2 = 1./(1.+qvf1)
518           qvf1 = qvf1*qvf2
519           grid%p(i,k,j) = - 0.5*(grid%mu_1(i,j)+qvf1*grid%mub(i,j))/grid%rdnw(k)/qvf2
520           grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(i,k,j)
521           qvf = 1. + rvovrd*moist(i,k,j,P_QV)
522           grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf*(((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
523           grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
524         ENDIF
525       ENDDO
526
527       DO k = kte-2, 1, -1
528       DO i = its, min(ite,ide-1)
529         IF ( grid%imask_nostag(i,j) .EQ. 1 ) THEN
530           qvf1 = 0.5*(moist(i,k,j,P_QV)+moist(i,k+1,j,P_QV))
531           qvf2 = 1./(1.+qvf1)
532           qvf1 = qvf1*qvf2
533           grid%p(i,k,j) = grid%p(i,k+1,j) - (grid%mu_1(i,j) + qvf1*grid%mub(i,j))/qvf2/grid%rdn(k+1)
534           grid%p_hyd(i,k,j) = grid%p(i,k,j) + grid%pb(i,k,j)
535           qvf = 1. + rvovrd*moist(i,k,j,P_QV)
536           grid%alt(i,k,j) = (r_d/p1000mb)*(grid%t_1(i,k,j)+t0)*qvf* &
537                        (((grid%p(i,k,j)+grid%pb(i,k,j))/p1000mb)**cvpm)
538           grid%al(i,k,j) = grid%alt(i,k,j) - grid%alb(i,k,j)
539         ENDIF
540       ENDDO
541       ENDDO
542
543     ENDDO
544#endif
545
546   ENDIF
547
548   IF ( grid%press_adj .and. ( grid%id .NE. 1 ) .AND. .NOT. ( config_flags%restart ) .AND. &
549       ( ( config_flags%input_from_hires ) .OR. ( config_flags%input_from_file ) ) ) THEN
550      DO j = jts, MIN(jte,jde-1)
551         DO i = its, MIN(ite,ide-1)
552            grid%mu_2(i,j) = grid%mu_2(i,j) + grid%al(i,1,j) / ( grid%alt(i,1,j) * grid%alb(i,1,j) ) * &
553                                    g * ( grid%ht(i,j) - grid%ht_fine(i,j) )
554         END DO
555       END DO
556       DO j = jts,min(jte,jde-1)
557       DO i = its, min(ite,ide-1)
558          grid%mu_1(i,j)=grid%mu_2(i,j)
559       ENDDO
560       ENDDO
561
562   END IF
563
564   IF ( first_trip_for_this_domain ) THEN
565
566   CALL wrf_debug ( 100 , 'start_domain_em: Before call to phy_init' )
567
568! namelist MPDT does not exist yet, so set it here
569! MPDT is the call frequency for microphysics in minutes (0 means every step)
570   MPDT = 0.
571
572! set GMT outside of phy_init because phy_init may not be called on this
573! process if, for example, it is a moving nest and if this part of the domain is not
574! being initialized (not the leading edge).
575   CALL domain_setgmtetc( grid, start_of_simulation )
576
577!-----------------------------------------------------------------------------
578! Adaptive time step: Added by T. Hutchinson, WSI  11/6/07
579!
580!
581
582   IF ( ( grid%use_adaptive_time_step ) .AND. &
583        ( ( grid%dfi_opt .EQ. DFI_NODFI ) .OR. ( grid%dfi_stage .EQ. DFI_FST ) ) ) THEN
584
585!  Calculate any variables that were not set
586
587      if (grid%starting_time_step == -1) then
588         grid%starting_time_step = NINT(6 * MIN(grid%dx,grid%dy) / 1000)
589      endif
590
591      if (grid%max_time_step == -1) then
592         grid%max_time_step = 3*grid%starting_time_step
593      endif
594
595      if (grid%min_time_step == -1) then
596         grid%min_time_step = 0.5*grid%starting_time_step
597      endif
598
599!     Set a starting timestep.
600
601      grid%dt = grid%starting_time_step
602
603!     Initialize max cfl values
604     
605      grid%last_max_vert_cfl = 0
606      grid%last_max_horiz_cfl = 0
607
608!     Check to assure that time_step_sound is to be dynamically set.
609
610      CALL nl_set_time_step_sound ( 1 , 0 )
611      grid%time_step_sound = 0
612
613      grid%max_msftx=MAXVAL(grid%msftx)
614      grid%max_msfty=MAXVAL(grid%msfty)
615#ifdef DM_PARALLEL
616      CALL wrf_dm_maxval(grid%max_msftx, idex, jdex)
617      CALL wrf_dm_maxval(grid%max_msfty, idex, jdex)
618#endif
619
620!     This first call just initializes variables.
621!     If a restart, get initialized variables from restart file
622
623      IF ( .NOT. ( config_flags%restart ) ) then
624         CALL adapt_timestep(grid, config_flags)
625      END IF
626
627
628   END IF
629
630! End of adaptive time step modifications
631!-----------------------------------------------------------------------------
632
633
634   CALL set_tiles ( grid , grid%imask_nostag, ims, ime, jms, jme, ips, ipe, jps, jpe )
635!
636! Phy init can do reads and broadcasts when initializing physics -- landuse for example. However, if
637! we're running on a reduced mesh (that is, some tasks don't have any work) we have to at least let them
638! pass through this code so the broadcasts don't hang on the other, active tasks.  Set the number of
639! tiles to a minimum of 1 and assume that the backwards patch ranges (ips=0, ipe=-1) will prevent
640! anything else from happening on the blank tasks.  JM 20080605
641!
642   if ( allowed_to_read ) grid%num_tiles = max(1,grid%num_tiles)
643!
644! Phy_init is not necessarily thread-safe; do not multi-thread this loop.
645! The tiling is to handle the fact that we may be masking off part of the computation.
646!
647   DO ij = 1, grid%num_tiles
648
649!tgs do not need physics initialization for backward DFI integration
650     IF ( grid%dfi_opt .NE. DFI_NODFI .and. grid%dfi_stage .EQ. DFI_FST) THEN     !tgs
651        grid%stepra=nint(grid%RADT*60./grid%DT)
652        grid%stepra=max(grid%stepra,1)
653        grid%stepbl=nint(grid%BLDT*60./grid%DT)
654        grid%stepbl=max(grid%stepbl,1)
655        grid%stepcu=nint(grid%CUDT*60./grid%DT)
656        grid%stepcu=max(grid%stepcu,1)
657        grid%stepfg=nint(grid%FGDT*60./grid%DT)
658        grid%stepfg=max(grid%stepfg,1)
659     ENDIF
660     
661    IF ( ( grid%dfi_opt .EQ. DFI_NODFI ) .or. &
662          ( ( grid%dfi_stage .NE. DFI_BCK ) .and. &
663          ( grid%dfi_stage .NE. DFI_STARTBCK ) ) ) THEN    !tgs, mods by tah
664
665     CALL phy_init (  grid%id , config_flags, grid%DT, grid%RESTART, grid%znw, grid%znu,   &
666                      grid%p_top, grid%tsk, grid%RADT,grid%BLDT,grid%CUDT, MPDT, &
667                      grid%rucuten, grid%rvcuten, grid%rthcuten,             &
668                      grid%rqvcuten, grid%rqrcuten, grid%rqccuten,           &
669                      grid%rqscuten, grid%rqicuten,                          &
670                      grid%rushten, grid%rvshten, grid%rthshten,             &
671                      grid%rqvshten, grid%rqrshten, grid%rqcshten,           &
672                      grid%rqsshten, grid%rqishten, grid%rqgshten,           &
673                      grid%rublten,grid%rvblten,grid%rthblten,               &
674                      grid%rqvblten,grid%rqcblten,grid%rqiblten,             &
675                      grid%rthraten,grid%rthratenlw,grid%rthratensw,         &
676                      grid%stepbl,grid%stepra,grid%stepcu,                   &
677                      grid%w0avg, grid%rainnc, grid%rainc, grid%raincv, grid%rainncv,  &
678                      grid%snownc, grid%snowncv, grid%graupelnc, grid%graupelncv,  &
679                      grid%nca,grid%swrad_scat,                    &
680                      grid%cldefi,grid%lowlyr,                          &
681                      grid%mass_flux,                              &
682                      grid%rthften, grid%rqvften,                       &
683                      grid%cldfra,                                      &
684#ifdef WRF_CHEM
685                      grid%cldfra_old,                                  &
686#endif
687#ifndef WRF_CHEM
688                      cldfra_old,                                  &
689#endif
690                      grid%glw,grid%gsw,grid%emiss,grid%embck,            &
691                      grid%lu_index,                                      &
692                      grid%landuse_ISICE, grid%landuse_LUCATS,            &
693                      grid%landuse_LUSEAS, grid%landuse_ISN,              &
694                      grid%lu_state,                                      &
695                      grid%xlat,grid%xlong,grid%albedo,grid%albbck,grid%GMT,grid%JULYR,grid%JULDAY,     &
696                      grid%levsiz, num_ozmixm, num_aerosolc, grid%paerlev,  &
697                      grid%tmn,grid%xland,grid%znt,grid%z0,grid%ust,grid%mol,grid%pblh,grid%tke_pbl,    &
698                      grid%exch_h,grid%thc,grid%snowc,grid%mavail,grid%hfx,grid%qfx,grid%rainbl, &
699                      grid%tslb,grid%zs,grid%dzs,config_flags%num_soil_layers,grid%warm_rain,  &
700                      grid%adv_moist_cond,                         &
701                      grid%apr_gr,grid%apr_w,grid%apr_mc,grid%apr_st,grid%apr_as,      &
702                      grid%apr_capma,grid%apr_capme,grid%apr_capmi,          &
703                      grid%xice,grid%xicem,grid%vegfra,grid%snow,grid%canwat,grid%smstav,         &
704                      grid%smstot, grid%sfcrunoff,grid%udrunoff,grid%grdflx,grid%acsnow,      &
705                      grid%acsnom,grid%ivgtyp,grid%isltyp, grid%sfcevp,grid%smois,     &
706                      grid%sh2o, grid%snowh, grid%smfr3d,                    &
707                      grid%snoalb,                 &
708                      grid%DX,grid%DY,grid%f_ice_phy,grid%f_rain_phy,grid%f_rimef_phy, &
709                      grid%mp_restart_state,grid%tbpvs_state,grid%tbpvs0_state,&
710                      allowed_to_read, grid%moved, start_of_simulation,               &
711                      grid%LAGDAY, &
712                      ids, ide, jds, jde, kds, kde,           &
713                      ims, ime, jms, jme, kms, kme,           &
714                      grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), kts, kte, &
715                      config_flags%num_urban_layers,                        & !multi-layer urban
716                      ozmixm,grid%pin,                             &     ! Optional
717                      grid%m_ps_1,grid%m_ps_2,grid%m_hybi,aerosolc_1,aerosolc_2,&  ! Optional
718                      grid%rundgdten,grid%rvndgdten,grid%rthndgdten,         &     ! Optional
719                      grid%rphndgdten,grid%rqvndgdten,grid%rmundgdten,       &     ! Optional
720                      grid%FGDT,grid%stepfg,                        &     ! Optional
721                      grid%cugd_tten,grid%cugd_ttens,grid%cugd_qvten,        &   ! Optional
722                      grid%cugd_qvtens,grid%cugd_qcten,                 &   ! Optional
723                      grid%DZR, grid%DZB, grid%DZG,                          & !Optional urban
724                      grid%TR_URB2D,grid%TB_URB2D,grid%TG_URB2D,grid%TC_URB2D,    & !Optional urban
725                      grid%QC_URB2D, grid%XXXR_URB2D,grid%XXXB_URB2D,        & !Optional urban
726                      grid%XXXG_URB2D, grid%XXXC_URB2D,                 & !Optional urban
727                      grid%TRL_URB3D, grid%TBL_URB3D, grid%TGL_URB3D,        & !Optional urban
728                      grid%SH_URB2D, grid%LH_URB2D, grid%G_URB2D, grid%RN_URB2D,  & !Optional urban
729                      grid%TS_URB2D, grid%FRC_URB2D, grid%UTYPE_URB2D,      & !Optional urban
730                      grid%TRB_URB4D,grid%TW1_URB4D,grid%TW2_URB4D,grid%TGB_URB4D,grid%TLEV_URB3D,  & !multi-layer urban
731                      grid%QLEV_URB3D,grid%TW1LEV_URB3D,grid%TW2LEV_URB3D,        & !multi-layer urban
732                      grid%TGLEV_URB3D,grid%TFLEV_URB3D,grid%SF_AC_URB3D,         & !multi-layer urban
733                      grid%LF_AC_URB3D,grid%CM_AC_URB3D,grid%SFVENT_URB3D,grid%LFVENT_URB3D,   & !multi-layer urban
734                      grid%SFWIN1_URB3D,grid%SFWIN2_URB3D, & !multi-layer urban
735                      grid%SFW1_URB3D,grid%SFW2_URB3D,grid%SFR_URB3D,grid%SFG_URB3D,     & !multi-layer urban
736                      grid%A_U_BEP,grid%A_V_BEP,grid%A_T_BEP,grid%A_Q_BEP,             & !multi-layer urban
737                      grid%A_E_BEP,grid%B_U_BEP,grid%B_V_BEP,grid%B_T_BEP,             & !multi-layer urban
738                      grid%B_Q_BEP,grid%B_E_BEP,grid%DLG_BEP,                          & !multi-layer urban
739                      grid%DL_U_BEP,grid%SF_BEP,grid%VL_BEP,                           & !multi-layer urban
740                      grid%TML,grid%T0ML,grid%HML,grid%H0ML,grid%HUML,grid%HVML,grid%TMOML,     & !Optional oml
741                      grid%itimestep, grid%fdob,            &
742                      t00, p00, a,                      & ! for obs_nudge base state
743                      grid%TYR, grid%TYRA, grid%TDLY, grid%TLAG, grid%NYEAR, grid%NDAY,grid%tmn_update, &
744                      grid%achfx, grid%aclhf, grid%acgrdflx                 &
745                      ,grid%te_temf,grid%cf3d_temf,grid%wm_temf        & ! WA
746                      )
747       ENDIF   !tgs
748
749   ENDDO
750
751   CALL wrf_debug ( 100 , 'start_domain_em: After call to phy_init' )
752
753   IF (config_flags%do_avgflx_em .EQ. 1) THEN
754      WRITE ( message , FMT = '("start_em: initializing avgflx on domain ",I3)' ) &
755           & grid%id
756      CALL wrf_message(trim(message))
757      grid%avgflx_count = 0
758      f_flux = config_flags%do_avgflx_cugd .EQ. 1  ! WA 9/24/10
759      DO ij = 1, grid%num_tiles
760         call wrf_debug(200,'In start_em, before zero_avgflx call')
761         if (.not. grid%restart) call zero_avgflx(grid%avgflx_rum,grid%avgflx_rvm,grid%avgflx_wwm, &
762              &   ids, ide, jds, jde, kds, kde,           &
763              &   ims, ime, jms, jme, kms, kme,           &
764              &   grid%i_start(ij), grid%i_end(ij), grid%j_start(ij), grid%j_end(ij), kts, kte, f_flux, &
765              &   grid%avgflx_cfu1,grid%avgflx_cfd1,grid%avgflx_dfu1, &
766              &   grid%avgflx_efu1,grid%avgflx_dfd1,grid%avgflx_efd1 )
767         call wrf_debug(200,'In start_em, after zero_avgflx call')
768      ENDDO
769   ENDIF
770
771#ifdef MCELIO
772   grid%LU_MASK = 0.
773   WHERE ( grid%lu_index .EQ. 16 ) grid%LU_MASK = 1.
774#endif
775
776   END IF
777
778#if 0
779#include "CYCLE_TEST.inc"
780#endif
781
782!
783!
784
785 !  set physical boundary conditions for all initialized variables
786
787!-----------------------------------------------------------------------
788!  Stencils for patch communications  (WCS, 29 June 2001)
789!  Note:  the size of this halo exchange reflects the
790!         fact that we are carrying the uncoupled variables
791!         as state variables in the mass coordinate model, as
792!         opposed to the coupled variables as in the height
793!         coordinate model.
794!
795!                           * * * * *
796!         *        * * *    * * * * *
797!       * + *      * + *    * * + * *
798!         *        * * *    * * * * *
799!                           * * * * *
800!
801!j  grid%u_1                          x
802!j  grid%u_2                          x
803!j  grid%v_1                          x
804!j  grid%v_2                          x
805!j  grid%w_1                          x
806!j  grid%w_2                          x
807!j  grid%t_1                          x
808!j  grid%t_2                          x
809!j grid%ph_1                          x
810!j grid%ph_2                          x
811!
812!j  grid%t_init                       x
813!
814!j  grid%phb   x
815!j  grid%ph0   x
816!j  grid%php   x
817!j   grid%pb   x
818!j   grid%al   x
819!j  grid%alt   x
820!j  grid%alb   x
821!
822!  the following are 2D (xy) variables
823!
824!j grid%mu_1                          x
825!j grid%mu_2                          x
826!j grid%mub   x
827!j grid%mu0   x
828!j grid%ht    x
829!j grid%msftx x
830!j grid%msfty x
831!j grid%msfux x
832!j grid%msfuy x
833!j grid%msfvx x
834!j grid%msfvy x
835!j grid%sina  x
836!j grid%cosa  x
837!j grid%e     x
838!j grid%f     x
839!
840!  4D variables
841!
842! moist                        x
843!  chem                        x
844!scalar                        x
845
846!--------------------------------------------------------------
847
848#ifdef DM_PARALLEL
849#  include "HALO_EM_INIT_1.inc"
850#  include "HALO_EM_INIT_2.inc"
851#  include "HALO_EM_INIT_3.inc"
852#  include "HALO_EM_INIT_4.inc"
853#  include "HALO_EM_INIT_5.inc"
854#  include "PERIOD_BDY_EM_INIT.inc"
855#  include "PERIOD_BDY_EM_MOIST.inc"
856#  include "PERIOD_BDY_EM_TKE.inc"
857#  include "PERIOD_BDY_EM_SCALAR.inc"
858#  include "PERIOD_BDY_EM_CHEM.inc"
859#endif
860
861
862   CALL set_physical_bc3d( grid%u_1 , 'U' , config_flags ,                  &
863                         ids , ide , jds , jde , kds , kde ,        &
864                         ims , ime , jms , jme , kms , kme ,        &
865                         its , ite , jts , jte , kts , kte ,        &
866                         its , ite , jts , jte , kts , kte )
867   CALL set_physical_bc3d( grid%u_2 , 'U' , config_flags ,                  &
868                         ids , ide , jds , jde , kds , kde ,        &
869                         ims , ime , jms , jme , kms , kme ,        &
870                         its , ite , jts , jte , kts , kte ,        &
871                         its , ite , jts , jte , kts , kte )
872
873   CALL set_physical_bc3d( grid%v_1 , 'V' , config_flags ,                  &
874                         ids , ide , jds , jde , kds , kde ,        &
875                         ims , ime , jms , jme , kms , kme ,        &
876                         its , ite , jts , jte , kts , kte ,        &
877                         its , ite , jts , jte , kts , kte )
878   CALL set_physical_bc3d( grid%v_2 , 'V' , config_flags ,                  &
879                         ids , ide , jds , jde , kds , kde ,        &
880                         ims , ime , jms , jme , kms , kme ,        &
881                         its , ite , jts , jte , kts , kte ,        &
882                         its , ite , jts , jte , kts , kte )
883
884! set kinematic condition for w
885
886   CALL set_physical_bc2d( grid%ht , 'r' , config_flags ,                &
887                           ids , ide , jds , jde , &
888                           ims , ime , jms , jme , &
889                           its , ite , jts , jte , &
890                           its , ite , jts , jte   )
891
892! Set the w at the surface.  If this is the start of a forecast, or if this is a cycled run
893! and the start of that forecast, we define the w field.  However, if this a restart, then
894! we leave w alone.  Initial value is V dot grad(topo) at surface, then smoothly decaying
895! above that.
896
897   IF ( ( start_of_simulation .OR. config_flags%cycling ) .AND. ( .NOT. config_flags%restart ) ) THEN
898   fill_w_flag = .true.
899   CALL set_w_surface(  config_flags, grid%znw, fill_w_flag,             &
900                        grid%w_1, grid%ht, grid%u_1, grid%v_1, grid%cf1, &
901                        grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
902                        ids, ide, jds, jde, kds, kde,                    &
903                        ims, ime, jms, jme, kms, kme,                    &
904                        its, ite, jts, jte, kts, kte                     )
905   CALL set_w_surface(  config_flags, grid%znw, fill_w_flag,             &
906                        grid%w_2, grid%ht, grid%u_2, grid%v_2, grid%cf1, &
907                        grid%cf2, grid%cf3, grid%rdx, grid%rdy, grid%msftx, grid%msfty, &
908                        ids, ide, jds, jde, kds, kde,                    &
909                        ims, ime, jms, jme, kms, kme,                    &
910                        its, ite, jts, jte, kts, kte                     )
911
912! finished setting kinematic condition for w at the surface
913
914! set up slope-radiation constant arrays based on topography
915       DO j = jts,min(jte,jde-1)
916       DO i = its, min(ite,ide-1)
917          im1 = max(i-1,ids)
918          ip1 = min(i+1,ide-1)
919          jm1 = max(j-1,jds)
920          jp1 = min(j+1,jde-1)
921          grid%toposlpx(i,j)=(grid%ht(ip1,j)-grid%ht(im1,j))*grid%msftx(i,j)*grid%rdx/(ip1-im1)
922          grid%toposlpy(i,j)=(grid%ht(i,jp1)-grid%ht(i,jm1))*grid%msfty(i,j)*grid%rdy/(jp1-jm1)
923             hx = grid%toposlpx(i,j)
924             hy = grid%toposlpy(i,j)
925             pi = 4.*atan(1.)
926             grid%slope(i,j) = atan((hx**2+hy**2)**.5)
927             if (grid%slope(i,j).lt.1.e-4) then
928               grid%slope(i,j) = 0.
929               grid%slp_azi(i,j) = 0.
930             else
931               grid%slp_azi(i,j) = atan2(hx,hy)+pi
932
933! Rotate slope azimuth to lat-lon grid
934               if (grid%cosa(i,j).ge.0) then
935                 grid%slp_azi(i,j) = grid%slp_azi(i,j) - asin(grid%sina(i,j))
936               else
937                 grid%slp_azi(i,j) = grid%slp_azi(i,j) - (pi - asin(grid%sina(i,j)))
938               endif
939             endif
940       ENDDO
941       ENDDO
942   END IF
943
944   CALL set_physical_bc3d( grid%w_1 , 'W' , config_flags ,                 &
945                         ids , ide , jds , jde , kds , kde ,        &
946                         ims , ime , jms , jme , kms , kme ,        &
947                         its , ite , jts , jte , kts , kte ,        &
948                         its , ite , jts , jte , kts , kte )
949   CALL set_physical_bc3d( grid%w_2 , 'W' , config_flags ,                 &
950                         ids , ide , jds , jde , kds , kde ,        &
951                         ims , ime , jms , jme , kms , kme ,        &
952                         its , ite , jts , jte , kts , kte ,        &
953                         its , ite , jts , jte , kts , kte )
954
955   CALL set_physical_bc3d( grid%ph_1 , 'W' , config_flags ,                 &
956                         ids , ide , jds , jde , kds , kde ,        &
957                         ims , ime , jms , jme , kms , kme ,        &
958                         its , ite , jts , jte , kts , kte ,        &
959                         its , ite , jts , jte , kts , kte )
960
961   CALL set_physical_bc3d( grid%ph_2 , 'W' , config_flags ,                 &
962                         ids , ide , jds , jde , kds , kde ,        &
963                         ims , ime , jms , jme , kms , kme ,        &
964                         its , ite , jts , jte , kts , kte ,        &
965                         its , ite , jts , jte , kts , kte )
966
967   CALL set_physical_bc3d( grid%t_1 , 't' , config_flags ,                 &
968                         ids , ide , jds , jde , kds , kde ,        &
969                         ims , ime , jms , jme , kms , kme ,        &
970                         its , ite , jts , jte , kts , kte ,        &
971                         its , ite , jts , jte , kts , kte )
972
973   CALL set_physical_bc3d( grid%t_2 , 't' , config_flags ,                 &
974                         ids , ide , jds , jde , kds , kde ,        &
975                         ims , ime , jms , jme , kms , kme ,        &
976                         its , ite , jts , jte , kts , kte ,        &
977                         its , ite , jts , jte , kts , kte )
978
979   CALL set_physical_bc2d( grid%mu_1, 't' , config_flags ,   &
980                           ids , ide , jds , jde ,  &
981                           ims , ime , jms , jme ,  &
982                           its , ite , jts , jte ,  &
983                           its , ite , jts , jte   )
984   CALL set_physical_bc2d( grid%mu_2, 't' , config_flags ,   &
985                           ids , ide , jds , jde ,  &
986                           ims , ime , jms , jme ,  &
987                           its , ite , jts , jte ,  &
988                           its , ite , jts , jte   )
989   CALL set_physical_bc2d( grid%mub , 't' , config_flags ,   &
990                           ids , ide , jds , jde ,  &
991                           ims , ime , jms , jme ,  &
992                           its , ite , jts , jte ,  &
993                           its , ite , jts , jte   )
994!  CALL set_physical_bc2d( grid%mu0 , 't' , config_flags ,   &
995!                          ids , ide , jds , jde ,  &
996!                          ims , ime , jms , jme ,  &
997!                          its , ite , jts , jte ,  &
998!                          its , ite , jts , jte   )
999
1000
1001   CALL set_physical_bc3d( grid%phb , 'W' , config_flags ,                 &
1002                         ids , ide , jds , jde , kds , kde ,        &
1003                         ims , ime , jms , jme , kms , kme ,        &
1004                         its , ite , jts , jte , kts , kte ,        &
1005                         its , ite , jts , jte , kts , kte )
1006   CALL set_physical_bc3d( grid%ph0 , 'W' , config_flags ,                 &
1007                         ids , ide , jds , jde , kds , kde ,        &
1008                         ims , ime , jms , jme , kms , kme ,        &
1009                         its , ite , jts , jte , kts , kte ,        &
1010                         its , ite , jts , jte , kts , kte )
1011   CALL set_physical_bc3d( grid%php , 'W' , config_flags ,                 &
1012                         ids , ide , jds , jde , kds , kde ,        &
1013                         ims , ime , jms , jme , kms , kme ,        &
1014                         its , ite , jts , jte , kts , kte ,        &
1015                         its , ite , jts , jte , kts , kte )
1016
1017   CALL set_physical_bc3d( grid%pb , 't' , config_flags ,                 &
1018                         ids , ide , jds , jde , kds , kde ,        &
1019                         ims , ime , jms , jme , kms , kme ,        &
1020                         its , ite , jts , jte , kts , kte ,        &
1021                         its , ite , jts , jte , kts , kte )
1022   CALL set_physical_bc3d( grid%al , 't' , config_flags ,                 &
1023                         ids , ide , jds , jde , kds , kde ,        &
1024                         ims , ime , jms , jme , kms , kme ,        &
1025                         its , ite , jts , jte , kts , kte ,        &
1026                         its , ite , jts , jte , kts , kte )
1027   CALL set_physical_bc3d( grid%alt , 't' , config_flags ,                 &
1028                         ids , ide , jds , jde , kds , kde ,        &
1029                         ims , ime , jms , jme , kms , kme ,        &
1030                         its , ite , jts , jte , kts , kte ,        &
1031                         its , ite , jts , jte , kts , kte )
1032   CALL set_physical_bc3d( grid%alb , 't' , config_flags ,                 &
1033                         ids , ide , jds , jde , kds , kde ,        &
1034                         ims , ime , jms , jme , kms , kme ,        &
1035                         its , ite , jts , jte , kts , kte ,        &
1036                         its , ite , jts , jte , kts , kte )
1037   CALL set_physical_bc3d(grid%t_init, 't' , config_flags ,                 &
1038                         ids , ide , jds , jde , kds , kde ,        &
1039                         ims , ime , jms , jme , kms , kme ,        &
1040                         its , ite , jts , jte , kts , kte ,        &
1041                         its , ite , jts , jte , kts , kte )
1042   CALL set_physical_bc3d(grid%tke_2, 't' , config_flags ,                 &
1043                         ids , ide , jds , jde , kds , kde ,        &
1044                         ims , ime , jms , jme , kms , kme ,        &
1045                         its , ite , jts , jte , kts , kte ,        &
1046                         its , ite , jts , jte , kts , kte )
1047
1048   IF (num_moist > 0) THEN
1049
1050! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
1051
1052      loop_3d_m   : DO loop = 1 , num_moist
1053         CALL set_physical_bc3d( moist(:,:,:,loop) , 'r' , config_flags ,                 &
1054                                 ids , ide , jds , jde , kds , kde ,        &
1055                                 ims , ime , jms , jme , kms , kme ,        &
1056                                 its , ite , jts , jte , kts , kte ,        &
1057                                 its , ite , jts , jte , kts , kte )
1058      END DO loop_3d_m
1059
1060   ENDIF
1061
1062   IF (num_scalar > 0) THEN
1063
1064! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
1065
1066      loop_3d_s   : DO loop = 1 , num_scalar
1067         CALL set_physical_bc3d( scalar(:,:,:,loop) , 'r' , config_flags ,                 &
1068                                 ids , ide , jds , jde , kds , kde ,        &
1069                                 ims , ime , jms , jme , kms , kme ,        &
1070                                 its , ite , jts , jte , kts , kte ,        &
1071                                 its , ite , jts , jte , kts , kte )
1072      END DO loop_3d_s
1073
1074   ENDIF
1075
1076
1077#ifdef WRF_CHEM
1078        if(config_flags%tracer_opt > 0 )then
1079           call initialize_tracer (tracer,config_flags%chem_in_opt, &
1080                               config_flags%tracer_opt,num_tracer,  &
1081                               ids,ide, jds,jde, kds,kde,      & ! domain dims
1082                               ims,ime, jms,jme, kms,kme,      & ! memory dims
1083                               ips,ipe, jps,jpe, kps,kpe,      & ! patch  dims
1084                               its,ite, jts,jte, kts,kte )
1085
1086        endif
1087!
1088!   we do this here, so we only have one chem_init routine for either core....
1089!
1090         do j=jts,min(jte,jde-1)
1091            do i=its,min(ite,ide-1)
1092               do k=kts,kte
1093                 z_at_w(i,k,j)=(grid%ph_2(i,k,j)+grid%phb(i,k,j))/g
1094               enddo
1095               do k=kts,min(kte,kde-1)
1096                 tempfac=(grid%t_1(i,k,j) + t0)*((grid%p(i,k,j) + grid%pb(i,k,j))/p1000mb)**rcp
1097                 convfac(i,k,j) = (grid%p(i,k,j)+grid%pb(i,k,j))/rgasuniv/tempfac
1098               enddo
1099            enddo
1100         enddo
1101
1102         CALL chem_init (grid%id,chem,emis_ant,scalar,grid%dt,grid%bioemdt,grid%photdt,     &
1103                         grid%chemdt,                                       &
1104                         grid%stepbioe,grid%stepphot,grid%stepchem,grid%stepfirepl,      &
1105                         grid%plumerisefire_frq,z_at_w,grid%xlat,grid%xlong,g,           &
1106                         grid%aerwrf,config_flags,grid,            &
1107                         grid%alt,grid%t_1,grid%p,convfac,grid%ttday,grid%tcosz,         &
1108                         grid%julday,grid%gmt,&
1109                         grid%gd_cloud, grid%gd_cloud2,grid%raincv_a,grid%raincv_b,     &
1110                         grid%gd_cloud_a, grid%gd_cloud2_a,                              &
1111                         grid%gd_cloud_b, grid%gd_cloud2_b,                              &
1112                         grid%tauaer1,grid%tauaer2,grid%tauaer3,grid%tauaer4,                   &
1113                         grid%gaer1,grid%gaer2,grid%gaer3,grid%gaer4,                           &
1114                         grid%waer1,grid%waer2,grid%waer3,grid%waer4,                           &
1115                         grid%l2aer,grid%l3aer,grid%l4aer,grid%l5aer,grid%l6aer,grid%l7aer,     &
1116                         grid%extaerlw1,grid%extaerlw2,grid%extaerlw3,grid%extaerlw4,           &
1117                         grid%extaerlw5,grid%extaerlw6,grid%extaerlw7,grid%extaerlw8,           &
1118                         grid%extaerlw9,grid%extaerlw10,grid%extaerlw11,grid%extaerlw12,        &
1119                         grid%extaerlw13,grid%extaerlw14,grid%extaerlw15,grid%extaerlw16,       &
1120                         grid%tauaerlw1,grid%tauaerlw2,grid%tauaerlw3,grid%tauaerlw4,           &
1121                         grid%tauaerlw5,grid%tauaerlw6,grid%tauaerlw7,grid%tauaerlw8,           &
1122                         grid%tauaerlw9,grid%tauaerlw10,grid%tauaerlw11,grid%tauaerlw12,        &
1123                         grid%tauaerlw13,grid%tauaerlw14,grid%tauaerlw15,grid%tauaerlw16,       &
1124                         grid%pm2_5_dry,grid%pm2_5_water,grid%pm2_5_dry_ec,                &
1125                         grid%last_chem_time_year,grid%last_chem_time_month,               &
1126                         grid%last_chem_time_day,grid%last_chem_time_hour,                 &
1127                         grid%last_chem_time_minute,grid%last_chem_time_second,            &
1128                         grid%chem_in_opt,grid%kemit,                       &
1129                         ids , ide , jds , jde , kds , kde ,                &
1130                         ims , ime , jms , jme , kms , kme ,                &
1131                         its , ite , jts , jte , kts , kte                  )
1132
1133!
1134! calculate initial pm
1135!
1136!       print *,'calculating initial pm'
1137        select case (config_flags%chem_opt)
1138        case (GOCART_SIMPLE,GOCARTRACM_KPP,GOCARTRADM2,GOCARTRADM2_KPP)
1139           call sum_pm_gocart (                                             &
1140                grid%alt, chem, grid%pm2_5_dry, grid%pm2_5_dry_ec,grid%pm10,&
1141                ids,ide, jds,jde, kds,kde,                                  &
1142                ims,ime, jms,jme, kms,kme,                                  &
1143                its,ite, jts,jte, kts,kte-1                                 )
1144        case (RADM2SORG,RADM2SORG_AQ,RADM2SORG_KPP,RACMSORG_AQ,RACMSORG_KPP)
1145           call sum_pm_sorgam (                                             &
1146                grid%alt, chem, grid%h2oaj, grid%h2oai,                                    &
1147                grid%pm2_5_dry, grid%pm2_5_water, grid%pm2_5_dry_ec, grid%pm10,                 &
1148                grid%dust_opt,ids,ide, jds,jde, kds,kde,                    &
1149                ims,ime, jms,jme, kms,kme,                                  &
1150                its,ite, jts,jte, kts,kte-1                                 )
1151
1152        case (CBMZ_MOSAIC_4BIN,CBMZ_MOSAIC_8BIN,CBMZ_MOSAIC_4BIN_AQ,CBMZ_MOSAIC_8BIN_AQ, &
1153              CBMZ_MOSAIC_DMS_4BIN, CBMZ_MOSAIC_DMS_8BIN, CBMZ_MOSAIC_DMS_4BIN_AQ,       &
1154              CBMZ_MOSAIC_DMS_8BIN_AQ)
1155           call sum_pm_mosaic (                                             &
1156                grid%alt, chem,                                                  &
1157                grid%pm2_5_dry, grid%pm2_5_water, grid%pm2_5_dry_ec, grid%pm10,                 &
1158                ids,ide, jds,jde, kds,kde,                                  &
1159                ims,ime, jms,jme, kms,kme,                                  &
1160                its,ite, jts,jte, kts,kte-1                                 )
1161
1162        case default
1163           do j=jts,min(jte,jde-1)
1164              do k=kts,min(kte,kde-1)
1165                 do i=its,min(ite,ide-1)
1166                    grid%pm2_5_dry(i,k,j)    = 0.
1167                    grid%pm2_5_water(i,k,j)  = 0.
1168                    grid%pm2_5_dry_ec(i,k,j) = 0.
1169                    grid%pm10(i,k,j)         = 0.
1170                 enddo
1171              enddo
1172           enddo
1173        end select
1174       
1175        ! initialize advective tendency diagnostics for chem
1176        if ( grid%itimestep .eq. 0 .and. config_flags%chemdiag .eq. USECHEMDIAG ) then
1177          grid%advh_ct(:,:,:,:) = 0.
1178          grid%advz_ct(:,:,:,:) = 0.
1179        endif
1180#endif
1181
1182        ! initialize advective tendency diagnostics for non-chem
1183        if ( grid%itimestep .eq. 0 .and. config_flags%tenddiag .eq. USETENDDIAG ) then
1184          grid%advh_t(:,:,:,:) = 0.
1185          grid%advz_t(:,:,:,:) = 0.
1186        endif
1187
1188   IF (num_chem >= PARAM_FIRST_SCALAR ) THEN
1189! use of (:,:,:,loop) not efficient on DEC, but (ims,kms,jms,loop) not portable to SGI/Cray
1190
1191      loop_3d_c   : DO loop = PARAM_FIRST_SCALAR , num_chem
1192         CALL set_physical_bc3d( chem(:,:,:,loop) , 'r' , config_flags ,                 &
1193                                 ids , ide , jds , jde , kds , kde ,        &
1194                                 ims , ime , jms , jme , kms , kme ,        &
1195                                 its , ite , jts , jte , kts , kte ,        &
1196                                 its , ite , jts , jte , kts , kte )
1197      END DO loop_3d_c
1198
1199   ENDIF
1200
1201   CALL set_physical_bc2d( grid%msftx , 'r' , config_flags ,              &
1202                         ids , ide , jds , jde , &
1203                         ims , ime , jms , jme , &
1204                         its , ite , jts , jte , &
1205                         its , ite , jts , jte   )
1206   CALL set_physical_bc2d( grid%msfty , 'r' , config_flags ,              &
1207                         ids , ide , jds , jde , &
1208                         ims , ime , jms , jme , &
1209                         its , ite , jts , jte , &
1210                         its , ite , jts , jte   )
1211   CALL set_physical_bc2d( grid%msfux , 'x' , config_flags ,              &
1212                         ids , ide , jds , jde , &
1213                         ims , ime , jms , jme , &
1214                         its , ite , jts , jte , &
1215                         its , ite , jts , jte   )
1216   CALL set_physical_bc2d( grid%msfuy , 'x' , config_flags ,              &
1217                         ids , ide , jds , jde , &
1218                         ims , ime , jms , jme , &
1219                         its , ite , jts , jte , &
1220                         its , ite , jts , jte   )
1221   CALL set_physical_bc2d( grid%msfvx , 'y' , config_flags ,              &
1222                         ids , ide , jds , jde , &
1223                         ims , ime , jms , jme , &
1224                         its , ite , jts , jte , &
1225                         its , ite , jts , jte   )
1226   CALL set_physical_bc2d( grid%msfvy , 'y' , config_flags ,              &
1227                         ids , ide , jds , jde , &
1228                         ims , ime , jms , jme , &
1229                         its , ite , jts , jte , &
1230                         its , ite , jts , jte   )
1231   CALL set_physical_bc2d( grid%sina , 'r' , config_flags ,              &
1232                         ids , ide , jds , jde , &
1233                         ims , ime , jms , jme , &
1234                         its , ite , jts , jte , &
1235                         its , ite , jts , jte   )
1236   CALL set_physical_bc2d( grid%cosa , 'r' , config_flags ,              &
1237                         ids , ide , jds , jde , &
1238                         ims , ime , jms , jme , &
1239                         its , ite , jts , jte , &
1240                         its , ite , jts , jte   )
1241   CALL set_physical_bc2d( grid%e , 'r' , config_flags ,                 &
1242                         ids , ide , jds , jde , &
1243                         ims , ime , jms , jme , &
1244                         its , ite , jts , jte , &
1245                         its , ite , jts , jte   )
1246   CALL set_physical_bc2d( grid%f , 'r' , config_flags ,                 &
1247                         ids , ide , jds , jde , &
1248                         ims , ime , jms , jme , &
1249                         its , ite , jts , jte , &
1250                         its , ite , jts , jte   )
1251
1252#ifndef WRF_CHEM
1253      DEALLOCATE(CLDFRA_OLD)
1254#endif
1255#ifdef DM_PARALLEL
1256#   include "HALO_EM_INIT_1.inc"
1257#   include "HALO_EM_INIT_2.inc"
1258#   include "HALO_EM_INIT_3.inc"
1259#   include "HALO_EM_INIT_4.inc"
1260#   include "HALO_EM_INIT_5.inc"
1261#  include "PERIOD_BDY_EM_INIT.inc"
1262#  include "PERIOD_BDY_EM_MOIST.inc"
1263#  include "PERIOD_BDY_EM_TKE.inc"
1264#  include "PERIOD_BDY_EM_SCALAR.inc"
1265#  include "PERIOD_BDY_EM_CHEM.inc"
1266#endif
1267
1268!  FIRE
1269if(config_flags%ifire.eq.2)then
1270
1271   call sfire_driver_em_init ( grid , config_flags    &
1272            ,ids,ide, kds,kde, jds,jde                &
1273            ,ims,ime, kms,kme, jms,jme                &
1274            ,ips,ipe, kps,kpe, jps,jpe )
1275
1276   CALL wrf_debug ( 100 , 'start_domain_em: After call to sfire_driver_em_init' )
1277endif
1278
1279
1280     CALL wrf_debug ( 100 , 'start_domain_em: Returning' )
1281
1282     RETURN
1283
1284   END SUBROUTINE start_domain_em
1285
Note: See TracBrowser for help on using the repository browser.