Line data Source code
1 : module scamMod
2 : !----------------------------------------------------------------------
3 : !
4 : ! This module provides control variables and namelist functionality
5 : ! for SCAM.
6 : !
7 : ! As with global CAM, SCAM is initialized with state information
8 : ! from the initial and boundary files. For each succeeding timestep
9 : ! SCAM calculates the physics tendencies and combines these with
10 : ! the advective tendencies provided by the IOP file to produce
11 : ! a forcast. Generally, the control variables in this module
12 : ! determine what data and parameterizations will be used to make
13 : ! the forecast. For instance, many of the control variables in
14 : ! this module provide flexibility to affect the forecast by overriding
15 : ! parameterization prognosed tendencies with observed tendencies
16 : ! of a particular field program recorded on the IOP file.
17 : !
18 : ! Public functions/subroutines:
19 : ! scam_readnl
20 : !-----------------------------------------------------------------------
21 :
22 : use shr_kind_mod, only: r8 => shr_kind_r8, cl => shr_kind_cl
23 : use spmd_utils, only: masterproc,npes
24 : use pmgrid, only: plon, plat, plev, plevp
25 : use constituents, only: cnst_get_ind, pcnst, cnst_name
26 : use netcdf, only: NF90_NOERR,NF90_CLOSE,NF90_GET_VAR,NF90_INQUIRE_DIMENSION, &
27 : NF90_INQ_DIMID, NF90_INQ_VARID, NF90_NOWRITE, NF90_OPEN, &
28 : NF90_GET_ATT,NF90_GLOBAL,NF90_INQUIRE_ATTRIBUTE, &
29 : NF90_INQUIRE_VARIABLE, NF90_MAX_VAR_DIMS, nf90_get_var
30 : use shr_scam_mod, only: shr_scam_getCloseLatLon
31 : use cam_logfile, only: iulog
32 : use cam_abortutils, only: endrun
33 : use time_manager, only: get_curr_date, get_nstep,is_first_step,get_start_date,timemgr_time_inc
34 : use error_messages, only: handle_ncerr
35 :
36 :
37 : implicit none
38 : private
39 :
40 : ! PUBLIC INTERFACES:
41 :
42 : public :: scam_readnl ! read SCAM namelist options
43 : public :: readiopdata ! read iop boundary data
44 : public :: setiopupdate ! find index in iopboundary data for current time
45 : public :: plevs0 ! Define the pressures of the interfaces and midpoints
46 : public :: scmiop_flbc_inti
47 : public :: setiopupdate_init
48 :
49 : ! PUBLIC MODULE DATA:
50 :
51 : real(r8), public :: pressure_levels(plev)
52 : real(r8), public :: scmlat ! input namelist latitude for scam
53 : real(r8), public :: scmlon ! input namelist longitude for scam
54 : real(r8), public :: closeioplat ! closest iop latitude for scam
55 : real(r8), public :: closeioplon ! closest iop longitude for scam
56 : integer, public :: closeioplatidx ! file array index of closest iop latitude for scam
57 : integer, public :: closeioplonidx ! file array index closest iop longitude for scam
58 :
59 :
60 : integer, parameter :: num_switches = 20
61 : integer, parameter :: max_path_len = 128
62 :
63 : logical, public :: single_column ! Using IOP file or not
64 : logical, public :: use_iop ! Using IOP file or not
65 : logical, public :: use_pert_init ! perturb initial values
66 : logical, public :: use_pert_frc ! perturb forcing
67 : logical, public :: switch(num_switches) ! Logical flag settings from GUI
68 : logical, public :: l_uvphys ! If true, update u/v after TPHYS
69 : logical, public :: l_uvadvect ! If true, T, U & V will be passed to SLT
70 : logical, public :: l_conv ! use flux divergence terms for T and q?
71 : logical, public :: l_divtr ! use flux divergence terms for constituents?
72 : logical, public :: l_diag ! do we want available diagnostics?
73 :
74 : integer, public :: error_code ! Error code from netCDF reads
75 : integer, public :: initTimeIdx
76 : integer, public :: seedval
77 : integer :: bdate, last_date, last_sec
78 :
79 : character(len=max_path_len), public :: modelfile
80 : character(len=max_path_len), public :: analysisfile
81 : character(len=max_path_len), public :: sicfile
82 : character(len=max_path_len), public :: userfile
83 : character(len=max_path_len), public :: sstfile
84 : character(len=max_path_len), public :: lsmpftfile
85 : character(len=max_path_len), public :: pressfile
86 : character(len=max_path_len), public :: topofile
87 : character(len=max_path_len), public :: ozonefile
88 : character(len=max_path_len), public :: iopfile
89 : character(len=max_path_len), public :: absemsfile
90 : character(len=max_path_len), public :: aermassfile
91 : character(len=max_path_len), public :: aeropticsfile
92 : character(len=max_path_len), public :: timeinvfile
93 : character(len=max_path_len), public :: lsmsurffile
94 : character(len=max_path_len), public :: lsminifile
95 :
96 : ! note that scm_zadv_q is set to slt to be consistent with CAM BFB testing
97 :
98 :
99 : character(len=16), public :: scm_zadv_T = 'eulc '
100 : character(len=16), public :: scm_zadv_q = 'slt '
101 : character(len=16), public :: scm_zadv_uv = 'eulc '
102 :
103 : real(r8), public :: fixmascam
104 : real(r8), public :: betacam
105 : real(r8), public :: alphacam(pcnst)
106 : real(r8), public :: dqfxcam(plon,plev,pcnst)
107 :
108 : real(r8), public :: divq3d(plev,pcnst) ! 3D q advection
109 : real(r8), public :: divt3d(plev) ! 3D T advection
110 : real(r8), public :: divu3d(plev) ! 3D U advection
111 : real(r8), public :: divv3d(plev) ! 3D V advection
112 : real(r8), public :: vertdivq(plev,pcnst)! vertical q advection
113 : real(r8), public :: vertdivt(plev) ! vertical T advection
114 : real(r8), public :: vertdivu(plev) ! vertical T advection
115 : real(r8), public :: vertdivv(plev) ! vertical T advection
116 : real(r8), public :: ptend ! surface pressure tendency
117 : real(r8), public :: qdiff(plev) ! model minus observed humidity
118 : real(r8), public :: qobs(plev) ! actual W.V. Mixing ratio
119 : real(r8), public :: qinitobs(plev,pcnst)! initial tracer field
120 : real(r8), public :: cldliqobs(plev) ! actual W.V. Mixing ratio
121 : real(r8), public :: cldiceobs(plev) ! actual W.V. Mixing ratio
122 : real(r8), public :: numliqobs(plev) ! actual
123 : real(r8), public :: numiceobs(plev) ! actual
124 : real(r8), public :: precobs(1) ! observed precipitation
125 : real(r8), public :: lhflxobs(1) ! observed surface latent heat flux
126 : real(r8), public :: heat_glob_scm(1) ! observed heat total
127 : real(r8), public :: shflxobs(1) ! observed surface sensible heat flux
128 : real(r8), public :: q1obs(plev) ! observed apparent heat source
129 : real(r8), public :: q2obs(plev) ! observed apparent heat sink
130 : real(r8), public :: tdiff(plev) ! model minus observed temp
131 : real(r8), public :: tground(1) ! ground temperature
132 : real(r8), public :: psobs ! observed surface pressure
133 : real(r8), public :: tobs(plev) ! observed temperature
134 : real(r8), public :: tsair(1) ! air temperature at the surface
135 : real(r8), public :: udiff(plev) ! model minus observed uwind
136 : real(r8), public :: uobs(plev) ! actual u wind
137 : real(r8), public :: vdiff(plev) ! model minus observed vwind
138 : real(r8), public :: vobs(plev) ! actual v wind
139 : real(r8), public :: cldobs(plev) ! observed cld
140 : real(r8), public :: clwpobs(plev) ! observed clwp
141 : real(r8), public :: aldirobs(1) ! observed aldir
142 : real(r8), public :: aldifobs(1) ! observed aldif
143 : real(r8), public :: asdirobs(1) ! observed asdir
144 : real(r8), public :: asdifobs(1) ! observed asdif
145 :
146 : real(r8), public :: co2vmrobs(1) ! observed co2vmr
147 : real(r8), public :: ch4vmrobs(1) ! observed ch3vmr
148 : real(r8), public :: n2ovmrobs(1) ! observed n2ovmr
149 : real(r8), public :: f11vmrobs(1) ! observed f11vmr
150 : real(r8), public :: f12vmrobs(1) ! observed f12vmr
151 : real(r8), public :: soltsiobs(1) ! observed solar
152 :
153 : real(r8), public :: wfld(plev) ! Vertical motion (slt)
154 : real(r8), public :: wfldh(plevp) ! Vertical motion (slt)
155 : real(r8), public :: divq(plev,pcnst) ! Divergence of moisture
156 : real(r8), public :: divt(plev) ! Divergence of temperature
157 : real(r8), public :: divu(plev) ! Horiz Divergence of E/W
158 : real(r8), public :: divv(plev) ! Horiz Divergence of N/S
159 : ! mo_drydep algorithm
160 : real(r8), public, pointer :: loniop(:)
161 : real(r8), public, pointer :: latiop(:)
162 :
163 : integer, public :: iopTimeIdx ! index into iop dataset
164 : integer, public :: steplength ! Length of time-step
165 : integer, public :: base_date ! Date in (yyyymmdd) of start time
166 : integer, public :: base_secs ! Time of day of start time (sec)
167 :
168 : ! SCAM public data defaults
169 :
170 : logical, public :: doiopupdate = .false. ! do we need to read next iop timepoint
171 : logical, public :: have_lhflx = .false. ! dataset contains lhflx
172 : logical, public :: have_shflx = .false. ! dataset contains shflx
173 : logical, public :: have_heat_glob = .false. ! dataset contains heat total
174 : logical, public :: have_tg = .false. ! dataset contains tg
175 : logical, public :: have_tsair = .false. ! dataset contains tsair
176 : logical, public :: have_divq = .false. ! dataset contains divq
177 : logical, public :: have_divt = .false. ! dataset contains divt
178 : logical, public :: have_divq3d = .false. ! dataset contains divq3d
179 : logical, public :: have_vertdivu = .false. ! dataset contains vertdivu
180 : logical, public :: have_vertdivv = .false. ! dataset contains vertdivv
181 : logical, public :: have_vertdivt = .false. ! dataset contains vertdivt
182 : logical, public :: have_vertdivq = .false. ! dataset contains vertdivq
183 : logical, public :: have_divt3d = .false. ! dataset contains divt3d
184 : logical, public :: have_divu3d = .false. ! dataset contains divu3d
185 : logical, public :: have_divv3d = .false. ! dataset contains divv3d
186 : logical, public :: have_divu = .false. ! dataset contains divu
187 : logical, public :: have_divv = .false. ! dataset contains divv
188 : logical, public :: have_omega = .false. ! dataset contains omega
189 : logical, public :: have_phis = .false. ! dataset contains phis
190 : logical, public :: have_ptend = .false. ! dataset contains ptend
191 : logical, public :: have_ps = .false. ! dataset contains ps
192 : logical, public :: have_q = .false. ! dataset contains q
193 : logical, public :: have_q1 = .false. ! dataset contains Q1
194 : logical, public :: have_q2 = .false. ! dataset contains Q2
195 : logical, public :: have_prec = .false. ! dataset contains prec
196 : logical, public :: have_t = .false. ! dataset contains t
197 : logical, public :: have_u = .false. ! dataset contains u
198 : logical, public :: have_v = .false. ! dataset contains v
199 : logical, public :: have_cld = .false. ! dataset contains cld
200 : logical, public :: have_cldliq = .false. ! dataset contains cldliq
201 : logical, public :: have_cldice = .false. ! dataset contains cldice
202 : logical, public :: have_numliq = .false. ! dataset contains numliq
203 : logical, public :: have_numice = .false. ! dataset contains numice
204 : logical, public :: have_clwp = .false. ! dataset contains clwp
205 : logical, public :: have_aldir = .false. ! dataset contains aldir
206 : logical, public :: have_aldif = .false. ! dataset contains aldif
207 : logical, public :: have_asdir = .false. ! dataset contains asdir
208 : logical, public :: have_asdif = .false. ! dataset contains asdif
209 : logical, public :: use_camiop = .false. ! use cam generated forcing
210 : logical, public :: use_3dfrc = .false. ! use 3d forcing
211 : logical, public :: isrestart = .false. ! If this is a restart step or not
212 :
213 : ! SCAM namelist defaults
214 :
215 : logical, public :: scm_backfill_iop_w_init = .false. ! Backfill missing IOP data from initial file
216 : logical, public :: scm_relaxation = .false. ! Use relaxation
217 : logical, public :: scm_crm_mode = .false. ! Use column radiation mode
218 : logical, public :: scm_cambfb_mode = .false. ! Use extra CAM IOP fields to assure bit for bit match with CAM run
219 : logical, public :: scm_use_obs_T = .false. ! Use the SCAM-IOP observed T at each timestep instead of forecasting.
220 : logical, public :: scm_force_latlon = .false. ! force scam to use the lat lon fields specified in the namelist not closest
221 : real(r8), public :: scm_relaxation_low ! lowest level to apply relaxation
222 : real(r8), public :: scm_relaxation_high ! highest level to apply relaxation
223 : real(r8), public :: scm_relax_top_p = 0._r8 ! upper bound for scm relaxation
224 : real(r8), public :: scm_relax_bot_p = huge(1._r8) ! lower bound for scm relaxation
225 : real(r8), public :: scm_relax_tau_sec = 10800._r8 ! relaxation time constant (sec)
226 :
227 : ! +++BPM:
228 : ! modification... allow a linear ramp in relaxation time scale:
229 : logical, public :: scm_relax_linear = .false.
230 : real(r8), public :: scm_relax_tau_bot_sec = 10800._r8
231 : real(r8), public :: scm_relax_tau_top_sec = 10800._r8
232 : character(len=26), public :: scm_relax_fincl(pcnst)
233 :
234 : !
235 : ! note that scm_use_obs_uv is set to true to be consistent with CAM BFB testing
236 : !
237 :
238 : logical, public :: scm_use_obs_uv = .true. ! Use the SCAM-IOP observed u,v at each time step instead of forecasting.
239 :
240 : logical, public :: scm_use_obs_qv = .false. ! Use the SCAM-IOP observed qv at each time step instead of forecasting.
241 : logical, public :: scm_use_3dfrc = .false. ! Use CAMIOP 3d forcing if true, else use dycore vertical plus horizontal
242 : logical, public :: scm_iop_lhflxshflxTg = .false. !turn off LW rad
243 : logical, public :: scm_iop_Tg = .false. !turn off LW rad
244 :
245 : character(len=200), public :: scm_clubb_iop_name ! IOP name for CLUBB
246 :
247 : integer, allocatable, public :: tsec(:)
248 : integer, public :: ntime
249 :
250 : !=======================================================================
251 : contains
252 : !=======================================================================
253 :
254 1536 : subroutine scam_readnl(nlfile,single_column_in,scmlat_in,scmlon_in)
255 :
256 : use namelist_utils, only: find_group_name
257 : use units, only: getunit, freeunit
258 : use dycore, only: dycore_is
259 : use wrap_nf, only: wrap_open
260 :
261 :
262 : !---------------------------Arguments-----------------------------------
263 :
264 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input (nlfile=atm_in)
265 : logical, intent(in) :: single_column_in
266 : real(r8), intent(in) :: scmlat_in
267 : real(r8), intent(in) :: scmlon_in
268 :
269 : ! Local variables
270 : character(len=*), parameter :: sub = 'scam_readnl'
271 : integer :: unitn, ierr, i
272 : integer :: ncid
273 : integer :: iatt
274 : logical :: adv
275 :
276 : ! this list should include any variable that you might want to include in the namelist
277 : namelist /scam_nl/ iopfile, scm_iop_lhflxshflxTg, scm_iop_Tg, scm_relaxation, &
278 : scm_relax_top_p,scm_relax_bot_p,scm_relax_tau_sec, &
279 : scm_cambfb_mode,scm_crm_mode,scm_zadv_uv,scm_zadv_T,scm_zadv_q,&
280 : scm_use_obs_T, scm_use_obs_uv, scm_use_obs_qv, scm_use_3dfrc, &
281 : scm_relax_linear, scm_relax_tau_top_sec, &
282 : scm_relax_tau_bot_sec, scm_force_latlon, scm_relax_fincl, &
283 : scm_backfill_iop_w_init
284 :
285 1536 : single_column=single_column_in
286 :
287 1536 : iopfile = ' '
288 1536 : scm_clubb_iop_name = ' '
289 6144 : scm_relax_fincl(:) = ' '
290 1536 : if( single_column ) then
291 0 : if( npes>1) call endrun('SCAM_READNL: SCAM doesnt support using more than 1 pe.')
292 :
293 0 : if ( .not. (dycore_is('EUL') .or. dycore_is('SE')) .or. plon /= 1 .or. plat /=1 ) then
294 0 : call endrun('SCAM_SETOPTS: must compile model for SCAM mode when namelist parameter single_column is .true.')
295 : endif
296 :
297 0 : scmlat=scmlat_in
298 0 : scmlon=scmlon_in
299 :
300 0 : if( scmlat < -90._r8 .or. scmlat > 90._r8 ) then
301 0 : call endrun('SCAM_READNL: SCMLAT must be between -90. and 90. degrees.')
302 0 : elseif( scmlon < 0._r8 .or. scmlon > 360._r8 ) then
303 0 : call endrun('SCAM_READNL: SCMLON must be between 0. and 360. degrees.')
304 : end if
305 :
306 : ! Read namelist
307 0 : if (masterproc) then
308 0 : unitn = getunit()
309 0 : open( unitn, file=trim(nlfile), status='old' )
310 0 : call find_group_name(unitn, 'scam_nl', status=ierr)
311 0 : if (ierr == 0) then
312 0 : read(unitn, scam_nl, iostat=ierr)
313 0 : if (ierr /= 0) then
314 0 : call endrun(sub // ':: ERROR reading namelist')
315 : end if
316 : end if
317 0 : close(unitn)
318 0 : call freeunit(unitn)
319 : end if
320 :
321 : ! Error checking:
322 :
323 0 : iopfile = trim(iopfile)
324 0 : if( iopfile /= "" ) then
325 0 : use_iop = .true.
326 : else
327 0 : call endrun('SCAM_READNL: must specify IOP file for single column mode')
328 : endif
329 :
330 0 : call wrap_open( iopfile, NF90_NOWRITE, ncid )
331 :
332 0 : if( nf90_inquire_attribute( ncid, NF90_GLOBAL, 'CAM_GENERATED_FORCING', iatt ) == NF90_NOERR ) then
333 0 : use_camiop = .true.
334 : else
335 0 : use_camiop = .false.
336 : endif
337 :
338 : ! If we are not forcing the lat and lon from the namelist use the closest lat and lon that is found in the IOP file.
339 0 : if (.not.scm_force_latlon) then
340 0 : call shr_scam_GetCloseLatLon( ncid, scmlat, scmlon, closeioplat, closeioplon, closeioplatidx, closeioplonidx )
341 0 : write(iulog,*) 'SCAM_READNL: using closest IOP column to lat/lon specified in drv_in'
342 0 : write(iulog,*) ' requested lat,lon =',scmlat,', ',scmlon
343 0 : write(iulog,*) ' closest IOP lat,lon =',closeioplat,', ',closeioplon
344 0 : scmlat = closeioplat
345 0 : scmlon = closeioplon
346 : end if
347 :
348 0 : if (masterproc) then
349 0 : write (iulog,*) 'Single Column Model Options: '
350 0 : write (iulog,*) '============================='
351 0 : write (iulog,*) ' iopfile = ',trim(iopfile)
352 0 : write (iulog,*) ' scm_backfill_iop_w_init = ',scm_backfill_iop_w_init
353 0 : write (iulog,*) ' scm_cambfb_mode = ',scm_cambfb_mode
354 0 : write (iulog,*) ' scm_crm_mode = ',scm_crm_mode
355 0 : write (iulog,*) ' scm_force_latlon = ',scm_force_latlon
356 0 : write (iulog,*) ' scm_iop_Tg = ',scm_iop_Tg
357 0 : write (iulog,*) ' scm_iop_lhflxshflxTg = ',scm_iop_lhflxshflxTg
358 0 : write (iulog,*) ' scm_relaxation = ',scm_relaxation
359 0 : write (iulog,*) ' scm_relax_bot_p = ',scm_relax_bot_p
360 0 : write (iulog,*) ' scm_relax_linear = ',scm_relax_linear
361 0 : write (iulog,*) ' scm_relax_tau_bot_sec = ',scm_relax_tau_bot_sec
362 0 : write (iulog,*) ' scm_relax_tau_sec = ',scm_relax_tau_sec
363 0 : write (iulog,*) ' scm_relax_tau_top_sec = ',scm_relax_tau_top_sec
364 0 : write (iulog,*) ' scm_relax_top_p = ',scm_relax_top_p
365 0 : write (iulog,*) ' scm_use_obs_T = ',scm_use_obs_T
366 0 : write (iulog,*) ' scm_use_3dfrc = ',scm_use_3dfrc
367 0 : write (iulog,*) ' scm_use_obs_qv = ',scm_use_obs_qv
368 0 : write (iulog,*) ' scm_use_obs_uv = ',scm_use_obs_uv
369 0 : write (iulog,*) ' scm_zadv_T = ',trim(scm_zadv_T)
370 0 : write (iulog,*) ' scm_zadv_q = ',trim(scm_zadv_q)
371 0 : write (iulog,*) ' scm_zadv_uv = ',trim(scm_zadv_uv)
372 0 : write (iulog,*) ' scm_relax_finc: '
373 : ! output scm_relax_fincl character array
374 0 : do i=1,pcnst
375 0 : if (scm_relax_fincl(i) /= '') then
376 0 : adv = mod(i,4)==0
377 0 : if (adv) then
378 0 : write (iulog, "(A18)") "'"//trim(scm_relax_fincl(i))//"',"
379 : else
380 0 : write (iulog, "(A18)", ADVANCE="NO") "'"//trim(scm_relax_fincl(i))//"',"
381 : end if
382 : else
383 : exit
384 : end if
385 : end do
386 0 : print *
387 : end if
388 : end if
389 :
390 1536 : end subroutine scam_readnl
391 0 : subroutine readiopdata(hyam, hybm, hyai, hybi, ps0)
392 : !-----------------------------------------------------------------------
393 : !
394 : ! Open and read netCDF file containing initial IOP conditions
395 : !
396 : !---------------------------Code history--------------------------------
397 : !
398 : ! Written by J. Truesdale August, 1996, revised January, 1998
399 : !
400 : !-----------------------------------------------------------------------
401 : use getinterpnetcdfdata, only: getinterpncdata
402 : use string_utils, only: to_lower
403 : use wrap_nf, only: wrap_inq_dimid,wrap_get_vara_realx
404 : !-----------------------------------------------------------------------
405 : implicit none
406 :
407 : character(len=*), parameter :: sub = "read_iop_data"
408 : !
409 : !------------------------------Input Arguments--------------------------
410 : !
411 : real(r8),intent(in) :: hyam(plev),hybm(plev),hyai(plevp),hybi(plevp),ps0
412 : !
413 : !------------------------------Locals-----------------------------------
414 : !
415 : integer :: NCID, status
416 : integer :: time_dimID, lev_dimID, lev_varID, varid
417 : integer :: i,j
418 : integer :: nlev
419 : integer :: total_levs
420 : integer :: u_attlen
421 :
422 : integer :: k, m
423 : integer :: icldliq,icldice
424 : integer :: inumliq,inumice
425 :
426 : logical :: have_srf ! value at surface is available
427 : logical :: fill_ends !
428 : logical :: have_cnst(pcnst)
429 : real(r8) :: dummy
430 : real(r8) :: srf(1) ! value at surface
431 : real(r8) :: hyamiop(plev) ! a hybrid coef midpoint
432 : real(r8) :: hybmiop(plev) ! b hybrid coef midpoint
433 : real(r8) :: pmid(plev) ! pressure at model levels (time n)
434 : real(r8) :: pint(plevp) ! pressure at model interfaces (n )
435 : real(r8) :: pdel(plev) ! pdel(k) = pint (k+1)-pint (k)
436 : real(r8) :: weight
437 : real(r8) :: tmpdata(1)
438 : real(r8) :: coldata(plev)
439 0 : real(r8), allocatable :: dplevs( : )
440 : integer :: strt4(4),cnt4(4)
441 : integer :: nstep
442 : integer :: ios
443 : character(len=128) :: units ! Units
444 :
445 0 : nstep = get_nstep()
446 0 : fill_ends= .false.
447 :
448 : !
449 : ! Open IOP dataset
450 : !
451 : call handle_ncerr( nf90_open (iopfile, 0, ncid),&
452 0 : 'ERROR - scamMod.F90:readiopdata', __LINE__)
453 :
454 : !
455 : ! if the dataset is a CAM generated dataset set use_camiop to true
456 : ! CAM IOP datasets have a global attribute called CAM_GENERATED_IOP
457 : !
458 0 : if ( nf90_inquire_attribute( ncid, NF90_GLOBAL, 'CAM_GENERATED_FORCING', attnum=i )== NF90_NOERR ) then
459 0 : use_camiop = .true.
460 : else
461 0 : use_camiop = .false.
462 : endif
463 :
464 : !=====================================================================
465 : !
466 : ! Read time variables
467 :
468 :
469 0 : status = nf90_inq_dimid (ncid, 'time', time_dimID )
470 0 : if (status /= NF90_NOERR) then
471 0 : status = nf90_inq_dimid (ncid, 'tsec', time_dimID )
472 0 : if (status /= NF90_NOERR) then
473 0 : if (masterproc) write(iulog,*) sub//':ERROR - Could not find dimension ID for time/tsec'
474 0 : status = NF90_CLOSE ( ncid )
475 0 : call endrun(sub // ':ERROR - time/tsec must be present on the IOP file.')
476 : end if
477 : end if
478 :
479 : call handle_ncerr( nf90_inquire_dimension( ncid, time_dimID, len=ntime ),&
480 0 : 'Error - scamMod.F90:readiopdata unable to find time dimension', __LINE__)
481 :
482 : !
483 : !======================================================
484 : ! read level data
485 : !
486 0 : status = NF90_INQ_DIMID( ncid, 'lev', lev_dimID )
487 0 : if ( status /= nf90_noerr ) then
488 0 : if (masterproc) write(iulog,*) sub//':ERROR - Could not find variable dim ID for lev'
489 0 : status = NF90_CLOSE ( ncid )
490 0 : call endrun(sub // ':ERROR - Could not find variable dim ID for lev')
491 : end if
492 :
493 : call handle_ncerr( nf90_inquire_dimension( ncid, lev_dimID, len=nlev ),&
494 0 : 'Error - scamMod.f90:readiopdata unable to find level dimension', __LINE__)
495 :
496 0 : allocate(dplevs(nlev+1),stat=ios)
497 0 : if( ios /= 0 ) then
498 0 : write(iulog,*) sub//':ERROR: failed to allocate dplevs; error = ',ios
499 0 : call endrun(sub//':ERROR:readiopdata failed to allocate dplevs')
500 : end if
501 :
502 0 : status = NF90_INQ_VARID( ncid, 'lev', lev_varID )
503 0 : if ( status /= nf90_noerr ) then
504 0 : if (masterproc) write(iulog,*) sub//':ERROR - scamMod.F90:readiopdata:Could not find variable ID for lev'
505 0 : status = NF90_CLOSE ( ncid )
506 0 : call endrun(sub//':ERROR:ould not find variable ID for lev')
507 : end if
508 :
509 0 : call handle_ncerr( nf90_get_var (ncid, lev_varID, dplevs(:nlev)),&
510 0 : 'Error - scamMod.F90:readiopdata unable to read pressure levels', __LINE__)
511 : !
512 : !CAM generated forcing already has pressure on millibars convert standard IOP if needed.
513 : !
514 : call handle_ncerr(nf90_inquire_attribute(ncid, lev_varID, 'units', len=u_attlen),&
515 0 : 'Error - scamMod.F90:readiopdata unable to find units attribute', __LINE__)
516 : call handle_ncerr(nf90_get_att(ncid, lev_varID, 'units', units),&
517 0 : 'Error - scamMod.F90:readiopdata unable to read units attribute', __LINE__)
518 0 : units=trim(to_lower(units(1:u_attlen)))
519 :
520 0 : if ( units=='pa' .or. units=='pascal' .or. units=='pascals' ) then
521 : !
522 : ! convert pressure from Pascals to Millibars ( lev is expressed in pascals in iop datasets )
523 : !
524 0 : do i=1,nlev
525 0 : dplevs( i ) = dplevs( i )/100._r8
526 : end do
527 : endif
528 :
529 0 : status = nf90_inq_varid( ncid, 'Ps', varid )
530 0 : if ( status /= nf90_noerr ) then
531 0 : have_ps= .false.
532 0 : if (masterproc) write(iulog,*) sub//':Could not find variable Ps'
533 0 : if ( .not. scm_backfill_iop_w_init ) then
534 0 : status = NF90_CLOSE( ncid )
535 0 : call endrun(sub//':ERROR :IOP file must contain Surface Pressure (Ps) variable')
536 : else
537 0 : if ( is_first_step() .and. masterproc) write(iulog,*) 'Using surface pressure value from IC file if present'
538 : endif
539 : else
540 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
541 0 : status = nf90_get_var(ncid, varid, psobs, strt4)
542 0 : have_ps = .true.
543 : endif
544 :
545 :
546 : ! If the IOP dataset has hyam,hybm,etc it is assumed to be a hybrid level
547 : ! dataset
548 :
549 0 : status = nf90_inq_varid( ncid, 'hyam', varid )
550 0 : if ( status == nf90_noerr .and. have_ps) then
551 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
552 0 : status = nf90_get_var(ncid, varid, hyamiop, strt4)
553 0 : status = nf90_inq_varid( ncid, 'hybm', varid )
554 0 : status = nf90_get_var(ncid, varid, hybmiop, strt4)
555 0 : do i = 1, nlev
556 0 : dplevs( i ) = 1000.0_r8 * hyamiop( i ) + psobs * hybmiop( i ) / 100.0_r8
557 : end do
558 : endif
559 :
560 : ! add the surface pressure to the pressure level data, so that
561 : ! surface boundary condition will be set properly,
562 : ! making sure that it is the highest pressure in the array.
563 : !
564 :
565 0 : total_levs = nlev+1
566 0 : dplevs(nlev+1) = psobs/100.0_r8 ! ps is expressed in pascals
567 0 : do i= nlev, 1, -1
568 0 : if ( dplevs(i) > psobs/100.0_r8) then
569 0 : total_levs = i
570 0 : dplevs(i) = psobs/100.0_r8
571 : end if
572 : end do
573 0 : if (.not. use_camiop ) then
574 0 : nlev = total_levs
575 : endif
576 0 : if ( nlev == 1 ) then
577 0 : if (masterproc) write(iulog,*) sub//':Error - scamMod.F90:readiopdata: Ps too low!'
578 0 : call endrun(sub//':ERROR:Ps value on datasets is incongurent with levs data - mismatch in units?')
579 : endif
580 :
581 : !=====================================================================
582 : !get global vmrs from camiop file
583 0 : status = nf90_inq_varid( ncid, 'co2vmr', varid )
584 0 : if ( status == nf90_noerr) then
585 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
586 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,co2vmrobs)
587 : else
588 0 : if (is_first_step()) write(iulog,*)'using column value of co2vmr from boundary data as global volume mixing ratio'
589 : end if
590 0 : status = nf90_inq_varid( ncid, 'ch4vmr', varid )
591 0 : if ( status == nf90_noerr) then
592 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
593 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,ch4vmrobs)
594 : else
595 0 : if (is_first_step()) write(iulog,*)'using column value of ch4vmr from boundary data as global volume mixing ratio'
596 : end if
597 0 : status = nf90_inq_varid( ncid, 'n2ovmr', varid )
598 0 : if ( status == nf90_noerr) then
599 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
600 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,n2ovmrobs)
601 : else
602 0 : if (is_first_step()) write(iulog,*)'using column value of n2ovmr from boundary data as global volume mixing ratio'
603 : end if
604 0 : status = nf90_inq_varid( ncid, 'f11vmr', varid )
605 0 : if ( status == nf90_noerr) then
606 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
607 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,f11vmrobs)
608 : else
609 0 : if (is_first_step()) write(iulog,*)'using column value of f11vmr from boundary data as global volume mixing ratio'
610 : end if
611 0 : status = nf90_inq_varid( ncid, 'f12vmr', varid )
612 0 : if ( status == nf90_noerr) then
613 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
614 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,f12vmrobs)
615 : else
616 0 : if (is_first_step()) write(iulog,*)'using column value of f12vmr from boundary data as global volume mixing ratio'
617 : end if
618 0 : status = nf90_inq_varid( ncid, 'soltsi', varid )
619 0 : if ( status == nf90_noerr) then
620 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
621 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,soltsiobs)
622 : else
623 0 : if (is_first_step()) write(iulog,*)'using column value of soltsi from boundary data as global solar tsi'
624 : end if
625 : !=====================================================================
626 : !get state variables from camiop file
627 :
628 0 : status = nf90_inq_varid( ncid, 'Tsair', varid )
629 0 : if ( status /= nf90_noerr ) then
630 0 : have_tsair = .false.
631 : else
632 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
633 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,tsair)
634 0 : have_tsair = .true.
635 : endif
636 : !
637 : ! read in Tobs For cam generated iop readin small t to avoid confusion
638 : ! with capital T defined in cam
639 : !
640 0 : tobs(:)= 0._r8
641 :
642 0 : if ( use_camiop ) then
643 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx,'t', have_tsair, &
644 : tsair(1), fill_ends, scm_crm_mode, &
645 0 : dplevs, nlev,psobs, hyam, hybm,tobs, status )
646 : else
647 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx,'T', have_tsair, &
648 : tsair(1), fill_ends, scm_crm_mode, &
649 0 : dplevs, nlev,psobs, hyam, hybm, tobs, status )
650 : endif
651 0 : if ( status /= nf90_noerr ) then
652 0 : have_t = .false.
653 0 : if (masterproc) write(iulog,*) sub//':Could not find variable T on IOP file'
654 0 : if ( scm_backfill_iop_w_init ) then
655 0 : if (masterproc) write(iulog,*) sub//':Using value of T(tobs) from IC file if it exists'
656 : else
657 0 : if (masterproc) write(iulog,*) sub//':set tobs to 0.'
658 : endif
659 : !
660 : ! set T3 to Tobs on first time step
661 : !
662 : else
663 0 : have_t = .true.
664 : endif
665 :
666 0 : status = nf90_inq_varid( ncid, 'Tg', varid )
667 0 : if (status /= nf90_noerr) then
668 0 : if (masterproc) write(iulog,*) sub//':Could not find variable Tg on IOP dataset'
669 0 : if ( have_tsair ) then
670 0 : if (masterproc) write(iulog,*) sub//':Using Tsair'
671 0 : tground = tsair ! use surface value from T field
672 0 : have_Tg = .true.
673 : else
674 0 : have_Tg = .true.
675 0 : if (masterproc) write(iulog,*) sub//':Using T at lowest level from IOP dataset'
676 0 : tground = tobs(plev)
677 : endif
678 : else
679 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
680 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,tground)
681 0 : have_Tg = .true.
682 : endif
683 :
684 0 : status = nf90_inq_varid( ncid, 'qsrf', varid )
685 :
686 0 : if ( status /= nf90_noerr ) then
687 0 : have_srf = .false.
688 : else
689 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
690 0 : status = nf90_get_var(ncid, varid, srf(1), strt4)
691 0 : have_srf = .true.
692 : endif
693 :
694 0 : qobs(:)= 0._r8
695 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'q', have_srf, &
696 : srf(1), fill_ends, scm_crm_mode, &
697 0 : dplevs, nlev,psobs, hyam, hybm, qobs, status )
698 0 : if ( status /= nf90_noerr ) then
699 0 : have_q = .false.
700 0 : if (masterproc) write(iulog,*) sub//':Could not find variable q on IOP file'
701 0 : if ( scm_backfill_iop_w_init ) then
702 0 : if (masterproc) write(iulog,*) sub//':Using values for q from IC file if available'
703 : else
704 0 : if (masterproc) write(iulog,*) sub//':Setting qobs to 0.'
705 : endif
706 : else
707 0 : have_q = .true.
708 : endif
709 :
710 0 : cldobs = 0._r8
711 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'cld', .false., &
712 0 : dummy, fill_ends, scm_crm_mode, dplevs, nlev,psobs, hyam, hybm, cldobs, status )
713 0 : if ( status /= nf90_noerr ) then
714 0 : have_cld = .false.
715 : else
716 0 : have_cld = .true.
717 : endif
718 :
719 0 : clwpobs = 0._r8
720 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'clwp', .false., &
721 0 : dummy, fill_ends, scm_crm_mode, dplevs, nlev,psobs, hyam, hybm, clwpobs, status )
722 0 : if ( status /= nf90_noerr ) then
723 0 : have_clwp = .false.
724 : else
725 0 : have_clwp = .true.
726 : endif
727 :
728 : !
729 : ! read divq (horizontal advection)
730 : !
731 0 : status = nf90_inq_varid( ncid, 'divqsrf', varid )
732 0 : if ( status /= nf90_noerr ) then
733 0 : have_srf = .false.
734 : else
735 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
736 0 : status = nf90_get_var(ncid, varid, srf(1), strt4)
737 0 : have_srf = .true.
738 : endif
739 :
740 0 : divq(:,:)=0._r8
741 :
742 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, &
743 : 'divq', have_srf, srf(1), fill_ends, scm_crm_mode, &
744 0 : dplevs, nlev,psobs, hyam, hybm, divq(:,1), status )
745 0 : if ( status /= nf90_noerr ) then
746 0 : have_divq = .false.
747 : else
748 0 : have_divq = .true.
749 : endif
750 :
751 : !
752 : ! read vertdivq if available
753 : !
754 0 : status = nf90_inq_varid( ncid, 'vertdivqsrf', varid )
755 0 : if ( status /= nf90_noerr ) then
756 0 : have_srf = .false.
757 : else
758 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
759 0 : status = nf90_get_var(ncid, varid, srf(1), strt4)
760 0 : have_srf = .true.
761 : endif
762 :
763 0 : vertdivq=0._r8
764 :
765 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'vertdivq', &
766 : have_srf, srf(1), fill_ends, scm_crm_mode, &
767 0 : dplevs, nlev,psobs, hyam, hybm, vertdivq(:,1), status )
768 0 : if ( status /= nf90_noerr ) then
769 0 : have_vertdivq = .false.
770 : else
771 0 : have_vertdivq = .true.
772 : endif
773 :
774 0 : status = nf90_inq_varid( ncid, 'vertdivqsrf', varid )
775 0 : if ( status /= nf90_noerr ) then
776 0 : have_srf = .false.
777 : else
778 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
779 0 : status = nf90_get_var(ncid, varid, srf(1), strt4)
780 0 : have_srf = .true.
781 : endif
782 : !
783 : ! add calls to get dynamics tendencies for all prognostic consts
784 : !
785 0 : divq3d=0._r8
786 :
787 0 : do m = 1, pcnst
788 0 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, trim(cnst_name(m))//'_dten', &
789 : have_srf, srf(1), fill_ends, scm_crm_mode, &
790 0 : dplevs, nlev,psobs, hyam, hybm, divq3d(:,m), status )
791 0 : write(iulog,*)'checking ',trim(cnst_name(m))//'_dten',status
792 0 : if ( status /= nf90_noerr ) then
793 0 : have_cnst(m) = .false.
794 0 : divq3d(1:,m)=0._r8
795 : else
796 0 : if (m==1) have_divq3d = .true.
797 0 : have_cnst(m) = .true.
798 : endif
799 :
800 0 : coldata = 0._r8
801 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, trim(cnst_name(m))//'_dqfx', &
802 : have_srf, srf(1), fill_ends, scm_crm_mode, &
803 0 : dplevs, nlev,psobs, hyam, hybm, coldata, status )
804 0 : if ( STATUS /= NF90_NOERR ) then
805 0 : dqfxcam(1,:,m)=0._r8
806 : else
807 0 : dqfxcam(1,:,m)=coldata(:)
808 : endif
809 :
810 0 : tmpdata = 0._r8
811 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, trim(cnst_name(m))//'_alph', &
812 : have_srf, srf(1), fill_ends, scm_crm_mode, &
813 0 : dplevs, nlev,psobs, hyam, hybm, tmpdata, status )
814 0 : if ( status /= nf90_noerr ) then
815 0 : alphacam(m)=0._r8
816 : else
817 0 : alphacam(m)=tmpdata(1)
818 : endif
819 :
820 : end do
821 :
822 :
823 0 : numliqobs = 0._r8
824 0 : call cnst_get_ind('NUMLIQ', inumliq, abort=.false.)
825 0 : if ( inumliq > 0 ) then
826 0 : have_srf = .false.
827 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'NUMLIQ', &
828 : have_srf, srf(1), fill_ends, scm_crm_mode, &
829 0 : dplevs, nlev,psobs, hyam, hybm, numliqobs, status )
830 0 : if ( status /= nf90_noerr ) then
831 0 : have_numliq = .false.
832 : else
833 0 : have_numliq = .true.
834 : endif
835 : else
836 0 : have_numliq = .false.
837 : end if
838 :
839 0 : have_srf = .false.
840 :
841 0 : cldliqobs = 0._r8
842 0 : call cnst_get_ind('CLDLIQ', icldliq, abort=.false.)
843 0 : if ( icldliq > 0 ) then
844 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'CLDLIQ', &
845 : have_srf, srf(1), fill_ends, scm_crm_mode, &
846 0 : dplevs, nlev,psobs, hyam, hybm, cldliqobs, status )
847 0 : if ( status /= nf90_noerr ) then
848 0 : have_cldliq = .false.
849 : else
850 0 : have_cldliq = .true.
851 : endif
852 : else
853 0 : have_cldliq = .false.
854 : endif
855 :
856 0 : cldiceobs = 0._r8
857 0 : call cnst_get_ind('CLDICE', icldice, abort=.false.)
858 0 : if ( icldice > 0 ) then
859 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'CLDICE', &
860 : have_srf, srf(1), fill_ends, scm_crm_mode, &
861 0 : dplevs, nlev,psobs, hyam, hybm, cldiceobs, status )
862 0 : if ( status /= nf90_noerr ) then
863 0 : have_cldice = .false.
864 : else
865 0 : have_cldice = .true.
866 : endif
867 : else
868 0 : have_cldice = .false.
869 : endif
870 :
871 0 : numiceobs = 0._r8
872 0 : call cnst_get_ind('NUMICE', inumice, abort=.false.)
873 0 : if ( inumice > 0 ) then
874 0 : have_srf = .false.
875 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'NUMICE', &
876 : have_srf, srf(1), fill_ends, scm_crm_mode, &
877 0 : dplevs, nlev,psobs, hyam, hybm, numiceobs, status )
878 0 : if ( status /= nf90_noerr ) then
879 0 : have_numice = .false.
880 : else
881 0 : have_numice = .true.
882 : endif
883 : else
884 0 : have_numice = .false.
885 : end if
886 :
887 : !
888 : ! read divu (optional field)
889 : !
890 0 : status = nf90_inq_varid( ncid, 'divusrf', varid )
891 0 : if ( status /= nf90_noerr ) then
892 0 : have_srf = .false.
893 : else
894 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
895 0 : status = nf90_get_var(ncid, varid, srf(1), strt4)
896 0 : have_srf = .true.
897 : endif
898 :
899 0 : divu = 0._r8
900 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divu', &
901 : have_srf, srf(1), fill_ends, scm_crm_mode, &
902 0 : dplevs, nlev,psobs, hyam, hybm, divu, status )
903 0 : if ( status /= nf90_noerr ) then
904 0 : have_divu = .false.
905 : else
906 0 : have_divu = .true.
907 : endif
908 : !
909 : ! read divv (optional field)
910 : !
911 0 : status = nf90_inq_varid( ncid, 'divvsrf', varid )
912 0 : if ( status /= nf90_noerr ) then
913 0 : have_srf = .false.
914 : else
915 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
916 0 : status = nf90_get_var(ncid, varid, srf(1), strt4)
917 0 : have_srf = .true.
918 : endif
919 :
920 0 : divv = 0._r8
921 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divv', &
922 : have_srf, srf(1), fill_ends, scm_crm_mode, &
923 0 : dplevs, nlev,psobs, hyam, hybm, divv, status )
924 0 : if ( status /= nf90_noerr ) then
925 0 : have_divv = .false.
926 : else
927 0 : have_divv = .true.
928 : endif
929 : !
930 : ! read divt (optional field)
931 : !
932 0 : status = nf90_inq_varid( ncid, 'divtsrf', varid )
933 0 : if ( status /= nf90_noerr ) then
934 0 : have_srf = .false.
935 : else
936 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
937 0 : status = nf90_get_var(ncid, varid, srf(1), strt4)
938 0 : have_srf = .true.
939 : endif
940 :
941 0 : divt=0._r8
942 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, &
943 : 'divT', have_srf, srf(1), fill_ends, scm_crm_mode, &
944 0 : dplevs, nlev,psobs, hyam, hybm, divt, status )
945 0 : if ( status /= nf90_noerr ) then
946 0 : have_divt = .false.
947 : else
948 0 : have_divt = .true.
949 : endif
950 :
951 : !
952 : ! read vertdivt if available
953 : !
954 0 : status = nf90_inq_varid( ncid, 'vertdivTsrf', varid )
955 0 : if ( status /= nf90_noerr ) then
956 0 : have_srf = .false.
957 : else
958 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
959 0 : status = nf90_get_var(ncid, varid, srf(1), strt4)
960 0 : have_srf = .true.
961 : endif
962 :
963 0 : vertdivt=0._r8
964 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'vertdivTx', &
965 : have_srf, srf(1), fill_ends, scm_crm_mode, &
966 0 : dplevs, nlev,psobs, hyam, hybm, vertdivt, status )
967 0 : if ( status /= nf90_noerr ) then
968 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'vertdivT', &
969 : have_srf, srf(1), fill_ends, scm_crm_mode, &
970 0 : dplevs, nlev,psobs, hyam, hybm, vertdivt, status )
971 0 : if ( status /= nf90_noerr ) then
972 0 : have_vertdivt = .false.
973 : else
974 0 : have_vertdivt = .true.
975 : endif
976 : else
977 0 : have_vertdivt = .true.
978 : endif
979 : !
980 : ! read divt3d (combined vertical/horizontal advection)
981 : ! (optional field)
982 :
983 0 : status = nf90_inq_varid( ncid, 'divT3dsrf', varid )
984 0 : if ( status /= nf90_noerr ) then
985 0 : have_srf = .false.
986 : else
987 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
988 0 : status = nf90_get_var(ncid, varid, srf(1), strt4)
989 0 : have_srf = .true.
990 : endif
991 :
992 0 : divT3d = 0._r8
993 :
994 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divT3d', &
995 : have_srf, srf(1), fill_ends, scm_crm_mode, &
996 0 : dplevs, nlev,psobs, hyam, hybm, divt3d, status )
997 0 : write(iulog,*)'checking divT3d:',status,nf90_noerr
998 0 : if ( status /= nf90_noerr ) then
999 0 : have_divt3d = .false.
1000 : else
1001 0 : have_divt3d = .true.
1002 : endif
1003 :
1004 0 : divU3d = 0._r8
1005 :
1006 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divU3d', &
1007 : have_srf, srf(1), fill_ends, scm_crm_mode, &
1008 0 : dplevs, nlev,psobs, hyam, hybm, divu3d, status )
1009 0 : if ( status /= nf90_noerr ) then
1010 0 : have_divu3d = .false.
1011 : else
1012 0 : have_divu3d = .true.
1013 : endif
1014 :
1015 0 : divV3d = 0._r8
1016 :
1017 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divV3d', &
1018 : have_srf, srf(1), fill_ends, scm_crm_mode, &
1019 0 : dplevs, nlev,psobs, hyam, hybm, divv3d, status )
1020 0 : if ( status /= nf90_noerr ) then
1021 0 : have_divv3d = .false.
1022 : else
1023 0 : have_divv3d = .true.
1024 : endif
1025 :
1026 0 : status = nf90_inq_varid( ncid, 'Ptend', varid )
1027 0 : if ( status /= nf90_noerr ) then
1028 0 : have_ptend = .false.
1029 0 : if (masterproc) write(iulog,*) sub//':Could not find variable Ptend. Setting to zero'
1030 0 : ptend = 0.0_r8
1031 : else
1032 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
1033 0 : status = nf90_get_var(ncid, varid, srf(1), strt4)
1034 0 : have_ptend = .true.
1035 0 : ptend= srf(1)
1036 : endif
1037 :
1038 0 : wfld=0._r8
1039 :
1040 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, &
1041 : 'omega', .true., ptend, fill_ends, scm_crm_mode, &
1042 0 : dplevs, nlev,psobs, hyam, hybm, wfld, status )
1043 0 : if ( status /= nf90_noerr ) then
1044 0 : have_omega = .false.
1045 0 : if (masterproc) write(iulog,*) sub//':Could not find variable omega on IOP'
1046 0 : if ( scm_backfill_iop_w_init ) then
1047 0 : if (masterproc) write(iulog,*) sub//'Using omega from IC file'
1048 : else
1049 0 : if (masterproc) write(iulog,*) sub//'setting Omega to 0. throughout the column'
1050 : endif
1051 : else
1052 0 : have_omega = .true.
1053 : endif
1054 0 : call plevs0(plev, psobs, ps0, hyam, hybm, hyai, hybi, pint, pmid ,pdel)
1055 : !
1056 : ! Build interface vector for the specified omega profile
1057 : ! (weighted average in pressure of specified level values)
1058 : !
1059 0 : wfldh(:) = 0.0_r8
1060 :
1061 0 : do k=2,plev
1062 0 : weight = (pint(k) - pmid(k-1))/(pmid(k) - pmid(k-1))
1063 0 : wfldh(k) = (1.0_r8 - weight)*wfld(k-1) + weight*wfld(k)
1064 : end do
1065 :
1066 0 : status = nf90_inq_varid( ncid, 'usrf', varid )
1067 0 : if ( status /= nf90_noerr ) then
1068 0 : have_srf = .false.
1069 : else
1070 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
1071 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,srf)
1072 0 : have_srf = .true.
1073 : endif
1074 :
1075 0 : uobs=0._r8
1076 :
1077 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, &
1078 : 'u', have_srf, srf(1), fill_ends, scm_crm_mode, &
1079 0 : dplevs, nlev,psobs, hyam, hybm, uobs, status )
1080 0 : if ( status /= nf90_noerr ) then
1081 0 : have_u = .false.
1082 : else
1083 0 : have_u = .true.
1084 : endif
1085 :
1086 0 : status = nf90_inq_varid( ncid, 'vsrf', varid )
1087 0 : if ( status /= nf90_noerr ) then
1088 0 : have_srf = .false.
1089 : else
1090 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
1091 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,srf)
1092 0 : have_srf = .true.
1093 : endif
1094 :
1095 0 : vobs=0._r8
1096 :
1097 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, &
1098 : 'v', have_srf, srf(1), fill_ends, scm_crm_mode, &
1099 0 : dplevs, nlev,psobs, hyam, hybm, vobs, status )
1100 0 : if ( status /= nf90_noerr ) then
1101 0 : have_v = .false.
1102 : else
1103 0 : have_v = .true.
1104 : endif
1105 :
1106 0 : status = nf90_inq_varid( ncid, 'Prec', varid )
1107 0 : if ( status /= nf90_noerr ) then
1108 0 : have_prec = .false.
1109 : else
1110 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
1111 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,precobs)
1112 0 : have_prec = .true.
1113 : endif
1114 :
1115 0 : q1obs = 0._r8
1116 :
1117 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'Q1', &
1118 : .false., dummy, fill_ends, scm_crm_mode, & ! datasets don't contain Q1 at surface
1119 0 : dplevs, nlev,psobs, hyam, hybm, q1obs, status )
1120 0 : if ( status /= nf90_noerr ) then
1121 0 : have_q1 = .false.
1122 : else
1123 0 : have_q1 = .true.
1124 : endif
1125 :
1126 0 : q1obs = 0._r8
1127 :
1128 : call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'Q2', &
1129 : .false., dummy, fill_ends, scm_crm_mode, & ! datasets don't contain Q2 at surface
1130 0 : dplevs, nlev,psobs, hyam, hybm, q1obs, status )
1131 0 : if ( status /= nf90_noerr ) then
1132 0 : have_q2 = .false.
1133 : else
1134 0 : have_q2 = .true.
1135 : endif
1136 :
1137 : ! Test for BOTH 'lhflx' and 'lh' without overwriting 'have_lhflx'.
1138 : ! Analagous changes made for the surface heat flux
1139 :
1140 0 : status = nf90_inq_varid( ncid, 'lhflx', varid )
1141 0 : if ( status /= nf90_noerr ) then
1142 0 : status = nf90_inq_varid( ncid, 'lh', varid )
1143 0 : if ( status /= nf90_noerr ) then
1144 0 : have_lhflx = .false.
1145 : else
1146 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
1147 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,lhflxobs)
1148 0 : have_lhflx = .true.
1149 : endif
1150 : else
1151 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
1152 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,lhflxobs)
1153 0 : have_lhflx = .true.
1154 : endif
1155 :
1156 0 : status = nf90_inq_varid( ncid, 'shflx', varid )
1157 0 : if ( status /= nf90_noerr ) then
1158 0 : status = nf90_inq_varid( ncid, 'sh', varid )
1159 0 : if ( status /= nf90_noerr ) then
1160 0 : have_shflx = .false.
1161 : else
1162 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
1163 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,shflxobs)
1164 0 : have_shflx = .true.
1165 : endif
1166 : else
1167 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
1168 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,shflxobs)
1169 0 : have_shflx = .true.
1170 : endif
1171 :
1172 : ! If REPLAY is used, then need to read in the global
1173 : ! energy fixer
1174 0 : status = nf90_inq_varid( ncid, 'heat_glob', varid )
1175 0 : if (status /= nf90_noerr) then
1176 0 : have_heat_glob = .false.
1177 : else
1178 0 : call wrap_get_vara_realx (ncid,varid,strt4,cnt4,heat_glob_scm)
1179 0 : have_heat_glob = .true.
1180 : endif
1181 :
1182 : !
1183 : ! fill in 3d forcing variables if we have both horizontal
1184 : ! and vertical components, but not the 3d
1185 : !
1186 0 : if ( .not. have_cnst(1) .and. have_divq .and. have_vertdivq ) then
1187 0 : do k=1,plev
1188 0 : do m=1,pcnst
1189 0 : divq3d(k,m) = divq(k,m) + vertdivq(k,m)
1190 : enddo
1191 : enddo
1192 0 : have_divq3d = .true.
1193 : endif
1194 :
1195 0 : if ( .not. have_divt3d .and. have_divt .and. have_vertdivt ) then
1196 0 : if (masterproc) write(iulog,*) sub//'Don''t have divt3d - using divt and vertdivt'
1197 0 : do k=1,plev
1198 0 : divt3d(k) = divt(k) + vertdivt(k)
1199 : enddo
1200 0 : have_divt3d = .true.
1201 : endif
1202 : !
1203 : ! make sure that use_3dfrc flag is set to true if we only have
1204 : ! 3d forcing available
1205 : !
1206 0 : if (scm_use_3dfrc) then
1207 0 : if (have_divt3d .and. have_divq3d) then
1208 0 : use_3dfrc = .true.
1209 : else
1210 0 : call endrun(sub//':ERROR :IOP file must have both divt3d and divq3d forcing when scm_use_3dfrc is set to .true.')
1211 : endif
1212 : endif
1213 :
1214 0 : status = nf90_inq_varid( ncid, 'beta', varid )
1215 0 : if ( status /= nf90_noerr ) then
1216 0 : betacam = 0._r8
1217 : else
1218 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
1219 0 : status = nf90_get_var(ncid, varid, srf(1), strt4)
1220 0 : betacam=srf(1)
1221 : endif
1222 :
1223 0 : status = nf90_inq_varid( ncid, 'fixmas', varid )
1224 0 : if ( status /= nf90_noerr ) then
1225 0 : fixmascam=1.0_r8
1226 : else
1227 0 : call get_start_count(ncid, varid, scmlat, scmlon, ioptimeidx, strt4, cnt4)
1228 0 : status = nf90_get_var(ncid, varid, srf(1), strt4)
1229 0 : fixmascam=srf(1)
1230 : endif
1231 :
1232 0 : status = nf90_close( ncid )
1233 :
1234 0 : deallocate(dplevs)
1235 :
1236 0 : end subroutine readiopdata
1237 :
1238 0 : subroutine setiopupdate
1239 :
1240 : !-----------------------------------------------------------------------
1241 : !
1242 : ! Open and read netCDF file to extract time information
1243 : !
1244 : !---------------------------Code history--------------------------------
1245 : !
1246 : ! Written by John Truesdale August, 1996
1247 : !
1248 : !-----------------------------------------------------------------------
1249 : implicit none
1250 :
1251 : character(len=*), parameter :: sub = "setiopupdate"
1252 :
1253 : !------------------------------Locals-----------------------------------
1254 :
1255 : integer :: next_date, next_sec
1256 : integer :: ncsec,ncdate ! current time of day,date
1257 : integer :: yr, mon, day ! year, month, and day component
1258 : !------------------------------------------------------------------------------
1259 :
1260 0 : call get_curr_date(yr,mon,day,ncsec)
1261 0 : ncdate=yr*10000 + mon*100 + day
1262 :
1263 : !------------------------------------------------------------------------------
1264 : ! Check if iop data needs to be updated and set doiopupdate accordingly
1265 : !------------------------------------------------------------------------------
1266 :
1267 0 : if ( is_first_step() ) then
1268 0 : doiopupdate = .true.
1269 :
1270 : else
1271 :
1272 0 : call timemgr_time_inc(bdate, 0, next_date, next_sec, inc_s=tsec(iopTimeIdx+1))
1273 0 : if ( ncdate > next_date .or. (ncdate == next_date &
1274 : .and. ncsec >= next_sec)) then
1275 0 : doiopupdate = .true.
1276 : ! check to see if we need to move iopindex ahead more than 1 step
1277 0 : do while ( ncdate > next_date .or. (ncdate == next_date .and. ncsec >= next_sec))
1278 0 : iopTimeIdx = iopTimeIdx + 1
1279 0 : call timemgr_time_inc(bdate, 0, next_date, next_sec, inc_s=tsec(iopTimeIdx+1))
1280 : end do
1281 : #if DEBUG > 2
1282 : if (masterproc) write(iulog,*) sub//'nstep = ',get_nstep()
1283 : if (masterproc) write(iulog,*) sub//'ncdate=',ncdate,' ncsec=',ncsec
1284 : if (masterproc) write(iulog,*) sub//'next_date=',next_date,' next_sec=',next_sec
1285 : if (masterproc) write(iulog,*) sub//':******* do iop update'
1286 : #endif
1287 : else
1288 0 : doiopupdate = .false.
1289 : end if
1290 : endif ! if (endstep = 1 )
1291 : !
1292 : ! make sure we're
1293 : ! not going past end of iop data
1294 : !
1295 0 : if ( ncdate > last_date .or. (ncdate == last_date &
1296 : .and. ncsec > last_sec)) then
1297 0 : call endrun(sub//':ERROR: Reached the end of the time varient dataset')
1298 : endif
1299 :
1300 : #if DEBUG > 1
1301 : if (masterproc) write(iulog,*) sub//':iop time index = ' , ioptimeidx
1302 : #endif
1303 :
1304 0 : end subroutine setiopupdate
1305 :
1306 : !===============================================================================
1307 :
1308 0 : subroutine plevs0 (nver, ps, ps0, hyam, hybm, hyai, hybi, pint ,pmid ,pdel)
1309 :
1310 : !-----------------------------------------------------------------------
1311 : !
1312 : ! Purpose:
1313 : ! Define the pressures of the interfaces and midpoints from the
1314 : ! coordinate definitions and the surface pressure.
1315 : !
1316 : ! Author: B. Boville
1317 : !
1318 : !-----------------------------------------------------------------------
1319 : implicit none
1320 :
1321 :
1322 : !-----------------------------------------------------------------------
1323 : integer , intent(in) :: nver ! vertical dimension
1324 : real(r8), intent(in) :: ps ! Surface pressure (pascals)
1325 : real(r8), intent(in) :: ps0 ! reference pressure (pascals)
1326 : real(r8), intent(in) :: hyam(plev) ! hybrid midpoint coef
1327 : real(r8), intent(in) :: hybm(plev) ! hybrid midpoint coef
1328 : real(r8), intent(in) :: hyai(plevp) ! hybrid interface coef
1329 : real(r8), intent(in) :: hybi(plevp) ! hybrid interface coef
1330 : real(r8), intent(out) :: pint(nver+1) ! Pressure at model interfaces
1331 : real(r8), intent(out) :: pmid(nver) ! Pressure at model levels
1332 : real(r8), intent(out) :: pdel(nver) ! Layer thickness (pint(k+1) - pint(k))
1333 : !-----------------------------------------------------------------------
1334 :
1335 : !---------------------------Local workspace-----------------------------
1336 : integer :: k ! Longitude, level indices
1337 : !-----------------------------------------------------------------------
1338 : !
1339 : ! Set interface pressures
1340 : !
1341 : !$OMP PARALLEL DO PRIVATE (K)
1342 0 : do k=1,nver+1
1343 0 : pint(k) = hyai(k)*ps0 + hybi(k)*ps
1344 : end do
1345 : !
1346 : ! Set midpoint pressures and layer thicknesses
1347 : !
1348 : !$OMP PARALLEL DO PRIVATE (K)
1349 0 : do k=1,nver
1350 0 : pmid(k) = hyam(k)*ps0 + hybm(k)*ps
1351 0 : pdel(k) = pint(k+1) - pint(k)
1352 : end do
1353 :
1354 0 : end subroutine plevs0
1355 :
1356 0 : subroutine scmiop_flbc_inti ( co2vmr, ch4vmr, n2ovmr, f11vmr, f12vmr )
1357 : !-----------------------------------------------------------------------
1358 : !
1359 : ! Purpose:
1360 : ! Get start count for variable
1361 : !
1362 : !-----------------------------------------------------------------------
1363 :
1364 : implicit none
1365 :
1366 : real(r8), intent(out) :: co2vmr, ch4vmr, n2ovmr, f11vmr, f12vmr
1367 :
1368 : !-----------------------------------------------------------------------
1369 :
1370 0 : co2vmr=co2vmrobs(1)
1371 0 : ch4vmr=ch4vmrobs(1)
1372 0 : n2ovmr=n2ovmrobs(1)
1373 0 : f11vmr=f11vmrobs(1)
1374 0 : f12vmr=f12vmrobs(1)
1375 0 : end subroutine scmiop_flbc_inti
1376 :
1377 0 : subroutine get_start_count (ncid ,varid ,scmlat, scmlon, timeidx, start ,count)
1378 :
1379 : !-----------------------------------------------------------------------
1380 : !
1381 : ! Purpose:
1382 : ! set global lower boundary conditions
1383 : !
1384 : !-----------------------------------------------------------------------
1385 :
1386 : implicit none
1387 :
1388 : character(len=*), parameter :: sub = "get_start_count"
1389 :
1390 : !-----------------------------------------------------------------------
1391 : integer , intent(in) :: ncid ! file id
1392 : integer , intent(in) :: varid ! variable id
1393 : integer , intent(in) :: TimeIdx ! time index
1394 : real(r8), intent(in) :: scmlat,scmlon! scm lat/lon
1395 : integer , intent(out) :: start(:),count(:)
1396 :
1397 : !---------------------------Local workspace-----------------------------
1398 : integer :: dims_set,nlev,var_ndims
1399 : logical :: usable_var
1400 : character(len=cl) :: dim_name
1401 : integer :: var_dimIDs( NF90_MAX_VAR_DIMS )
1402 : real(r8) :: closelat,closelon
1403 : integer :: latidx,lonidx,status,i
1404 : !-----------------------------------------------------------------------
1405 :
1406 0 : call shr_scam_GetCloseLatLon(ncid,scmlat,scmlon,closelat,closelon,latidx,lonidx)
1407 :
1408 0 : STATUS = NF90_INQUIRE_VARIABLE( NCID, varID, ndims=var_ndims )
1409 : !
1410 : ! surface variables
1411 : !
1412 0 : if ( var_ndims == 0 ) then
1413 0 : call endrun(sub//':ERROR: var_ndims is 0 for varid:',varid)
1414 : endif
1415 :
1416 0 : STATUS = NF90_INQUIRE_VARIABLE( NCID, varID, dimids=var_dimIDs)
1417 0 : if ( STATUS /= NF90_NOERR ) then
1418 0 : write(iulog,* ) sub//'ERROR - Cant get dimension IDs for varid', varid
1419 0 : call endrun(sub//':ERROR: Cant get dimension IDs for varid',varid)
1420 : endif
1421 : !
1422 : ! Initialize the start and count arrays
1423 : !
1424 0 : dims_set = 0
1425 0 : nlev = 1
1426 0 : do i = var_ndims, 1, -1
1427 :
1428 0 : usable_var = .false.
1429 0 : STATUS = NF90_INQUIRE_DIMENSION( NCID, var_dimIDs( i ), dim_name )
1430 :
1431 0 : if ( trim(dim_name) == 'lat' ) then
1432 0 : start( i ) = latIdx
1433 0 : count( i ) = 1 ! Extract a single value
1434 0 : dims_set = dims_set + 1
1435 0 : usable_var = .true.
1436 : endif
1437 :
1438 0 : if ( trim(dim_name) == 'lon' .or. trim(dim_name) == 'ncol' .or. trim(dim_name) == 'ncol_d' ) then
1439 0 : start( i ) = lonIdx
1440 0 : count( i ) = 1 ! Extract a single value
1441 0 : dims_set = dims_set + 1
1442 0 : usable_var = .true.
1443 : endif
1444 :
1445 0 : if ( trim(dim_name) == 'lev' ) then
1446 0 : STATUS = NF90_INQUIRE_DIMENSION( NCID, var_dimIDs( i ), len=nlev )
1447 0 : start( i ) = 1
1448 0 : count( i ) = nlev ! Extract all levels
1449 0 : dims_set = dims_set + 1
1450 0 : usable_var = .true.
1451 : endif
1452 :
1453 0 : if ( trim(dim_name) == 'ilev' ) then
1454 0 : STATUS = NF90_INQUIRE_DIMENSION( NCID, var_dimIDs( i ), len=nlev )
1455 0 : start( i ) = 1
1456 0 : count( i ) = nlev ! Extract all levels
1457 0 : dims_set = dims_set + 1
1458 0 : usable_var = .true.
1459 : endif
1460 :
1461 0 : if ( trim(dim_name) == 'time' .OR. trim(dim_name) == 'tsec' ) then
1462 0 : start( i ) = TimeIdx
1463 0 : count( i ) = 1 ! Extract a single value
1464 0 : dims_set = dims_set + 1
1465 0 : usable_var = .true.
1466 : endif
1467 : end do
1468 0 : end subroutine get_start_count
1469 :
1470 : !=========================================================================
1471 0 : subroutine setiopupdate_init
1472 :
1473 : !-----------------------------------------------------------------------
1474 : !
1475 : ! Open and read netCDF file to extract time information
1476 : ! This subroutine should be called at the first SCM time step
1477 : !
1478 : !---------------------------Code history--------------------------------
1479 : !
1480 : ! Written by John Truesdale August, 1996
1481 : ! Modified for E3SM by Peter Bogenschutz 2017 - onward
1482 : !
1483 : !-----------------------------------------------------------------------
1484 : implicit none
1485 :
1486 : !------------------------------Locals-----------------------------------
1487 :
1488 : integer :: NCID,i
1489 : integer :: tsec_varID, time_dimID
1490 : integer :: bdate_varID
1491 : integer :: STATUS
1492 : integer :: next_date, next_sec
1493 : integer :: ncsec,ncdate ! current time of day,date
1494 : integer :: yr, mon, day ! year, month, and day component
1495 : integer :: start_ymd,start_tod
1496 :
1497 : character(len=*), parameter :: sub = "setiopupdate_init"
1498 : !!------------------------------------------------------------------------------
1499 :
1500 : ! Open and read pertinent information from the IOP file
1501 :
1502 : call handle_ncerr( nf90_open (iopfile, 0, ncid),&
1503 0 : 'ERROR - scamMod.F90:setiopupdate_init Failed to open iop file', __LINE__)
1504 :
1505 : ! Read time (tsec) variable
1506 :
1507 0 : STATUS = NF90_INQ_VARID( NCID, 'tsec', tsec_varID )
1508 0 : if ( STATUS /= NF90_NOERR ) then
1509 0 : write(iulog,*)sub//':ERROR: Cant get variable ID for tsec'
1510 0 : STATUS = NF90_CLOSE ( NCID )
1511 0 : call endrun(sub//':ERROR: Cant get variable ID for tsec')
1512 : end if
1513 :
1514 0 : STATUS = NF90_INQ_VARID( NCID, 'bdate', bdate_varID )
1515 0 : if ( STATUS /= NF90_NOERR ) then
1516 0 : STATUS = NF90_INQ_VARID( NCID, 'basedate', bdate_varID )
1517 0 : if ( STATUS /= NF90_NOERR ) then
1518 0 : write(iulog,*)'ERROR - setiopupdate:Cant get variable ID for base date'
1519 0 : STATUS = NF90_CLOSE ( NCID )
1520 0 : call endrun(sub//':ERROR: Cant get variable ID for base date')
1521 : endif
1522 : endif
1523 :
1524 0 : STATUS = NF90_INQ_DIMID( NCID, 'time', time_dimID )
1525 0 : if ( STATUS /= NF90_NOERR ) then
1526 0 : STATUS = NF90_INQ_DIMID( NCID, 'tsec', time_dimID )
1527 0 : if ( STATUS /= NF90_NOERR ) then
1528 0 : write(iulog,* )'ERROR - setiopupdate:Could not find variable dim ID for time'
1529 0 : STATUS = NF90_CLOSE ( NCID )
1530 0 : call endrun(sub//':ERROR:Could not find variable dim ID for time')
1531 : end if
1532 : end if
1533 :
1534 0 : if ( STATUS /= NF90_NOERR ) &
1535 0 : write(iulog,*)'ERROR - setiopupdate:Cant get variable dim ID for time'
1536 :
1537 0 : STATUS = NF90_INQUIRE_DIMENSION( NCID, time_dimID, len=ntime )
1538 0 : if ( STATUS /= NF90_NOERR )then
1539 0 : write(iulog,*)'ERROR - setiopupdate:Cant get time dimlen'
1540 : endif
1541 :
1542 0 : if (.not.allocated(tsec)) allocate(tsec(ntime))
1543 :
1544 0 : STATUS = NF90_GET_VAR( NCID, tsec_varID, tsec )
1545 0 : if ( STATUS /= NF90_NOERR )then
1546 0 : write(iulog,*)'ERROR - setiopupdate:Cant get variable tsec'
1547 : endif
1548 0 : STATUS = NF90_GET_VAR( NCID, bdate_varID, bdate )
1549 0 : if ( STATUS /= NF90_NOERR )then
1550 0 : write(iulog,*)'ERROR - setiopupdate:Cant get variable bdate'
1551 : endif
1552 :
1553 : ! Close the netCDF file
1554 0 : STATUS = NF90_CLOSE( NCID )
1555 :
1556 : ! determine the last date in the iop dataset
1557 :
1558 0 : call timemgr_time_inc(bdate, 0, last_date, last_sec, inc_s=tsec(ntime))
1559 :
1560 : ! set the iop dataset index
1561 0 : iopTimeIdx=0
1562 0 : do i=1,ntime ! set the first ioptimeidx
1563 0 : call timemgr_time_inc(bdate, 0, next_date, next_sec, inc_s=tsec(i))
1564 0 : call get_start_date(yr,mon,day,start_tod)
1565 0 : start_ymd = yr*10000 + mon*100 + day
1566 :
1567 0 : if ( start_ymd > next_date .or. (start_ymd == next_date &
1568 0 : .and. start_tod >= next_sec)) then
1569 0 : iopTimeIdx = i
1570 : endif
1571 : enddo
1572 :
1573 0 : call get_curr_date(yr,mon,day,ncsec)
1574 0 : ncdate=yr*10000 + mon*100 + day
1575 :
1576 0 : if (iopTimeIdx == 0.or.iopTimeIdx >= ntime) then
1577 0 : call timemgr_time_inc(bdate, 0, next_date, next_sec, inc_s=tsec(1))
1578 0 : write(iulog,*) 'Error::setiopupdate: Current model time does not fall within IOP period'
1579 0 : write(iulog,*) ' Current CAM Date is ',ncdate,' and ',ncsec,' seconds'
1580 0 : write(iulog,*) ' IOP start is ',next_date,' and ',next_sec,'seconds'
1581 0 : write(iulog,*) ' IOP end is ',last_date,' and ',last_sec,'seconds'
1582 0 : call endrun(sub//':ERROR: Current model time does not fall within IOP period')
1583 : endif
1584 :
1585 0 : doiopupdate = .true.
1586 :
1587 0 : end subroutine setiopupdate_init
1588 :
1589 : end module scamMod
|
