-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathsibdrv_read_ecmwf.F90
230 lines (194 loc) · 9.73 KB
/
sibdrv_read_ecmwf.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
subroutine sibdrv_read_ecmwf( sib, time )
!****--------------------------------------------------------------------
! This subroutines reads the forcing surface meteorological data
! for the next time step.
! If required, it closes the current month's data file and opens the
! next month's data file.
! precip, snowfall and radiation data are aggregated data over 6 hours
! They have been stored such that the data at time x represent the
! aggregation over the preceeding 6 hours.
! Thus data at
! 0 hours are the sum of the values 18-24 hours of the previous day
! 6 hours are the sum of the values 0- 6 hours of the current day
! 12 hours are the sum of the values 6-12 hours of the current day
! 18 hours are the sum of the values 12-18 hours of the current day
!
! Consequently new data are read for all variables at the same time.
! For all point data variables, e.g. temp., the data are for
! now + 6 hours. For precip and radiation data the new data
! are the aggregation over the next six hours.
!
! Modifications:
! Kevin Schaefer moved conversion from pascals to millibars from sibdrv_interp to here (8/16/04)
! Kevin Schaefer moved conversion from dewpt to humidity here from sibdrv_interp (8/17/04)
!****--------------------------------------------------------------------
use sib_const_module, only: &
nsib, &
latsib, &
lonsib, &
subset, &
subcount
use sib_io_module, only: &
dr_format, &
driver_id
use physical_parameters, only: &
kapa => kappa, &
pi
use kinds
use sibtype
use timetype
use netcdf
use typeSizes
#include "nc_util.h"
type(sib_t), dimension(subcount), intent(inout) :: sib
type(time_struct), intent(in) :: time
integer(kind=int_kind) :: i, iyear, imon, iday, idoy, ihour, imin
character*80 filename
integer(kind=int_kind) :: nctimeid,ncyid,ncmid,nctdid,ncdoyid,nchid
integer(kind=int_kind), dimension(2) :: mstart,mcount
integer(kind=int_kind) :: nct2mid ! Total Cloud Cover
integer(kind=int_kind) :: nctccid ! Total Cloud Cover
integer(kind=int_kind) :: ncswdid ! Surface solar radiation downwards
integer(kind=int_kind) :: ncldwid ! Surface thermal radiation downwards
integer(kind=int_kind) :: ncuwdid ! U-wind at 10 m
integer(kind=int_kind) :: ncvwdid ! V-wind at 10 m
integer(kind=int_kind) :: ncdptid ! Dewpoint at 2 m
integer(kind=int_kind) :: ncsfpid ! Log Surface Pressure
integer(kind=int_kind) :: nclspid ! Large Scale Precipitation
integer(kind=int_kind) :: nccvpid ! Convective Precipitation
integer(kind=int_kind) :: ncsflid ! Snow Fall
real(kind=real_kind), dimension(nsib) :: t2m ! Total Cloud Cover
real(kind=real_kind), dimension(nsib) :: tcc ! Total Cloud Cover
real(kind=real_kind), dimension(nsib) :: swd ! Surface solar radiation downwards
real(kind=real_kind), dimension(nsib) :: ldw ! Surface thermal radiation downwards
real(kind=real_kind), dimension(nsib) :: dpt ! Dewpoint at 2 m
real(kind=dbl_kind), dimension(subcount) :: temp_dpt ! dew point
real(kind=real_kind), dimension(nsib) :: sfp ! Log Surface Pressure
real(kind=real_kind), dimension(nsib) :: lsp ! Large Scale Precipitation
real(kind=real_kind), dimension(nsib) :: cvp ! Convective Precipitation
integer(kind=int_kind) :: status
real(kind=real_kind) :: xtime,xyear,xmonth,xdoy,xday,xhour
real(kind=real_kind), dimension(nsib) :: uwd,vwd
real(kind=real_kind), dimension(nsib) :: xx
character(len=13) :: subname
data subname/'sibdrv_read '/
!*** Storing previous time steps data
do i=1,subcount
sib(i)%prog%ps1 = sib(i)%prog%ps2
sib(i)%prog%tm1 = sib(i)%prog%tm2
sib(i)%prog%tcc1 = sib(i)%prog%tcc2
sib(i)%prog%sh1 = sib(i)%prog%sh2
sib(i)%prog%spdm1 = sib(i)%prog%spdm2
sib(i)%prog%lspr1 = sib(i)%prog%lspr2
sib(i)%prog%cupr1 = sib(i)%prog%cupr2
sib(i)%prog%dlwbot1 = sib(i)%prog%dlwbot2
sib(i)%prog%sw_dwn1 = sib(i)%prog%sw_dwn2
enddo
! switch files if needed
if ( time%switch_driver ) then
status = nf90_close( driver_id )
write( filename, dr_format ) time%driver_year, time%driver_month
CHECK( nf90_open( trim(filename), nf90_nowrite, driver_id ) )
print *, 'switch driver to ',trim(filename)
endif
! read new driver data
! print*, subname,tau,nextsecond,nextday,nextdoy, &
! nmonth,nyear,nextmonth,nextyear
! check time values in driver data file
ENSURE_VAR( driver_id, 'time', nctimeid )
ENSURE_VAR( driver_id, 'year', ncyid )
ENSURE_VAR( driver_id, 'month',ncmid )
ENSURE_VAR( driver_id, 'doy', ncdoyid )
ENSURE_VAR( driver_id, 'day', nctdid )
ENSURE_VAR( driver_id, 'hour', nchid )
print*,subname,nctimeid,ncyid,ncmid,ncdoyid,nctdid,nchid,time%driver_recnum
! read time
mstart(1) = time%driver_recnum
CHECK( nf90_get_var( driver_id, nctimeid, xtime, mstart(1:1) ) )
CHECK( nf90_get_var( driver_id, ncyid, xyear, mstart(1:1) ) )
CHECK( nf90_get_var( driver_id, ncmid, xmonth, mstart(1:1) ) )
CHECK( nf90_get_var( driver_id, ncdoyid, xdoy, mstart(1:1) ) )
CHECK( nf90_get_var( driver_id, nctdid, xday, mstart(1:1) ) )
CHECK( nf90_get_var( driver_id, nchid, xhour, mstart(1:1) ) )
ihour=xhour
iday =xday
idoy =xdoy
imon =xmonth
iyear=xyear
imin=0
print*,subname,'Time level in file: ',ihour,iday,idoy,imon,iyear
if( time%driver_day /= iday .or. time%driver_hour /= ihour) then
! print*,subname,nyear0,nmonthofyear,ndayofyear,nsecondofday/3600
! print*,subname,nextsecond,nexthour,nextday
! print*,subname,'Not the right data'
! print *, time%driver_day, iday, time%driver_hour, ihour
! stop
endif
! get variable id's
ENSURE_VAR( driver_id, 't2m', nct2mid ) ! Temperature at 2 m
ENSURE_VAR( driver_id, 'tcc', nctccid ) ! Total Cloud Cover
ENSURE_VAR( driver_id, 'swd', ncswdid ) ! Surface solar rad downwards
ENSURE_VAR( driver_id, 'ldw', ncldwid ) ! Surface thermal rad down
ENSURE_VAR( driver_id, 'uwd', ncuwdid ) ! U-wind at 10 m
ENSURE_VAR( driver_id, 'vwd', ncvwdid ) ! V-wind at 10 m
ENSURE_VAR( driver_id, 'dpt', ncdptid ) ! Dewpoint at 2 m
ENSURE_VAR( driver_id, 'sfp', ncsfpid ) ! Log Surface Pressure
ENSURE_VAR( driver_id, 'lsp', nclspid ) ! Large Scale Precipitation
ENSURE_VAR( driver_id, 'cvp', nccvpid ) ! Convective Precipitation
ENSURE_VAR( driver_id, 'sfl', ncsflid ) ! Snow Fall
! get data
mstart=(/1,time%driver_recnum/); mcount=(/nsib,1/)
CHECK( nf90_get_var( driver_id, nct2mid, t2m, mstart, mcount ) ) !Temperature at 2 m
CHECK( nf90_get_var( driver_id, nctccid, tcc, mstart, mcount ) ) !Total Cloud Cover
CHECK( nf90_get_var( driver_id, ncswdid, swd, mstart, mcount ) ) !Surface solar rad downwards
CHECK( nf90_get_var( driver_id, ncldwid, ldw, mstart, mcount ) ) !Surface thermal rad downwards
CHECK( nf90_get_var( driver_id, ncuwdid, uwd, mstart, mcount ) ) ! U-wind at 10 m
CHECK( nf90_get_var( driver_id, ncvwdid, vwd, mstart, mcount ) ) ! V-wind at 10 m
CHECK( nf90_get_var( driver_id, ncdptid, dpt, mstart, mcount ) ) ! Dewpoint at 2 m
CHECK( nf90_get_var( driver_id, ncsfpid, sfp, mstart, mcount ) ) ! Surface Pressure
CHECK( nf90_get_var( driver_id, nclspid, lsp, mstart, mcount ) ) ! Large Scale Precipitation
CHECK( nf90_get_var( driver_id, nccvpid, cvp, mstart, mcount ) ) ! Convective Precipitation
CHECK( nf90_get_var(driver_id, ncsflid, xx, mstart, mcount ) ) ! Snow Fall
do i=1,subcount
! pull out landpoints in subdomain
sib(i)%prog%tm2 = t2m(subset(i))
sib(i)%prog%tcc2 = tcc(subset(i))
sib(i)%prog%sw_dwn2 = swd(subset(i))
sib(i)%prog%dlwbot2 = ldw(subset(i))
temp_dpt(i) = dpt(subset(i))
sib(i)%prog%ps2 = sfp(subset(i))
sib(i)%prog%lspr2 = lsp(subset(i))
sib(i)%prog%cupr2 = cvp(subset(i))
! scale radiation to w/m2
sib(i)%prog%sw_dwn2 = sib(i)%prog%sw_dwn2/time%driver_step
sib(i)%prog%dlwbot2 = sib(i)%prog%dlwbot2/time%driver_step
! 10 m wind
sib(i)%prog%spdm2=SQRT(uwd(subset(i))*uwd(subset(i))+vwd(subset(i))*vwd(subset(i)))
! add snowfall to large scale precip and let SiB decide about snow.
sib(i)%prog%lspr2 = sib(i)%prog%lspr2+xx(subset(i))
! convert to mm
sib(i)%prog%lspr2 = sib(i)%prog%lspr2*1000
if ( sib(i)%prog%lspr2 <= 1.e-3 ) sib(i)%prog%lspr2 = 0.0_dbl_kind
sib(i)%prog%cupr2 = sib(i)%prog%cupr2*1000
! KS Conversion ln(pascals) to millibars
sib(i)%prog%ps2 = exp(sib(i)%prog%ps2)*0.01
enddo
!
! KS convert dew point to specific humidity
call qsat_eau(subcount,sib%prog%ps2*100.0,temp_dpt,sib%prog%sh2)
print*,subname,'New driver data read ',ihour,iday,imon,iyear
print*,'------------------------------------------------------'
print*,'Extrema of new input data'
print*, minval(sib%prog%tm2 ),maxval(sib%prog%tm2 ),' Temperature'
print*, minval(sib%prog%tcc2 ),maxval(sib%prog%tcc2),' Total cloudiness'
print*, minval(sib%prog%sh2 ),maxval(sib%prog%sh2 ),' humidity'
print*, minval(sib%prog%spdm2 ),maxval(sib%prog%spdm2),' Surface wind'
print*, minval(sib%prog%ps2 ),maxval(sib%prog%ps2 ),' Pressure'
print*, minval(sib%prog%dlwbot2 ),maxval(sib%prog%dlwbot2), &
' Long wave down'
print*, minval(sib%prog%lspr2 ),maxval(sib%prog%lspr2),' Large sc precip'
print*, minval(sib%prog%cupr2 ),maxval(sib%prog%cupr2),' Convective '
print*, minval(sib%prog%sw_dwn2),maxval(sib%prog%sw_dwn2), &
' Short wave down'
print*,'-----------------------------------------------------'
end subroutine sibdrv_read_ecmwf