Merge pull request #22 from wirc-sjsu/develop-w21
[WRF-Fire-merge.git] / share / wrf_timeseries.F
blobcab5a0aa10c131ec041a7becd6c2819646027ceb
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 ! This routine prints out the current value of variables at all specified
3 !   time series locations that are within the current patch.
5 ! Michael G. Duda -- 25 August 2005
6 ! vertical profiles added by Torge Lorenz -- June 2012
7 ! ability to output at either i/j or lat/lon locations, and ability to 
8 !     output W, added by Pat Hawbecker -- Jan 2019
9 ! ability to output O3 (chem), added by Xin Zhang -- Feb 2020
10
11 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
12 SUBROUTINE calc_ts_locations( grid )
14    USE module_domain, ONLY : domain, get_ijk_from_grid
15    USE module_configure, ONLY : model_config_rec, grid_config_rec_type, model_to_grid_config_rec
16    USE module_dm, ONLY : wrf_dm_min_real
17    USE module_llxy
18    USE module_state_description
19    USE module_model_constants
21    IMPLICIT NONE
23    ! Arguments
24    TYPE (domain), INTENT(INOUT) :: grid
26    ! Externals
27    LOGICAL, EXTERNAL :: wrf_dm_on_monitor
28    INTEGER, EXTERNAL :: get_unused_unit
30    ! Local variables
31    INTEGER :: ntsloc_temp
32    INTEGER :: i, j, k, iunit
33    REAL :: ts_rx, ts_ry, ts_xlat, ts_xlong, ts_hgt
34    REAL :: known_lat, known_lon
35    CHARACTER (LEN=132) :: message
36    CHARACTER (LEN=24) :: ts_profile_filename
37 #if (WRF_CHEM == 1)
38    INTEGER, PARAMETER :: TS_FIELDS = 8
39    CHARACTER (LEN=2), DIMENSION(TS_FIELDS) :: &
40       ts_file_endings = (/ 'UU', 'VV', 'PH', 'TH', 'QV' ,'WW', 'PR', 'O3'/)
41 #else
42    INTEGER, PARAMETER :: TS_FIELDS = 7
43    CHARACTER (LEN=2), DIMENSION(TS_FIELDS) :: &
44       ts_file_endings = (/ 'UU', 'VV', 'PH', 'TH', 'QV' ,'WW', 'PR'/) 
45 #endif
46    INTEGER   ierr
47    CHARACTER (len=19) simulation_start_date
48    INTEGER simulation_start_year   , &
49            simulation_start_month  , &
50            simulation_start_day    , &
51            simulation_start_hour   , &
52            simulation_start_minute , &
53            simulation_start_second
55    TYPE (PROJ_INFO) :: ts_proj
56    TYPE (grid_config_rec_type) :: config_flags
58    INTEGER :: ids, ide, jds, jde, kds, kde,        &
59               ims, ime, jms, jme, kms, kme,        &
60               ips, ipe, jps, jpe, kps, kpe,        &
61               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
67    IF ( grid%ntsloc .LE. 0 ) RETURN
69 #if ((EM_CORE == 1) && (DA_CORE != 1))
70    IF ( grid%dfi_stage == DFI_FST ) THEN
71 #endif
72       CALL get_ijk_from_grid ( grid ,                               &
73                                ids, ide, jds, jde, kds, kde,        &
74                                ims, ime, jms, jme, kms, kme,        &
75                                ips, ipe, jps, jpe, kps, kpe,        &
76                                imsx, imex, jmsx, jmex, kmsx, kmex,  &
77                                ipsx, ipex, jpsx, jpex, kpsx, kpex,  &
78                                imsy, imey, jmsy, jmey, kmsy, kmey,  &
79                                ipsy, ipey, jpsy, jpey, kpsy, kpey )
80    
81       CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
82    
83       ! Set up map transformation structure
84       CALL map_init(ts_proj)
85    
86       IF (ips <= 1 .AND. 1 <= ipe .AND. &
87           jps <= 1 .AND. 1 <= jpe) THEN
88          known_lat = grid%xlat(1,1)
89          known_lon = grid%xlong(1,1)
90       ELSE
91          known_lat = 9999.
92          known_lon = 9999.
93       END IF
94       known_lat = wrf_dm_min_real(known_lat)
95       known_lon = wrf_dm_min_real(known_lon)
96    
97       ! Mercator
98       IF (config_flags%map_proj == PROJ_MERC) THEN
99          CALL map_set(PROJ_MERC, ts_proj,               &
100                       truelat1 = config_flags%truelat1, &
101                       lat1     = known_lat,             &
102                       lon1     = known_lon,             &
103                       knowni   = 1.,                    &
104                       knownj   = 1.,                    &
105                       dx       = config_flags%dx)
106    
107       ! Lambert conformal
108       ELSE IF (config_flags%map_proj == PROJ_LC) THEN
109          CALL map_set(PROJ_LC, ts_proj,                  &
110                       truelat1 = config_flags%truelat1,  &
111                       truelat2 = config_flags%truelat2,  &
112                       stdlon   = config_flags%stand_lon, &
113                       lat1     = known_lat,              &
114                       lon1     = known_lon,              &
115                       knowni   = 1.,                     &
116                       knownj   = 1.,                     &
117                       dx       = config_flags%dx)
118    
119       ! Polar stereographic
120       ELSE IF (config_flags%map_proj == PROJ_PS) THEN
121          CALL map_set(PROJ_PS, ts_proj,                  &
122                       truelat1 = config_flags%truelat1,  &
123                       stdlon   = config_flags%stand_lon, &
124                       lat1     = known_lat,              &
125                       lon1     = known_lon,              &
126                       knowni   = 1.,                     &
127                       knownj   = 1.,                     &
128                       dx       = config_flags%dx)
129    
130 #if (EM_CORE == 1)
131       ! Cassini (global ARW)
132       ELSE IF (config_flags%map_proj == PROJ_CASSINI) THEN
133          CALL map_set(PROJ_CASSINI, ts_proj,                            &
134                       latinc   = grid%dy*360.0/(2.0*EARTH_RADIUS_M*PI), &
135                       loninc   = grid%dx*360.0/(2.0*EARTH_RADIUS_M*PI), & 
136                       lat1     = known_lat,                             &
137                       lon1     = known_lon,                             &
138                       lat0     = config_flags%pole_lat,                 &
139                       lon0     = config_flags%pole_lon,                 &
140                       knowni   = 1.,                                    &
141                       knownj   = 1.,                                    &
142                       stdlon   = config_flags%stand_lon)
143 #endif
145       ! Rotated latitude-longitude
146       ELSE IF (config_flags%map_proj == PROJ_ROTLL) THEN
147          CALL map_set(PROJ_ROTLL, ts_proj,                      &
148 ! I have no idea how this should work for NMM nested domains
149                       ixdim    = grid%e_we-1,                   &
150                       jydim    = grid%e_sn-1,                   &
151                       phi      = real(grid%e_sn-2)*grid%dy/2.0, &
152                       lambda   = real(grid%e_we-2)*grid%dx,     &
153                       lat1     = config_flags%cen_lat,          &
154                       lon1     = config_flags%cen_lon,          &
155                       latinc   = grid%dy,                       &
156                       loninc   = grid%dx,                       &
157                       stagger  = HH)
158    
159       END IF
160    
161       ! Determine simulation start time
162       ierr = 0
163       CALL nl_get_simulation_start_year   ( 1 , simulation_start_year   )
164       CALL nl_get_simulation_start_month  ( 1 , simulation_start_month  )
165       CALL nl_get_simulation_start_day    ( 1 , simulation_start_day    )
166       CALL nl_get_simulation_start_hour   ( 1 , simulation_start_hour   )
167       CALL nl_get_simulation_start_minute ( 1 , simulation_start_minute )
168       CALL nl_get_simulation_start_second ( 1 , simulation_start_second )
169       WRITE ( simulation_start_date , FMT = '(I4.4,"-",I2.2,"-",I2.2,"_",I2.2,":",I2.2,":",I2.2)' ) &
170               simulation_start_year,simulation_start_month,simulation_start_day,simulation_start_hour, &
171               simulation_start_minute,simulation_start_second
172 !     WRITE(message,*)'wrf_timeseries: SIMULATION_START_DATE = ',simulation_start_date(1:19)
173 !     CALL wrf_debug ( 0, TRIM( message ) )
174       
175       ! Determine time series locations for domain
176       IF (.NOT. grid%have_calculated_tslocs) THEN
177          grid%have_calculated_tslocs = .TRUE.
178          WRITE(message, '(A43,I3,A15,A19)') 'Computing time series locations for domain ', grid%id, &
179                                             ' starting from ',simulation_start_date
180          CALL wrf_message(message)
181    
182          ntsloc_temp = 0
184          ! Boolean if/else for ij and lat/lon that sets ts_rx(y) to either ij or lat/lon value from tslist
185          DO k=1,grid%ntsloc
186             ! Ideal case (which has a cartesian coordinate) or specified (i,j) in tslist
187             IF (config_flags%map_proj == 0 .OR. grid%tslist_ij) THEN 
188                ts_rx = grid%itsloc(k)  
189                ts_ry = grid%jtsloc(k)
190             ! Real-data case with input locations provided as (lat,lon)
191             ELSE
192                CALL latlon_to_ij(ts_proj, grid%lattsloc(k), grid%lontsloc(k), ts_rx, ts_ry)
193             END IF
194             
196             ntsloc_temp = ntsloc_temp + 1
197             grid%itsloc(ntsloc_temp) = NINT(ts_rx)
198             grid%jtsloc(ntsloc_temp) = NINT(ts_ry)
199             grid%id_tsloc(ntsloc_temp) = k
200    
201             ! Is point outside of domain (or on the edge of domain)?
202             IF (grid%itsloc(ntsloc_temp) < ids .OR. grid%itsloc(ntsloc_temp) > ide .OR. &
203                 grid%jtsloc(ntsloc_temp) < jds .OR. grid%jtsloc(ntsloc_temp) > jde) THEN
204                ntsloc_temp = ntsloc_temp - 1
205    
206             END IF
207    
208          END DO
209    
210          grid%next_ts_time = 1
211    
212          grid%ntsloc_domain = ntsloc_temp
213    
214          DO k=1,grid%ntsloc_domain
216             ! If location is outside of patch, we need to get lat/lon of TS grid cell from another patch
217             IF (grid%itsloc(k) < ips .OR. grid%itsloc(k) > ipe .OR. &
218                 grid%jtsloc(k) < jps .OR. grid%jtsloc(k) > jpe) THEN
219                ts_xlat  = 1.E30
220                ts_xlong = 1.E30
221                ts_hgt   = 1.E30
222             ELSE
223                ts_xlat  = grid%xlat(grid%itsloc(k),grid%jtsloc(k))
224                ts_xlong = grid%xlong(grid%itsloc(k),grid%jtsloc(k))
225 #if (EM_CORE == 1)
226                ts_hgt   = grid%ht(grid%itsloc(k),grid%jtsloc(k))
227 #endif
228             END IF
229 #if DM_PARALLEL
230             ts_xlat  = wrf_dm_min_real(ts_xlat)
231             ts_xlong = wrf_dm_min_real(ts_xlong)
232             ts_hgt   = wrf_dm_min_real(ts_hgt)
233 #endif
234    
235             IF ( wrf_dm_on_monitor() ) THEN
236                iunit = get_unused_unit()
237                IF ( iunit <= 0 ) THEN
238                   CALL wrf_error_fatal('Error in calc_ts_locations: could not find a free Fortran unit.')
239                END IF
240                WRITE(grid%ts_filename(k),'(A)') TRIM(grid%nametsloc(grid%id_tsloc(k)))//'.d00.TS'
241                i = LEN_TRIM(grid%ts_filename(k))
242                WRITE(grid%ts_filename(k)(i-4:i-3),'(I2.2)') grid%id
243                OPEN(UNIT=iunit, FILE=TRIM(grid%ts_filename(k)), FORM='FORMATTED', STATUS='REPLACE')
244 #if (EM_CORE == 1)
245                IF ( .NOT. grid%tslist_ij ) THEN 
246                   WRITE(UNIT=iunit, &
247                         FMT='(A26,I2,I3,A6,A2,F7.3,A1,F8.3,A3,I4,A1,I4,A3,F7.3,A1,F8.3,A2,F6.1,A7,A2,A19)') &
248                         grid%desctsloc(grid%id_tsloc(k))//' ', grid%id, grid%id_tsloc(k), &
249                         ' '//grid%nametsloc(grid%id_tsloc(k)), &
250                         ' (', grid%lattsloc(grid%id_tsloc(k)), ',', grid%lontsloc(grid%id_tsloc(k)), ') (', &
251                         grid%itsloc(k), ',', grid%jtsloc(k), ') (', &
252                         ts_xlat, ',', ts_xlong, ') ', &
253                         ts_hgt,' meters','  ',simulation_start_date(1:19)
254                ELSE
255                   WRITE(UNIT=iunit, &
256                         FMT='(A26,I2,I3,A6,A2,F7.3,A1,F8.3,A3,I4,A1,I4,A3,F7.3,A1,F8.3,A2,F6.1,A7,A2,A19)') &
257                         grid%desctsloc(grid%id_tsloc(k))//' ', grid%id, grid%id_tsloc(k), &
258                         ' '//grid%nametsloc(grid%id_tsloc(k)), &
259                         ' (', ts_xlat, ',', ts_xlong, ') (', &
260                         grid%itsloc(k), ',', grid%jtsloc(k), ') (', &
261                         ts_xlat, ',', ts_xlong, ') ', &
262                         ts_hgt,' meters','  ',simulation_start_date(1:19)
263                END IF
264 #else
265                WRITE(UNIT=iunit, &
266                      FMT='(A26,I2,I3,A6,A2,F7.3,A1,F8.3,A3,I4,A1,I4,A3,F7.3,A1,F8.3,A2)') &
267                      grid%desctsloc(grid%id_tsloc(k))//' ', grid%id, grid%id_tsloc(k), &
268                      ' '//grid%nametsloc(grid%id_tsloc(k)), &
269                      ' (', grid%lattsloc(grid%id_tsloc(k)), ',', grid%lontsloc(grid%id_tsloc(k)), ') (', &
270                      grid%itsloc(k), ',', grid%jtsloc(k), ') (', &
271                      ts_xlat, ',', ts_xlong, ') '
272 #endif
273                CLOSE(UNIT=iunit)
275                ts_profile_filename = grid%ts_filename(k)
276                DO j=1,SIZE(ts_file_endings)
277                   ! Create the output files for the vertical profiles, one file for each variable
278                   iunit = get_unused_unit()
279                   IF ( iunit <= 0 ) THEN
280                      CALL wrf_error_fatal('Error in calc_ts_locations: could not find a free Fortran unit.')
281                   END IF
282                   i = LEN_TRIM(ts_profile_filename)
283                   WRITE(ts_profile_filename(i-1:i),'(A2)') ts_file_endings(j)
284                   OPEN(UNIT=iunit, FILE=TRIM(ts_profile_filename), FORM='FORMATTED', STATUS='REPLACE')
285 #if (EM_CORE == 1)
286                   IF ( .NOT. grid%tslist_ij ) THEN 
287                      WRITE(UNIT=iunit, &
288                            FMT='(A26,I2,I3,A6,A2,F7.3,A1,F8.3,A3,I4,A1,I4,A3,F7.3,A1,F8.3,A2,F6.1,A7,A2,A19)') &
289                            grid%desctsloc(grid%id_tsloc(k))//' ', grid%id, grid%id_tsloc(k), &
290                            ' '//grid%nametsloc(grid%id_tsloc(k)), &
291                            ' (', grid%lattsloc(grid%id_tsloc(k)), ',', grid%lontsloc(grid%id_tsloc(k)), ') (', &
292                            grid%itsloc(k), ',', grid%jtsloc(k), ') (', &
293                            ts_xlat, ',', ts_xlong, ') ', &
294                            ts_hgt,' meters','  ',simulation_start_date
295                   ELSE
296                      WRITE(UNIT=iunit, &
297                            FMT='(A26,I2,I3,A6,A2,F7.3,A1,F8.3,A3,I4,A1,I4,A3,F7.3,A1,F8.3,A2,F6.1,A7,A2,A19)') &
298                            grid%desctsloc(grid%id_tsloc(k))//' ', grid%id, grid%id_tsloc(k), &
299                            ' '//grid%nametsloc(grid%id_tsloc(k)), &
300                            ' (', ts_xlat, ',', ts_xlong, ') (', &
301                            grid%itsloc(k), ',', grid%jtsloc(k), ') (', &
302                            ts_xlat, ',', ts_xlong, ') ', &
303                            ts_hgt,' meters,','  ',simulation_start_date
304                  END IF
305 #else
306                   WRITE(UNIT=iunit, &
307                         FMT='(A26,I2,I3,A6,A2,F7.3,A1,F8.3,A3,I4,A1,I4,A3,F7.3,A1,F8.3,A2)') &
308                         grid%desctsloc(grid%id_tsloc(k))//' ', grid%id, grid%id_tsloc(k), &
309                         ' '//grid%nametsloc(grid%id_tsloc(k)), &
310                         ' (', grid%lattsloc(grid%id_tsloc(k)), ',', grid%lontsloc(grid%id_tsloc(k)), ') (', &
311                         grid%itsloc(k), ',', grid%jtsloc(k), ') (', &
312                         ts_xlat, ',', ts_xlong, ') '
313 #endif               
314                   CLOSE(UNIT=iunit)
315                END DO
316             END IF
317          END DO
318    
319       END IF
320 #if ((EM_CORE == 1) && (DA_CORE != 1))
321    END IF
322 #endif
324 END SUBROUTINE calc_ts_locations
327 SUBROUTINE calc_ts( grid )
329    USE module_domain
330    USE module_configure, ONLY : model_config_rec, grid_config_rec_type, model_to_grid_config_rec
331    USE module_model_constants
333    IMPLICIT NONE
335    ! Arguments
336    TYPE (domain), INTENT(INOUT) :: grid
338    LOGICAL, EXTERNAL :: wrf_dm_on_monitor
340    ! Local variables
341    INTEGER :: i, k, mm, n, ix, iy, rc
342    REAL :: earth_u, earth_v,                       &
343            output_t, output_q, clw, xtime_minutes
344    REAL, PARAMETER :: MISSING = -999.0
345    REAL, ALLOCATABLE, DIMENSION(:) :: p8w
346    REAL, ALLOCATABLE, DIMENSION(:) :: earth_u_profile, earth_v_profile
347    TYPE (grid_config_rec_type) :: config_flags
349    ! Parameter ts_model_level:  
350        ! TRUE to output T, Q, and wind at lowest model level
351        ! FALSE to output T and Q at 2-m and wind at 10-m diagnostic levels:
352    LOGICAL, PARAMETER :: ts_model_level = .FALSE. 
354    !Allocate the arrays for wind components
355 #if ( EM_CORE == 1 )   
356    ALLOCATE ( earth_u_profile(grid%max_ts_level), earth_v_profile(grid%max_ts_level) )
357 #endif
359    IF ( grid%ntsloc_domain .LE. 0 ) RETURN
361 #if ((EM_CORE == 1) && (DA_CORE != 1))
362    IF ( grid%dfi_opt /= DFI_NODFI .AND. grid%dfi_stage /= DFI_FST ) RETURN
363 #endif
365    n = grid%next_ts_time
367    ALLOCATE(p8w(grid%sm32:grid%em32))
368    CALL model_to_grid_config_rec ( grid%id , model_config_rec , config_flags )
370    DO i=1,grid%ntsloc_domain
372       ix = grid%itsloc(i)
373       iy = grid%jtsloc(i)
374   
375       IF (grid%sp31 <= ix .AND. ix <= grid%ep31 .AND. &
376           grid%sp33 <= iy .AND. iy <= grid%ep33) THEN
377        
378          IF (ts_model_level) THEN
379    
380             !
381             ! Output from the lowest model computational level:
382             !
383 #if (EM_CORE == 1)
384             earth_u = grid%u_2(ix,1,iy)*grid%cosa(ix,iy)-grid%v_2(ix,1,iy)*grid%sina(ix,iy)
385             earth_v = grid%v_2(ix,1,iy)*grid%cosa(ix,iy)+grid%u_2(ix,1,iy)*grid%sina(ix,iy)
386             IF (grid%use_theta_m == 1) THEN
387                output_t = (grid%t_2(ix,1,iy)+T0)/(1.+R_v/R_d*grid%moist(ix,1,iy,P_QV)) - T0
388             ELSE
389                output_t =  grid%t_2(ix,1,iy)
390             END IF
391 #else
392             earth_u = grid%u(ix,1,iy)
393             earth_v = grid%v(ix,1,iy)
394             output_t = grid%t(ix,1,iy)
395 #endif
396             output_q = grid%moist(ix,1,iy,P_QV)
397    
398          ELSE
399    
400             !
401             ! Output at 2-m and 10-m diagnostic levels:
402             !
403 #if (EM_CORE == 1)
404             DO k=1,grid%max_ts_level
405                 ! interpolation loop: if you want u,v on cell centers, tslist_unstagger_winds = .true.
406                 IF (config_flags%tslist_unstagger_winds) THEN
407                     earth_u_profile(k) = &
408                     ((grid%u_2(ix,k,iy)*grid%cosa(ix,iy)-grid%v_2(ix,k,iy)*grid%sina(ix,iy)) + &
409                     (grid%u_2(ix+1,k,iy)*grid%cosa(ix+1,iy)-grid%v_2(ix+1,k,iy)*grid%sina(ix+1,iy)))/2.0 
410                     earth_v_profile(k) = &
411                     ((grid%v_2(ix,k,iy)*grid%cosa(ix,iy)+grid%u_2(ix,k,iy)*grid%sina(ix,iy)) + &
412                     (grid%v_2(ix,k,iy+1)*grid%cosa(ix,iy+1)+grid%u_2(ix,k,iy+1)*grid%sina(ix,iy+1)))/2.0 
413                 ELSE
414                     earth_u_profile(k) = grid%u_2(ix,k,iy)*grid%cosa(ix,iy)-grid%v_2(ix,k,iy)*grid%sina(ix,iy) 
415                     earth_v_profile(k) = grid%v_2(ix,k,iy)*grid%cosa(ix,iy)+grid%u_2(ix,k,iy)*grid%sina(ix,iy)
416                 END IF
417             END DO
418             earth_u = grid%u10(ix,iy)*grid%cosa(ix,iy)-grid%v10(ix,iy)*grid%sina(ix,iy)
419             earth_v = grid%v10(ix,iy)*grid%cosa(ix,iy)+grid%u10(ix,iy)*grid%sina(ix,iy)
420             output_q = grid%q2(ix,iy)
421 #else
422             earth_u = grid%u10(ix,iy)
423             earth_v = grid%v10(ix,iy)
424             output_q = grid%qsfc(ix,iy)
425 #endif
426             output_t = grid%t2(ix,iy)
427    
428          END IF
429    
430 #if (EM_CORE == 1)
431          ! Calculate column-integrated liquid/ice  (kg/m^2 or mm)
432          CALL calc_p8w(grid, ix, iy, p8w, grid%sm32, grid%em32)
433          clw=0.
434          DO mm = 1, num_moist
435             IF ( (mm == P_QC) .OR. (mm == P_QR) .OR. (mm == P_QI) .OR. &
436                  (mm == P_QS) .OR. (mm == P_QG) ) THEN
437                DO k=grid%sm32,grid%em32-1
438                   clw=clw+grid%moist(ix,k,iy,mm)*(p8w(k)-p8w(k+1))
439                END DO
440            END IF
441          END DO
442          clw = clw / g
443 #endif
444    
445          CALL domain_clock_get( grid, minutesSinceSimulationStart=xtime_minutes )
446          grid%ts_hour(n,i) = xtime_minutes / 60.
447 #if (EM_CORE == 1)
448             !Create vertical profiles, from lowest model level up to desired level max_ts_level
449             DO k=1,grid%max_ts_level
450             grid%ts_u_profile(n,i,k)   = earth_u_profile(k)
451             grid%ts_v_profile(n,i,k)   = earth_v_profile(k)
452             grid%ts_w_profile(n,i,k)   = (grid%w_2(ix,k,iy)+grid%w_2(ix,k+1,iy))/2.0 ! w on cell center
453             grid%ts_gph_profile(n,i,k) = 0.5*((grid%phb(ix,k,iy)+grid%ph_2(ix,k,iy)) &
454                                              +(grid%phb(ix,k+1,iy)+grid%ph_2(ix,k+1,iy)))/9.81 
455             IF (grid%use_theta_m == 1) THEN
456                grid%ts_th_profile(n,i,k)  = (grid%t_2(ix,k,iy) + T0)/(1.+R_v/R_d*grid%moist(ix,k,iy,P_QV))
457             ELSE
458                grid%ts_th_profile(n,i,k)  =  grid%t_2(ix,k,iy) + T0
459             END IF
460             grid%ts_qv_profile(n,i,k)  = grid%moist(ix,k,iy,P_QV)
461             grid%ts_p_profile(n,i,k)   = grid%pb(ix,k,iy)+grid%p(ix,k,iy)
462             END DO
463 #endif
464 #if (WRF_CHEM == 1)
465             DO k=1,grid%max_ts_level
466             grid%ts_o3_profile(n,i,k)   = grid%chem(ix,k,iy,p_o3)
467             END DO
468 #endif
469          grid%ts_u(n,i)    = earth_u
470          grid%ts_v(n,i)    = earth_v
471          grid%ts_t(n,i)    = output_t
472          grid%ts_q(n,i)    = output_q
473          grid%ts_psfc(n,i) = grid%psfc(ix,iy)
474 #if (EM_CORE == 1)
475          grid%ts_glw(n,i)  = grid%glw(ix,iy)
476          grid%ts_gsw(n,i)  = grid%gsw(ix,iy)
477          grid%ts_hfx(n,i)  = grid%hfx(ix,iy)
478          grid%ts_lh(n,i)   = grid%lh(ix,iy)
479          grid%ts_clw(n,i)  = clw
480          grid%ts_rainc(n,i)  = grid%rainc(ix,iy)
481          grid%ts_rainnc(n,i) = grid%rainnc(ix,iy)
482          grid%ts_tsk(n,i)  = grid%tsk(ix,iy)
483          IF ( model_config_rec%process_time_series == 2 ) THEN
484             !!! Solar diagnostics
485             grid%ts_cldfrac2d(n,i) = grid%cldfrac2d(ix,iy)
486             grid%ts_wvp(n,i) = grid%wvp(ix,iy)
487             grid%ts_lwp(n,i) = grid%lwp(ix,iy)
488             grid%ts_iwp(n,i) = grid%iwp(ix,iy)
489             grid%ts_swp(n,i) = grid%swp(ix,iy)
490             grid%ts_wp_sum(n,i) = grid%wp_sum(ix,iy)
491             grid%ts_lwp_tot(n,i) = grid%lwp_tot(ix,iy)
492             grid%ts_iwp_tot(n,i) = grid%iwp_tot(ix,iy)
493             grid%ts_wp_tot_sum(n,i) = grid%wp_tot_sum(ix,iy)
494             grid%ts_re_qc(n,i) = grid%re_qc(ix,iy)
495             grid%ts_re_qi(n,i) = grid%re_qi(ix,iy)
496             grid%ts_re_qs(n,i) = grid%re_qs(ix,iy)
497             grid%ts_re_qc_tot(n,i) = grid%re_qc_tot(ix,iy)
498             grid%ts_re_qi_tot(n,i) = grid%re_qi_tot(ix,iy)
499             grid%ts_tau_qc(n,i) = grid%tau_qc(ix,iy)
500             grid%ts_tau_qi(n,i) = grid%tau_qi(ix,iy)
501             grid%ts_tau_qs(n,i) = grid%tau_qs(ix,iy)
502             grid%ts_tau_qc_tot(n,i) = grid%tau_qc_tot(ix,iy)
503             grid%ts_tau_qi_tot(n,i) = grid%tau_qi_tot(ix,iy)
504             grid%ts_cbaseht(n,i) = grid%cbaseht(ix,iy)
505             grid%ts_ctopht(n,i) = grid%ctopht(ix,iy)
506             grid%ts_cbaseht_tot(n,i) = grid%cbaseht_tot(ix,iy)
507             grid%ts_ctopht_tot(n,i) = grid%ctopht_tot(ix,iy)
508             grid%ts_clrnidx(n,i) = grid%clrnidx(ix,iy)
509             grid%ts_sza(n,i) = grid%sza(ix,iy)
510             grid%ts_ghi_accum(n,i) = grid%ghi_accum(ix,iy)
511             grid%ts_swdown(n,i) = grid%swdown(ix,iy)
512             grid%ts_swddni(n,i) = grid%swddni(ix,iy)
513             grid%ts_swddif(n,i) = grid%swddif(ix,iy)
514             ! Calc extra solar variables if FARMS or RRTMG/RRTMG FAST
515             if ( config_flags%swint_opt == 2 .or. &
516                  config_flags%ra_sw_physics == RRTMG_SWSCHEME .or. &
517                  config_flags%ra_sw_physics == RRTMG_SWSCHEME_FAST ) then
518               grid%ts_swdownc(n,i) = grid%swdownc(ix,iy)
519               grid%ts_swddnic(n,i) = grid%swddnic(ix,iy)
520               if ( config_flags%swint_opt == 2 ) then  ! FARMS
521                 grid%ts_swdown2(n,i) = grid%swdown2(ix,iy)
522                 grid%ts_swddni2(n,i) = grid%swddni2(ix,iy)
523                 grid%ts_swddif2(n,i) = grid%swddif2(ix,iy)
524                 grid%ts_swdownc2(n,i) = grid%swdownc2(ix,iy)
525                 grid%ts_swddnic2(n,i) = grid%swddnic2(ix,iy)
526               else
527                 grid%ts_swdown2(n,i) = MISSING
528                 grid%ts_swddni2(n,i) = MISSING
529                 grid%ts_swddif2(n,i) = MISSING
530                 grid%ts_swdownc2(n,i) = MISSING
531                 grid%ts_swddnic2(n,i) = MISSING
532               end if
533             else
534               grid%ts_swdownc(n,i) = MISSING
535               grid%ts_swddnic(n,i) = MISSING
536               grid%ts_swdown2(n,i) = MISSING
537               grid%ts_swddni2(n,i) = MISSING
538               grid%ts_swddif2(n,i) = MISSING
539               grid%ts_swdownc2(n,i) = MISSING
540               grid%ts_swddnic2(n,i) = MISSING
541             end if
542          END IF
543 #else
544          grid%ts_tsk(n,i)  = grid%nmm_tsk(ix,iy)
545 #endif
546          grid%ts_tslb(n,i) = grid%tslb(ix,1,iy)
548       ELSE
549 #if (EM_CORE == 1 )   
550          DO k=1,grid%max_ts_level
551          grid%ts_u_profile(n,i,k)     = 1.E30
552          grid%ts_v_profile(n,i,k)     = 1.E30
553          grid%ts_w_profile(n,i,k)     = 1.E30 
554          grid%ts_gph_profile(n,i,k)   = 1.E30
555          grid%ts_th_profile(n,i,k)    = 1.E30
556          grid%ts_qv_profile(n,i,k)    = 1.E30
557          grid%ts_p_profile(n,i,k)     = 1.E30
558          END DO
559 #endif   
560 #if (WRF_CHEM == 1)
561          DO k=1,grid%max_ts_level
562          grid%ts_o3_profile(n,i,k)    = 1.E30
563          END DO
564 #endif
565          grid%ts_hour(n,i) = 1.E30
566          grid%ts_u(n,i)    = 1.E30
567          grid%ts_v(n,i)    = 1.E30
568          grid%ts_t(n,i)    = 1.E30
569          grid%ts_q(n,i)    = 1.E30
570          grid%ts_psfc(n,i) = 1.E30
571 #if (EM_CORE == 1)
572          grid%ts_glw(n,i)  = 1.E30
573          grid%ts_gsw(n,i)  = 1.E30
574          grid%ts_hfx(n,i)  = 1.E30
575          grid%ts_lh(n,i)   = 1.E30
576          grid%ts_clw(n,i)  = 1.E30
577          grid%ts_rainc(n,i)  = 1.E30
578          grid%ts_rainnc(n,i) = 1.E30
579          IF ( model_config_rec%process_time_series == 2 ) THEN
580             !!! Solar diagnostics
581             grid%ts_cldfrac2d(n,i) = 1.E30
582             grid%ts_wvp(n,i) = 1.E30
583             grid%ts_lwp(n,i) = 1.E30
584             grid%ts_iwp(n,i) = 1.E30
585             grid%ts_swp(n,i) = 1.E30
586             grid%ts_wp_sum(n,i) = 1.E30
587             grid%ts_lwp_tot(n,i) = 1.E30
588             grid%ts_iwp_tot(n,i) = 1.E30
589             grid%ts_wp_tot_sum(n,i) = 1.E30
590             grid%ts_re_qc(n,i) = 1.E30
591             grid%ts_re_qi(n,i) = 1.E30
592             grid%ts_re_qs(n,i) = 1.E30
593             grid%ts_re_qc_tot(n,i) = 1.E30
594             grid%ts_re_qi_tot(n,i) = 1.E30
595             grid%ts_tau_qc(n,i) = 1.E30
596             grid%ts_tau_qi(n,i) = 1.E30
597             grid%ts_tau_qs(n,i) = 1.E30
598             grid%ts_tau_qc_tot(n,i) = 1.E30
599             grid%ts_tau_qi_tot(n,i) = 1.E30
600             grid%ts_cbaseht(n,i) = 1.E30
601             grid%ts_ctopht(n,i) = 1.E30
602             grid%ts_cbaseht_tot(n,i) = 1.E30
603             grid%ts_ctopht_tot(n,i) = 1.E30
604             grid%ts_clrnidx(n,i) = 1.E30
605             grid%ts_sza(n,i) = 1.E30
606             grid%ts_ghi_accum(n,i) = 1.E30
607             grid%ts_swdown(n,i) = 1.E30
608             grid%ts_swddni(n,i) = 1.E30
609             grid%ts_swddif(n,i) = 1.E30
610             grid%ts_swdownc(n,i) = 1.E30
611             grid%ts_swddnic(n,i) = 1.E30
612             grid%ts_swdown2(n,i) = 1.E30
613             grid%ts_swddni2(n,i) = 1.E30
614             grid%ts_swddif2(n,i) = 1.E30
615             grid%ts_swdownc2(n,i) = 1.E30
616             grid%ts_swddnic2(n,i) = 1.E30
617          END IF
618 #endif
619          grid%ts_tsk(n,i)  = 1.E30
620          grid%ts_tslb(n,i) = 1.E30
622       END IF
623    END DO
625 #if (EM_CORE == 1) 
626    DEALLOCATE(p8w, earth_u_profile, earth_v_profile)
627 #endif 
629    grid%next_ts_time = grid%next_ts_time + 1
631    IF ( grid%next_ts_time > grid%ts_buf_size ) CALL write_ts(grid)
633 END SUBROUTINE calc_ts
636 SUBROUTINE write_ts( grid )
638    USE module_domain, ONLY : domain
639    USE module_dm, ONLY : wrf_dm_min_reals
640    USE module_state_description
641    USE module_configure, ONLY : model_config_rec
643    IMPLICIT NONE
645    ! Arguments
646    TYPE (domain), INTENT(INOUT) :: grid
648    LOGICAL, EXTERNAL :: wrf_dm_on_monitor
649    INTEGER, EXTERNAL :: get_unused_unit
651    ! Local variables
652    INTEGER :: i, n, ix, iy, iunit, k
653    REAL, ALLOCATABLE, DIMENSION(:,:) :: ts_buf
654    CHARACTER (LEN=24) :: ts_profile_filename
655    CHARACTER (LEN=26) :: profile_format
657    IF ( grid%ntsloc_domain .LE. 0 ) RETURN
659 #if ((EM_CORE == 1) && (DA_CORE != 1))
660    IF ( grid%dfi_opt /= DFI_NODFI .AND. grid%dfi_stage /= DFI_FST ) RETURN
661 #endif
663 #ifdef DM_PARALLEL
664    ALLOCATE(ts_buf(grid%ts_buf_size,grid%max_ts_locs))
666    ts_buf(:,:) = grid%ts_hour(:,:)
667    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_hour(:,:),grid%ts_buf_size*grid%max_ts_locs)
669 #if (EM_CORE == 1) 
670    DO k=1,grid%max_ts_level
671    ts_buf(:,:) = grid%ts_u_profile(:,:,k)
672    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_u_profile(:,:,k),grid%ts_buf_size*grid%max_ts_locs)
673    END DO
675    DO k=1,grid%max_ts_level
676    ts_buf(:,:) = grid%ts_v_profile(:,:,k)
677    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_v_profile(:,:,k),grid%ts_buf_size*grid%max_ts_locs)
678    END DO
680    DO k=1,grid%max_ts_level
681    ts_buf(:,:) = grid%ts_w_profile(:,:,k)
682    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_w_profile(:,:,k),grid%ts_buf_size*grid%max_ts_locs)
683    END DO
685    DO k=1,grid%max_ts_level
686    ts_buf(:,:) = grid%ts_gph_profile(:,:,k)
687    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_gph_profile(:,:,k),grid%ts_buf_size*grid%max_ts_locs)
688    END DO
690    DO k=1,grid%max_ts_level
691    ts_buf(:,:) = grid%ts_th_profile(:,:,k)
692    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_th_profile(:,:,k),grid%ts_buf_size*grid%max_ts_locs)
693    END DO
695    DO k=1,grid%max_ts_level
696    ts_buf(:,:) = grid%ts_qv_profile(:,:,k)
697    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_qv_profile(:,:,k),grid%ts_buf_size*grid%max_ts_locs)
698    END DO
700    DO k=1,grid%max_ts_level
701    ts_buf(:,:) = grid%ts_p_profile(:,:,k)
702    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_p_profile(:,:,k),grid%ts_buf_size*grid%max_ts_locs)
703    END DO
704 #endif         
705 #if (WRF_CHEM == 1)
706    DO k=1,grid%max_ts_level
707    ts_buf(:,:) = grid%ts_o3_profile(:,:,k)
708    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_o3_profile(:,:,k),grid%ts_buf_size*grid%max_ts_locs)
709    END DO
710 #endif
711    ts_buf(:,:) = grid%ts_u(:,:)
712    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_u(:,:),grid%ts_buf_size*grid%max_ts_locs)
714    ts_buf(:,:) = grid%ts_v(:,:)
715    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_v(:,:),grid%ts_buf_size*grid%max_ts_locs)
717    ts_buf(:,:) = grid%ts_t(:,:)
718    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_t(:,:),grid%ts_buf_size*grid%max_ts_locs)
720    ts_buf(:,:) = grid%ts_q(:,:)
721    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_q(:,:),grid%ts_buf_size*grid%max_ts_locs)
723    ts_buf(:,:) = grid%ts_psfc(:,:)
724    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_psfc(:,:),grid%ts_buf_size*grid%max_ts_locs)
726 #if (EM_CORE == 1)
727    ts_buf(:,:) = grid%ts_glw(:,:)
728    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_glw(:,:),grid%ts_buf_size*grid%max_ts_locs)
730    ts_buf(:,:) = grid%ts_gsw(:,:)
731    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_gsw(:,:),grid%ts_buf_size*grid%max_ts_locs)
733    ts_buf(:,:) = grid%ts_hfx(:,:)
734    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_hfx(:,:),grid%ts_buf_size*grid%max_ts_locs)
736    ts_buf(:,:) = grid%ts_lh(:,:)
737    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_lh(:,:),grid%ts_buf_size*grid%max_ts_locs)
739    ts_buf(:,:) = grid%ts_clw(:,:)
740    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_clw(:,:),grid%ts_buf_size*grid%max_ts_locs)
742    ts_buf(:,:) = grid%ts_rainc(:,:)
743    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_rainc(:,:),grid%ts_buf_size*grid%max_ts_locs)
745    ts_buf(:,:) = grid%ts_rainnc(:,:)
746    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_rainnc(:,:),grid%ts_buf_size*grid%max_ts_locs)
748    IF ( model_config_rec%process_time_series == 2 ) THEN
749       !!! Solar diagnostics
750       ts_buf(:,:) = grid%ts_cldfrac2d(:,:)
751       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_cldfrac2d(:,:),grid%ts_buf_size*grid%max_ts_locs)
753       ts_buf(:,:) = grid%ts_wvp(:,:)
754       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_wvp(:,:),grid%ts_buf_size*grid%max_ts_locs)
756       ts_buf(:,:) = grid%ts_lwp(:,:)
757       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_lwp(:,:),grid%ts_buf_size*grid%max_ts_locs)
759       ts_buf(:,:) = grid%ts_iwp(:,:)
760       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_iwp(:,:),grid%ts_buf_size*grid%max_ts_locs)
762       ts_buf(:,:) = grid%ts_swp(:,:)
763       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swp(:,:),grid%ts_buf_size*grid%max_ts_locs)
765       ts_buf(:,:) = grid%ts_wp_sum(:,:)
766       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_wp_sum(:,:),grid%ts_buf_size*grid%max_ts_locs)
768       ts_buf(:,:) = grid%ts_lwp_tot(:,:)
769       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_lwp_tot(:,:),grid%ts_buf_size*grid%max_ts_locs)
771       ts_buf(:,:) = grid%ts_iwp_tot(:,:)
772       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_iwp_tot(:,:),grid%ts_buf_size*grid%max_ts_locs)
774       ts_buf(:,:) = grid%ts_wp_tot_sum(:,:)
775       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_wp_tot_sum(:,:),grid%ts_buf_size*grid%max_ts_locs)
777       ts_buf(:,:) = grid%ts_re_qc(:,:)
778       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_re_qc(:,:),grid%ts_buf_size*grid%max_ts_locs)
780       ts_buf(:,:) = grid%ts_re_qi(:,:)
781       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_re_qi(:,:),grid%ts_buf_size*grid%max_ts_locs)
783       ts_buf(:,:) = grid%ts_re_qs(:,:)
784       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_re_qs(:,:),grid%ts_buf_size*grid%max_ts_locs)
786       ts_buf(:,:) = grid%ts_re_qc_tot(:,:)
787       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_re_qc_tot(:,:),grid%ts_buf_size*grid%max_ts_locs)
789       ts_buf(:,:) = grid%ts_re_qi_tot(:,:)
790       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_re_qi_tot(:,:),grid%ts_buf_size*grid%max_ts_locs)
792       ts_buf(:,:) = grid%ts_tau_qc(:,:)
793       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_tau_qc(:,:),grid%ts_buf_size*grid%max_ts_locs)
795       ts_buf(:,:) = grid%ts_tau_qi(:,:)
796       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_tau_qi(:,:),grid%ts_buf_size*grid%max_ts_locs)
798       ts_buf(:,:) = grid%ts_tau_qs(:,:)
799       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_tau_qs(:,:),grid%ts_buf_size*grid%max_ts_locs)
801       ts_buf(:,:) = grid%ts_tau_qc_tot(:,:)
802       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_tau_qc_tot(:,:),grid%ts_buf_size*grid%max_ts_locs)
804       ts_buf(:,:) = grid%ts_tau_qi_tot(:,:)
805       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_tau_qi_tot(:,:),grid%ts_buf_size*grid%max_ts_locs)
807       ts_buf(:,:) = grid%ts_cbaseht(:,:)
808       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_cbaseht(:,:),grid%ts_buf_size*grid%max_ts_locs)
810       ts_buf(:,:) = grid%ts_ctopht(:,:)
811       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_ctopht(:,:),grid%ts_buf_size*grid%max_ts_locs)
813       ts_buf(:,:) = grid%ts_cbaseht_tot(:,:)
814       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_cbaseht_tot(:,:),grid%ts_buf_size*grid%max_ts_locs)
816       ts_buf(:,:) = grid%ts_ctopht_tot(:,:)
817       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_ctopht_tot(:,:),grid%ts_buf_size*grid%max_ts_locs)
819       ts_buf(:,:) = grid%ts_clrnidx(:,:)
820       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_clrnidx(:,:),grid%ts_buf_size*grid%max_ts_locs)
822       ts_buf(:,:) = grid%ts_sza(:,:)
823       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_sza(:,:),grid%ts_buf_size*grid%max_ts_locs)
825       ts_buf(:,:) = grid%ts_ghi_accum(:,:)
826       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_ghi_accum(:,:),grid%ts_buf_size*grid%max_ts_locs)
828       ts_buf(:,:) = grid%ts_swdown(:,:)
829       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swdown(:,:),grid%ts_buf_size*grid%max_ts_locs)
831       ts_buf(:,:) = grid%ts_swddni(:,:)
832       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swddni(:,:),grid%ts_buf_size*grid%max_ts_locs)
834       ts_buf(:,:) = grid%ts_swddif(:,:)
835       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swddif(:,:),grid%ts_buf_size*grid%max_ts_locs)
837       ts_buf(:,:) = grid%ts_swdownc(:,:)
838       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swdownc(:,:),grid%ts_buf_size*grid%max_ts_locs)
840       ts_buf(:,:) = grid%ts_swddnic(:,:)
841       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swddnic(:,:),grid%ts_buf_size*grid%max_ts_locs)
843       ts_buf(:,:) = grid%ts_swdown2(:,:)
844       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swdown2(:,:),grid%ts_buf_size*grid%max_ts_locs)
846       ts_buf(:,:) = grid%ts_swddni2(:,:)
847       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swddni2(:,:),grid%ts_buf_size*grid%max_ts_locs)
849       ts_buf(:,:) = grid%ts_swddif2(:,:)
850       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swddif2(:,:),grid%ts_buf_size*grid%max_ts_locs)
852       ts_buf(:,:) = grid%ts_swdownc2(:,:)
853       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swdownc2(:,:),grid%ts_buf_size*grid%max_ts_locs)
855       ts_buf(:,:) = grid%ts_swddnic2(:,:)
856       CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_swddnic2(:,:),grid%ts_buf_size*grid%max_ts_locs)
857    END IF
858 #endif
860    ts_buf(:,:) = grid%ts_tsk(:,:)
861    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_tsk(:,:),grid%ts_buf_size*grid%max_ts_locs)
863    ts_buf(:,:) = grid%ts_tslb(:,:)
864    CALL wrf_dm_min_reals(ts_buf(:,:),grid%ts_tslb(:,:),grid%ts_buf_size*grid%max_ts_locs)
866    DEALLOCATE(ts_buf)
867 #endif
869    IF ( wrf_dm_on_monitor() ) THEN
871       iunit = get_unused_unit()
872       IF ( iunit <= 0 ) THEN
873          CALL wrf_error_fatal('Error in write_ts: could not find a free Fortran unit.')
874       END IF
876       DO i=1,grid%ntsloc_domain
878          ix = grid%itsloc(i)
879          iy = grid%jtsloc(i)
881          OPEN(UNIT=iunit, FILE=TRIM(grid%ts_filename(i)), STATUS='unknown', POSITION='append', FORM='formatted')
883          DO n=1,grid%next_ts_time - 1
885 #if (EM_CORE == 1)
886             IF ( model_config_rec%process_time_series == 1 ) THEN
887                WRITE(UNIT=iunit,FMT='(i2,f13.6,i5,i5,i5,1x,14(f13.5,1x))')  &
888                                  grid%id, grid%ts_hour(n,i),                &
889                                  grid%id_tsloc(i), ix, iy,                  &
890                                  grid%ts_t(n,i),                            &
891                                  grid%ts_q(n,i),                            &
892                                  grid%ts_u(n,i),                            &
893                                  grid%ts_v(n,i),                            &
894                                  grid%ts_psfc(n,i),                         &
895                                  grid%ts_glw(n,i),                          &
896                                  grid%ts_gsw(n,i),                          &
897                                  grid%ts_hfx(n,i),                          &
898                                  grid%ts_lh(n,i),                           &
899                                  grid%ts_tsk(n,i),                          &
900                                  grid%ts_tslb(n,i),                         &
901                                  grid%ts_rainc(n,i),                        &
902                                  grid%ts_rainnc(n,i),                       &
903                                  grid%ts_clw(n,i)
904             ELSE
905                !!! WRF-Solar diagnostics
906                WRITE(UNIT=iunit,FMT='(i2,f13.6,i5,i5,i5,1x,49(f13.5,1x))')  &
907                                  grid%id, grid%ts_hour(n,i),                &
908                                  grid%id_tsloc(i), ix, iy,                  &
909                                  grid%ts_t(n,i),                            &
910                                  grid%ts_q(n,i),                            &
911                                  grid%ts_u(n,i),                            &
912                                  grid%ts_v(n,i),                            &
913                                  grid%ts_psfc(n,i),                         &
914                                  grid%ts_glw(n,i),                          &
915                                  grid%ts_gsw(n,i),                          &
916                                  grid%ts_hfx(n,i),                          &
917                                  grid%ts_lh(n,i),                           &
918                                  grid%ts_tsk(n,i),                          &
919                                  grid%ts_tslb(n,i),                         &
920                                  grid%ts_rainc(n,i),                        &
921                                  grid%ts_rainnc(n,i),                       &
922                                  grid%ts_clw(n,i),                          &
923                                  grid%ts_cldfrac2d(n,i),                    &
924                                  grid%ts_wvp(n,i),                          &
925                                  grid%ts_lwp(n,i),                          &
926                                  grid%ts_iwp(n,i),                          &
927                                  grid%ts_swp(n,i),                          &
928                                  grid%ts_wp_sum(n,i),                       &
929                                  grid%ts_lwp_tot(n,i),                      &
930                                  grid%ts_iwp_tot(n,i),                      &
931                                  grid%ts_wp_tot_sum(n,i),                   &
932                                  grid%ts_re_qc(n,i),                        &
933                                  grid%ts_re_qi(n,i),                        &
934                                  grid%ts_re_qs(n,i),                        &
935                                  grid%ts_re_qc_tot(n,i),                    &
936                                  grid%ts_re_qi_tot(n,i),                    &
937                                  grid%ts_tau_qc(n,i),                       &
938                                  grid%ts_tau_qi(n,i),                       &
939                                  grid%ts_tau_qs(n,i),                       &
940                                  grid%ts_tau_qc_tot(n,i),                   &
941                                  grid%ts_tau_qi_tot(n,i),                   &
942                                  grid%ts_cbaseht(n,i),                      &
943                                  grid%ts_ctopht(n,i),                       &
944                                  grid%ts_cbaseht_tot(n,i),                  &
945                                  grid%ts_ctopht_tot(n,i),                   &
946                                  grid%ts_clrnidx(n,i),                      &
947                                  grid%ts_sza(n,i),                          &
948                                  grid%ts_ghi_accum(n,i),                    &
949                                  grid%ts_swdown(n,i),                       &
950                                  grid%ts_swddni(n,i),                       &
951                                  grid%ts_swddif(n,i),                       &
952                                  grid%ts_swdownc(n,i),                      &
953                                  grid%ts_swddnic(n,i),                      &
954                                  grid%ts_swdown2(n,i),                      &
955                                  grid%ts_swddni2(n,i),                      &
956                                  grid%ts_swddif2(n,i),                      &
957                                  grid%ts_swdownc2(n,i),                     &
958                                  grid%ts_swddnic2(n,i)
959             END IF
960 #else
961             WRITE(UNIT=iunit,FMT='(i2,f13.6,i5,i5,i5,1x,7(f13.5,1x))')   &
962                               grid%id, grid%ts_hour(n,i),                &
963                               grid%id_tsloc(i), ix, iy,                  &
964                               grid%ts_t(n,i),                            &
965                               grid%ts_q(n,i),                            &
966                               grid%ts_u(n,i),                            &
967                               grid%ts_v(n,i),                            &
968                               grid%ts_psfc(n,i),                         &
969                               grid%ts_tsk(n,i),                          &
970                               grid%ts_tslb(n,i)
971 #endif
972          END DO
973          CLOSE(UNIT=iunit)
974          
975          !Set write format for vertical profiles, depending on the highest model level of interest
976 #if (EM_CORE == 1)   
977          profile_format = '(f13.6,1x,000(f13.5,1x))'
978          k= LEN_TRIM(profile_format)
979          WRITE(profile_format(12:14),'(I3.3)') grid%max_ts_level
981          !Write u wind component profile to separate file
982          iunit = get_unused_unit()
983             IF ( iunit <= 0 ) THEN
984             CALL wrf_error_fatal('Error in write_ts: could not find a free Fortran unit.')
985             END IF
986          !Recreate filename for u wind component profiles
987          WRITE(ts_profile_filename,'(A)') TRIM(grid%nametsloc(grid%id_tsloc(i)))//'.d00.TS'
988          k = LEN_TRIM(ts_profile_filename)
989          WRITE(ts_profile_filename(k-4:k-3),'(I2.2)') grid%id
990          WRITE(ts_profile_filename(k-1:k),'(A2)') 'UU'
991          
992          OPEN(UNIT=iunit, FILE=TRIM(ts_profile_filename), STATUS='unknown', POSITION='append', FORM='formatted')
994          DO n=1,grid%next_ts_time - 1
996             WRITE(UNIT=iunit,FMT=profile_format)           &
997                               grid%ts_hour(n,i),                      &
998                               grid%ts_u_profile(n,i,1:grid%max_ts_level)                                       
999          END DO
1000          CLOSE(UNIT=iunit)
1002          !Write v wind component profile to separate file
1003          iunit = get_unused_unit()
1004             IF ( iunit <= 0 ) THEN
1005             CALL wrf_error_fatal('Error in write_ts: could not find a free Fortran unit.')
1006             END IF
1007          !Recreate filename for v wind component profiles
1008          k = LEN_TRIM(ts_profile_filename)
1009          WRITE(ts_profile_filename(k-1:k),'(A2)') 'VV'
1010          
1011          OPEN(UNIT=iunit, FILE=TRIM(ts_profile_filename), STATUS='unknown', POSITION='append', FORM='formatted')
1013          DO n=1,grid%next_ts_time - 1
1015             WRITE(UNIT=iunit,FMT=profile_format)  &
1016                               grid%ts_hour(n,i),             & 
1017                               grid%ts_v_profile(n,i,1:grid%max_ts_level)                                       
1018          END DO
1019          CLOSE(UNIT=iunit)
1021          iunit = get_unused_unit()
1022             IF ( iunit <= 0 ) THEN
1023             CALL wrf_error_fatal('Error in write_ts: could not find a free Fortran unit.')
1024             END IF
1025          k = LEN_TRIM(ts_profile_filename)
1026          WRITE(ts_profile_filename(k-1:k),'(A2)') 'WW'
1027          
1028          OPEN(UNIT=iunit, FILE=TRIM(ts_profile_filename), STATUS='unknown', POSITION='append', FORM='formatted')
1030          DO n=1,grid%next_ts_time - 1
1032             WRITE(UNIT=iunit,FMT=profile_format)  &
1033                               grid%ts_hour(n,i),             & 
1034                               grid%ts_w_profile(n,i,1:grid%max_ts_level)                                       
1035          END DO
1036          CLOSE(UNIT=iunit)
1038          !Write geopotential height profile to separate file
1039          iunit = get_unused_unit()
1040             IF ( iunit <= 0 ) THEN
1041             CALL wrf_error_fatal('Error in write_ts: could not find a free Fortran unit.')
1042             END IF
1043          !Recreate filename for geopotential height profiles
1044          k = LEN_TRIM(ts_profile_filename)
1045          WRITE(ts_profile_filename(k-1:k),'(A2)') 'PH'
1046          
1047          OPEN(UNIT=iunit, FILE=TRIM(ts_profile_filename), STATUS='unknown', POSITION='append', FORM='formatted')
1049          DO n=1,grid%next_ts_time - 1
1051             WRITE(UNIT=iunit,FMT=profile_format)  &
1052                               grid%ts_hour(n,i),             &
1053                               grid%ts_gph_profile(n,i,1:grid%max_ts_level)                                       
1054          END DO
1055          CLOSE(UNIT=iunit)
1056          
1057          !Write potential temperature profile to separate file
1058          iunit = get_unused_unit()
1059             IF ( iunit <= 0 ) THEN
1060             CALL wrf_error_fatal('Error in write_ts: could not find a free Fortran unit.')
1061             END IF
1062          !Recreate filename for potential temperature profiles
1063          k = LEN_TRIM(ts_profile_filename)
1064          WRITE(ts_profile_filename(k-1:k),'(A2)') 'TH'
1065          
1066          OPEN(UNIT=iunit, FILE=TRIM(ts_profile_filename), STATUS='unknown', POSITION='append', FORM='formatted')
1068          DO n=1,grid%next_ts_time - 1
1070             WRITE(UNIT=iunit,FMT=profile_format)   &
1071                               grid%ts_hour(n,i),              &
1072                               grid%ts_th_profile(n,i,1:grid%max_ts_level)                                       
1073          END DO
1074          CLOSE(UNIT=iunit) 
1075        
1076          !Write water vapor mixing ratio profile to separate file
1077          iunit = get_unused_unit()
1078             IF ( iunit <= 0 ) THEN
1079             CALL wrf_error_fatal('Error in write_ts: could not find a free Fortran unit.')
1080             END IF
1081          !Recreate filename for water vapor mixing ratio profiles
1082          k = LEN_TRIM(ts_profile_filename)
1083          WRITE(ts_profile_filename(k-1:k),'(A2)') 'QV'
1084          
1085          OPEN(UNIT=iunit, FILE=TRIM(ts_profile_filename), STATUS='unknown', POSITION='append', FORM='formatted')
1087          DO n=1,grid%next_ts_time - 1
1089             WRITE(UNIT=iunit,FMT=profile_format)  &
1090                               grid%ts_hour(n,i),             &
1091                               grid%ts_qv_profile(n,i,1:grid%max_ts_level)                                       
1092          END DO
1093          CLOSE(UNIT=iunit)
1094        
1095          !Write pressure profile to separate file
1096          iunit = get_unused_unit()
1097             IF ( iunit <= 0 ) THEN
1098             CALL wrf_error_fatal('Error in write_ts: could not find a free Fortran unit.')
1099             END IF
1100          !Recreate filename for pressure profiles
1101          k = LEN_TRIM(ts_profile_filename)
1102          WRITE(ts_profile_filename(k-1:k),'(A2)') 'PR'
1103          
1104          OPEN(UNIT=iunit, FILE=TRIM(ts_profile_filename), STATUS='unknown', POSITION='append', FORM='formatted')
1106          DO n=1,grid%next_ts_time - 1
1108             WRITE(UNIT=iunit,FMT=profile_format)  &
1109                               grid%ts_hour(n,i),             &
1110                               grid%ts_p_profile(n,i,1:grid%max_ts_level)                                       
1111          END DO
1112          CLOSE(UNIT=iunit)
1114 #endif   
1115 #if (WRF_CHEM == 1)
1116          !Write O3 profile to separate file
1117          iunit = get_unused_unit()
1118             IF ( iunit <= 0 ) THEN
1119             CALL wrf_error_fatal('Error in write_ts: could not find a free Fortran unit.')
1120             END IF
1121          !Recreate filename for O3 profiles
1122          k = LEN_TRIM(ts_profile_filename)
1123          WRITE(ts_profile_filename(k-1:k),'(A2)') 'O3'
1125          OPEN(UNIT=iunit, FILE=TRIM(ts_profile_filename), STATUS='unknown', POSITION='append', FORM='formatted')
1127          DO n=1,grid%next_ts_time - 1
1129             WRITE(UNIT=iunit,FMT=profile_format)  &
1130                               grid%ts_hour(n,i),             &
1131                               grid%ts_o3_profile(n,i,1:grid%max_ts_level)
1132          END DO
1133          CLOSE(UNIT=iunit)
1134 #endif
1136       END DO
1138    END IF
1140    grid%next_ts_time = 1
1142 END SUBROUTINE write_ts
1145 #if (EM_CORE == 1)
1146 SUBROUTINE calc_p8w(grid, ix, iy, p8w, k_start, k_end)
1148    USE module_domain
1149    USE module_model_constants
1151    IMPLICIT NONE
1153    ! Arguments
1154    TYPE (domain), INTENT(IN) :: grid
1155    INTEGER, INTENT(IN) :: ix, iy, k_start, k_end
1156    REAL, DIMENSION(k_start:k_end), INTENT(OUT) :: p8w
1158    ! Local variables
1159    INTEGER :: k
1160    REAL    :: z0, z1, z2, w1, w2 
1161    REAL, DIMENSION(k_start:k_end)   :: z_at_w
1162    REAL, DIMENSION(k_start:k_end-1) :: z
1165    DO k = k_start, k_end
1166       z_at_w(k) = (grid%phb(ix,k,iy)+grid%ph_2(ix,k,iy))/g
1167    END DO
1169    DO k = k_start, k_end-1
1170       z(k) = 0.5*(z_at_w(k) + z_at_w(k+1))
1171    END DO
1173    DO k = k_start+1, k_end-1
1174       p8w(k) = grid%fnm(k)*(grid%p(ix,k,iy)+grid%pb(ix,k,iy)) + &
1175                grid%fnp(k)*(grid%p(ix,k-1,iy)+grid%pb(ix,k-1,iy))
1176    END DO
1178    z0 = z_at_w(k_start)
1179    z1 = z(k_start)
1180    z2 = z(k_start+1)
1181    w1 = (z0 - z2)/(z1 - z2)
1182    w2 = 1. - w1
1183    p8w(k_start) = w1*(grid%p(ix,k_start,iy)+grid%pb(ix,k_start,iy)) + &
1184                   w2*(grid%p(ix,k_start+1,iy)+grid%pb(ix,k_start+1,iy))
1186    z0 = z_at_w(k_end)
1187    z1 = z(k_end-1)
1188    z2 = z(k_end-2)
1189    w1 = (z0 - z2)/(z1 - z2)
1190    w2 = 1. - w1
1191    p8w(k_end) = exp(w1*log(grid%p(ix,k_end-1,iy)+grid%pb(ix,k_end-1,iy)) + &
1192                     w2*log(grid%p(ix,k_end-2,iy)+grid%pb(ix,k_end-2,iy)))
1194 END SUBROUTINE calc_p8w
1195 #endif