#include module mo_dyninp !--------------------------------------------------------------------- ! ... dynamics input data and code for ncep netcdf input !--------------------------------------------------------------------- use mo_grid, only : plon, plev, plevp, plat use mo_control, only : dyn_soilw implicit none save integer, parameter :: mxdynflds = 38 real, parameter :: qmin = 2.e-12 character(len=168) :: & lpath, & ! local path rpath, & ! remote pathname filename, & ! dynamics filename active_filename, & ! active dynamics filename last_filename ! prior dynamics filename logical :: & cosbdyn, & ! .true. => cos blocked, cray format rmdyn ! .true. => remove old dynamics files integer :: & lundyn, & ! fortran unit number attached to file nfld, & ! number of fields in data records nlon, & ! number of longitudes nlat, & ! number of latitudes nlev, & ! number of model layers dyndate1, & ! date previous to current date (in yymmdd format) dynsec1, & ! seconds relative to dyndate1 ndyn1, & ! index in data arrays for previous time sample dyndate2, & ! current date (in yymmdd format) dynsec2, & ! seconds relative to dyndate2 ndyn2, & ! index in data arrays for current time sample dynoffset, & ! offset (seconds) added to (date,sec) info mfilh ! time sample number of the one most recently read. ! (maintained internally... not read from integer header) integer :: & ncid, & ! netcdf file id vid(mxdynflds) = 0, & ! ids for dynamics fields datevid(2), & ! ids for date and datesec ntim, & ! number of time samples in current file curr_ts, & ! time index for most recently read data ndyn_grid ! dynamics xform grid index integer, allocatable, dimension(:) :: & date, & ! dates in input file datesec ! seconds relative to date integer :: lat_lims(2) real, allocatable :: & psin(:,:,:,:), & ! surface pressure (input, pa) ts_avg_in(:,:,:,:), & ! surface temperature fsds_avg_in(:,:,:,:), & ! surface fsds tsin(:,:,:,:), & ! surface temperature (input, k) tauxin(:,:,:,:), & ! surface zonal stress (input, n/m2) tauyin(:,:,:,:), & ! surface meridional stress (input, n/m2) hflxin(:,:,:,:), & ! surface heat flux (input, w/m2) qflxin(:,:,:,:), & ! surface moisture flux (input, kg/m2/s) tin(:,:,:,:,:), & ! temperature (input, k) qin(:,:,:,:,:), & ! specific humidity (input, kg/kg) uin(:,:,:,:,:), & ! u wind( m/s ) vin(:,:,:,:,:), & ! v wind( m/s ) snowin(:,:,:,:), & ! snow depth (m) fsdsin(:,:,:,:), & ! downward short-wave (input, W/m2) soilwin(:,:,:,:) ! soil moisture fraction logical :: & rdt, & ! read t rdoro, & ! read oro rdphis, & ! read phis calcedot, & ! f => read etadot rdomega, & ! read omega limqin, & ! true => limit input q s.t. it doesn''t exceed sat mr. arvdiff, & ! read cgs and kvh for archived vertical diffusion divdiff, & ! read variables used to diagnose vertical diffusion (taux,tauy,shflx,lhflx) arconvccm, & ! read zmu, zmd, zeu, conbeta, coneta for archived ccm convection arcldphys, & ! read fields for archived cloud physics artiedtke, & ! read pdmfuen, pdmfude, pdmfden, pdmfdde for archived tiedtke convection ditiedtke, & ! read qflx for diagnosed tiedtke convection addqrad, & ! t => add radiative heating tendency to the temperature field supplied to the convection routine initialization = .true., & force_synch = .false., & close_file = .true. logical :: & mandatory(mxdynflds) = .false. character(len=32) :: & nam(mxdynflds) ! names of dynamics fields: character(len=32) :: dyn_flnm_prefix = 'h' character(len=16) :: dyn_flnm_suffix = '.nc' character(len=16) :: dyn_flnm_date_frmt = 'NNNN' real :: & scfac(mxdynflds) = 1. ! scale factor to convert units in input files to the internal units. type date_type integer :: start(4) integer :: digits(4) end type type(date_type) :: date_template data nam / & 'U', & ! 1) zonal wind component 'V', & ! 2) meridional wind component 'PS', & ! 3) surface pressure 'ORO', & ! 4) orography 'PHIS', & ! 5) surface geopotential 'ETADOT', & ! 6) etadot (vertical wind component, d(eta)/dt ) 'OMEGA', & ! 7) omega (vertical wind component, d(pressure)/dt ) 'T', & ! 8) temperature 'Q', & ! 9) specific humidity 'TS', & ! 10) surface temperature (K) 'TAUX', & ! 11) surface zonal stress (N/m2) 'TAUY', & ! 12) surface meridional stress (N/m2) 'SHFLX', & ! 13) surface heat flux (W/m2) 'QFLX', & ! 14) surface moisture flux (kg/m2/s) 'KVH', & ! 15) vertical diffusion coefficient 'CGS', & ! 16) counter-gradient coefficient 'ZMU', & ! 17) zhang mu2 from conv_ccm, kg/m2/s 'ZMD', & ! 18) zhang md2 from conv_ccm, kg/m2/s 'ZEU', & ! 19) zhang eu2 from conv_ccm, 1/s 'CONETA', & ! 20) hack convective mass flux (kg/m^2/s) 'CONBETA', & ! 21) hack overshoot parameter (fraction) 'QRAD', & ! 22) radiative heating tendency (K/s) 'CMFDQR', & ! 23) dq/dt due to convective rainout (kg/kg/s) 'CLDFRC', & ! 24) total cloud fraction 'CONCLD', & ! 25) convective cloud fraction 'NRAIN', & ! 26) dq/dt due to stratiform precip (kg/kg/s) 'NEVAPR', & ! 27) rate of evaporation of strat precip (kg/kg/s) 'GBCW', & ! 28) grid box cloud water (kg/kg) 'PDMFUEN', & ! 29) tiedtke entrainment for updroaft 'PDMFUDE', & ! 30) tiedtke detrainment for updroaft 'PDMFDEN', & ! 31) tiedtke entrainment for downdroaft 'PDMFDDE', & ! 32) tiedtke detrainment for downdroaft 'SNOWH', & ! 33) snow height 'FSDS', & ! 34) srf downward direct shortwave radiation 'SOILW', & ! 35) soil moisture 'ICE', & ! 36) sea ice 'TS_avg', & ! 37) average surface temperature 'FSDS_avg' / ! 38) average FSDS private :: readdata, intp2d, lotim, nextdyn contains subroutine inidyn( xplon, xplat, xplev, xplevp, icdate, & icsec, offset, lun, remove, xrdt, & xrdoro, xrdphis, xcalcedot, xrdomega, xarvdiff, & xdivdiff, xarconvccm, xaddqrad, xarcldphys, xartiedtke, & xditiedtke, dynnm, dynsc, hyam, hybm, & hyai, hybi, phis, oro, plonl, platl, pplon ) !--------------------------------------------------------------------------- ! ... initialize module, dynamics arrays, and a few other variables. ! the dynamics fields will be linearly interpolated to the initial time (icdate,icsec). !--------------------------------------------------------------------------- use netcdf use mo_mpi use mo_constants, only : d2r, r2d, pi, phi, lam use mo_regrider, only : regrid_inti, regrid_lat_limits use mo_file_utils, only : open_netcdf_file use mo_mpi, only : masternode, base_lat, thisnode use mo_calendar, only : addsec2dat, diffdat use mo_control, only : dyn_flsp, xactive_drydep, xactive_emissions use mo_control, only : dyn_has_ts_avg, dyn_has_fsds_avg use mo_charutl, only : lastsl, incstr implicit none !--------------------------------------------------------------------------- ! ... dummy arguments !--------------------------------------------------------------------------- integer, intent(in) :: & xplon, xplat, xplev, xplevp integer, intent(in) :: & plonl, platl, pplon character(len=32), intent(in) :: & dynnm(mxdynflds) ! names of dynamics fields: integer, intent(in) :: & icdate, & ! date of initial conditions in yymmdd format icsec, & ! seconds relative to date for initial conditions offset, & ! add offset (seconds) to dates read from headers lun ! unit number for dynamics input logical, intent(in) :: & remove, & ! t => remove old dynamics files xrdt, & ! read t xrdoro, & ! read oro from the dynamics files xrdphis, & ! read phis xcalcedot, & ! f => read etadot xrdomega, & ! read omega xarvdiff, & ! read cgs and kvh for archived vertical diffusion xdivdiff, & ! read variables used to diagnose vertical diffusion (taux,tauy,shflx,lhflx) xarconvccm, & ! read zmu, zmd, zeu, conbeta, coneta for archived ccm convection xaddqrad, & ! t => read radiative heating tendency xarcldphys, & ! read fields for archived cloud physics xartiedtke, & ! read pdmfuen, pdmfude, pdmfden, pdmfdde for archived ! tiedtke convection xditiedtke ! read qflx for diagnosed tiedtke convection real, intent(in) :: & dynsc(mxdynflds) ! scale factor to convert units in input files to the internal units real, dimension(plev), intent(out) :: & hyam, & ! hybrid a coefficient at layer midpoints. hybm ! hybrid b coefficient at layer midpoints. real, dimension(plevp), intent(out) :: & hyai, & ! hybrid a coefficient at layer interfaces. hybi ! hybrid b coefficient at layer interfaces. real, dimension(plonl,platl,pplon), intent(out) :: & phis, & ! surface geopotential units: m^2/s^2 oro ! orography, values: 0 over ocean, 1 over land, 2 over sea ice. !--------------------------------------------------------------------------- ! ... local variables !--------------------------------------------------------------------------- integer :: & i, ii, j, k, & istat, & ! return status astat, & ! allocation status vid1, & slen, & wings ! number of "wing" latitudes integer :: did(4) integer :: bcast_ints(4) real :: & tmp real, allocatable :: & dyn_lats(:), dyn_lons(:), wrk(:) !--------------------------------------------------------------------------- ! ... set variables in module !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- ! ... if local path not specified, take filename from remote path ! and put it in subdirectory dyn !--------------------------------------------------------------------------- lpath = dyn_flsp%local_path rpath = dyn_flsp%remote_path filename = dyn_flsp%nl_filename lundyn = lun rmdyn = remove ndyn1 = 1 ndyn2 = 2 dynoffset = offset mfilh = 0 rdt = xrdt rdoro = xrdoro rdphis = xrdphis calcedot = xcalcedot rdomega = xrdomega arvdiff = xarvdiff divdiff = xdivdiff arconvccm = xarconvccm arcldphys = xarcldphys artiedtke = xartiedtke ditiedtke = xditiedtke addqrad = xaddqrad do i = 1,mxdynflds if( dynnm(i)(1:1) /= ' ' ) then nam(i) = dynnm(i) end if scfac(i) = dynsc(i) end do !--------------------------------------------------------------------------- ! ... allocate module arrays !--------------------------------------------------------------------------- allocate( psin(plonl,-3:platl+4,pplon,2), & ts_avg_in(plonl,1:platl,pplon,2), & fsds_avg_in(plonl,1:platl,pplon,2), & tsin(plonl,1:platl,pplon,2), & tauxin(plonl,1:platl,pplon,2), & tauyin(plonl,1:platl,pplon,2), & hflxin(plonl,1:platl,pplon,2), & qflxin(plonl,1:platl,pplon,2), & tin(plonl,plev,1:platl,pplon,2), & qin(plonl,plev,1:platl,pplon,2), & uin(plonl,plev,-3:platl+4,pplon,2), & vin(plonl,plev,-2:platl+3,pplon,2),stat=astat ) if( astat /= 0 ) then write(*,*) 'inidyn: failed to allocate psin .. vin; error = ',astat call endrun end if if( xactive_drydep ) then allocate( soilwin(plonl,1:platl,pplon,2),stat=astat ) if( astat /= 0 ) then write(*,*) 'inidyn: failed to allocate soilwin; error = ',astat call endrun end if soilwin(:,:,:,:) = 0. end if psin(:,:,:,:) = 0. tsin(:,:,:,:) = 0. tauxin(:,:,:,:) = 0. tauyin(:,:,:,:) = 0. hflxin(:,:,:,:) = 0. qflxin(:,:,:,:) = 0. tin(:,:,:,:,:) = 0. qin(:,:,:,:,:) = qmin uin(:,:,:,:,:) = 0. vin(:,:,:,:,:) = 0. !--------------------------------------------------------------------------- ! ... open netcdf file !--------------------------------------------------------------------------- if( masternode ) then write(*,*) thisnode,': inidyn: opening filename = ',trim(filename) write(*,*) thisnode,': inidyn: lpath = ',trim(lpath) write(*,*) thisnode,': inidyn: rpath = ',trim(rpath) ncid = open_netcdf_file( filename, lpath, rpath, masteronly=.true. ) write(*,*) thisnode,': inidyn: opened netcdf file ',trim(lpath) // trim(filename) active_filename = filename end if !--------------------------------------------------------------------------- ! ... Set mandatory flag ! u, v, ps, q are mandatory !--------------------------------------------------------------------------- mandatory(1:3) = .true. ! U, V, PS mandatory( 4) = rdoro ! ORO mandatory( 5) = rdphis ! PHIS mandatory( 7) = rdomega ! OMEGA mandatory( 8) = rdt ! T mandatory( 9) = .true. ! Q mandatory(11) = divdiff ! TAUX mandatory(12) = divdiff ! TAUY mandatory(13) = divdiff ! SHFLX mandatory(14) = divdiff ! QFLX mandatory(15) = arvdiff ! KVH mandatory(16) = arvdiff ! CGS mandatory(17) = arconvccm ! ZMU mandatory(18) = arconvccm ! ZMD mandatory(19) = arconvccm ! ZEU mandatory(20) = arconvccm ! CONETA mandatory(21) = arconvccm ! CONBETA mandatory(22) = addqrad ! QRAD mandatory(23) = arcldphys ! CMFDQR mandatory(24) = arcldphys ! CLDFRC mandatory(25) = arcldphys ! CONCLD mandatory(26) = arcldphys ! NRAIN mandatory(27) = arcldphys ! NEVAPR mandatory(28) = arcldphys ! GBCW mandatory(29) = artiedtke ! PDMFUEN mandatory(30) = artiedtke ! PDMFUDE mandatory(31) = artiedtke ! PDMFDEN mandatory(32) = artiedtke ! PDMFDDE mandatory(33) = .true. ! SNOWH if( xactive_drydep .or. xactive_emissions ) then mandatory(34) = .true. ! FSDS end if !--------------------------------------------------------------------------- ! ... get variable ids for fields required to be in the input data !--------------------------------------------------------------------------- if( masternode ) then do i = 1,mxdynflds if( mandatory(i) ) then write(*,*) 'inidyn: looking for mandatory variable ',nam(i) call handle_ncerr( nf_inq_varid( ncid, nam(i), vid(i) ), & 'inidyn: '// trim(nam(i)) // & ' not found in dynamics input file' ) else write(*,*) 'inidyn: looking for optional variable ',nam(i) istat = nf_inq_varid( ncid, nam(i), vid(i) ) if( istat == NF_NOERR ) then write(*,*) 'inidyn: Found optional variable '// trim(nam(i)) // & ' in dynamics input file' else write(*,*) 'inidyn: did not find optional variable '// trim(nam(i)) // & ' in dynamics input file' vid(i) = 0 end if end if end do endif #ifdef USE_MPI call mpi_bcast( vid, mxdynflds, mpi_integer, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': inidyn: failed to bcast vid; error = ',istat call endrun endif #endif dyn_soilw = xactive_drydep .and. vid(35) /= 0 dyn_has_ts_avg = vid(37) /= 0 dyn_has_fsds_avg = vid(38) /= 0 !--------------------------------------------------------------------------- ! ... snow !--------------------------------------------------------------------------- if( mandatory(33) ) then allocate( snowin(plonl,1:platl,pplon,2),stat=astat ) if( astat /= 0 ) then write(*,*) 'inidyn: failed to allocate snowin; error = ',astat call endrun end if snowin(:,:,:,:) = 0. end if !--------------------------------------------------------------------------- ! ... direct solar short wave !--------------------------------------------------------------------------- if( mandatory(34) ) then allocate( fsdsin(plonl,1:platl,pplon,2),stat=astat ) if( astat /= 0 ) then write(*,*) 'inidyn: failed to allocate fsdsin; error = ',astat call endrun end if fsdsin(:,:,:,:) = 0. end if !--------------------------------------------------------------------------- ! ... determine dimensions by looking at zonal velocity variable !--------------------------------------------------------------------------- if( masternode ) then call handle_ncerr( nf_inq_vardimid( ncid, vid(1), did ), & 'inidyn: getting dimids for zonal velocity' ) call handle_ncerr( nf_inq_dimlen( ncid, did(1), nlon ), & 'inidyn: getting nlon' ) call handle_ncerr( nf_inq_dimlen( ncid, did(2), nlat ), & 'inidyn: getting nlat' ) call handle_ncerr( nf_inq_dimlen( ncid, did(3), nlev ), & 'inidyn: getting nlev' ) call handle_ncerr( nf_inq_dimlen( ncid, did(4), ntim ), & 'inidyn: getting ntim' ) #ifdef USE_MPI bcast_ints(:) = (/ nlon, nlat, nlev, ntim /) #endif end if #ifdef USE_MPI call mpi_bcast( bcast_ints, 4, mpi_integer, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': inidyn: failed to bcast data dims; error = ',istat call endrun endif if( .not. masternode ) then nlon = bcast_ints(1) nlat = bcast_ints(2) nlev = bcast_ints(3) ntim = bcast_ints(4) endif #endif !--------------------------------------------------------------------------- ! ... check vertical level count against model !--------------------------------------------------------------------------- if( masternode ) then if( nlev < plev ) then write(*,*) 'inidyn: dynamical dataset vertical level count - ',nlev write(*,*) ' is < model level count - ',plev call endrun end if !--------------------------------------------------------------------------- ! ... get date info !--------------------------------------------------------------------------- call handle_ncerr( nf_inq_varid( ncid, 'date', datevid(1) ), & 'inidyn: date not found in dynamics input file' ) call handle_ncerr( nf_inq_varid( ncid, 'datesec', datevid(2) ), & 'inidyn: datesec not found in dynamics input file' ) end if allocate( date(ntim), stat=astat ) if( astat /= 0 ) then write(*,*) 'inidyn: failed to allocate date' call endrun end if allocate( datesec(ntim), stat=astat ) if( astat /= 0 ) then write(*,*) 'inidyn: failed to allocate datesec' call endrun end if if( masternode ) then call handle_ncerr( nf_get_var_int( ncid, datevid(1), date ), 'inidyn: getting date' ) call handle_ncerr( nf_get_var_int( ncid, datevid(2), datesec ), 'inidyn: getting datesec' ) write(*,*) thisnode,': inidyn: date' write(*,'(10i10)') date(:) write(*,*) thisnode,': inidyn: datesec' write(*,'(10i10)') datesec(:) !--------------------------------------------------------------------------- ! ... adjust dates by the specified offset !--------------------------------------------------------------------------- do i = 1,ntim call addsec2dat( dynoffset, date(i), datesec(i) ) end do end if #ifdef USE_MPI call mpi_bcast( date, ntim, mpi_integer, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': inidyn: failed to bcast date; error = ',istat call endrun endif call mpi_bcast( datesec, ntim, mpi_integer, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': inidyn: failed to bcast datesec; error = ',istat call endrun endif #endif !--------------------------------------------------------------------------- ! ... find latest date that is earlier than (icdate,icsec) !--------------------------------------------------------------------------- curr_ts = lotim( icdate, icsec ) if( curr_ts == 0 ) then write(*,*) thisnode,': inidyn: first time sample after ic: ', date(1), datesec(1) call endrun else write(*,*) thisnode,': inidyn: curr_ts = ',curr_ts end if !--------------------------------------------------------------------------- ! ... allocate dyn lats, lons and read !--------------------------------------------------------------------------- allocate( dyn_lats(nlat), stat=astat ) if( astat /= 0 ) then write(*,*) 'inidyn: failed to allocate dyn_lats vector' call endrun end if allocate( dyn_lons(nlon), stat=astat ) if( astat /= 0 ) then write(*,*) 'inidyn: failed to allocate dyn_lons vector' call endrun end if if( masternode ) then call handle_ncerr( nf_inq_varid( ncid, 'lat', vid1 ), 'inidyn: lat dim not in dyn file' ) call handle_ncerr( nf_get_var_double( ncid, vid1, dyn_lats ), 'inidyn: getting dyn lats' ) call handle_ncerr( nf_inq_varid( ncid, 'lon', vid1 ), 'inidyn: lon dim not in dyn file' ) call handle_ncerr( nf_get_var_double( ncid, vid1, dyn_lons ), 'inidyn: getting dyn lons' ) end if #ifdef USE_MPI call mpi_bcast( dyn_lons, nlon, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': inidyn: failed to bcast dyn_lons; error = ',istat call endrun endif call mpi_bcast( dyn_lats, nlat, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': inidyn: failed to bcast dyn_lats; error = ',istat call endrun endif #endif !--------------------------------------------------------------------------- ! ... set the transform from inputs lats to simulation lats !--------------------------------------------------------------------------- if( masternode ) then write(*,*) thisnode,': inidyn: source lats (deg)' write(*,'(10f8.2)') dyn_lats write(*,*) thisnode,': inidyn: target lats (deg)' write(*,'(10f8.2)') phi*r2d write(*,*) thisnode,': inidyn: source lons (deg)' write(*,'(10f8.2)') dyn_lons write(*,*) thisnode,': inidyn: target lons (deg)' write(*,'(10f8.2)') lam*r2d endif dyn_lats(:nlat) = d2r * dyn_lats(:nlat) dyn_lons(:nlon) = d2r * dyn_lons(:nlon) #ifdef USE_MPI wings = 4 #else wings = 0 #endif ndyn_grid = regrid_inti( nlat, plat, nlon, plon, dyn_lons, lam, dyn_lats, phi, wings, platl ) if( ndyn_grid < 0 ) then write(*,*) 'inidyn: regrid_inti error = ',ndyn_grid call endrun elseif( ndyn_grid == 0 ) then write(*,*) 'inidyn: no regridding needed, regrid_inti using grid number ',ndyn_grid else write(*,*) 'inidyn: regrid_inti using grid number ',ndyn_grid end if if( ndyn_grid > 0 ) then lat_lims = regrid_lat_limits( ndyn_grid ) else lat_lims = (/ base_lat + 1, base_lat + platl /) end if write(*,'(i3,'' : inidyn: lat_lims = '',2i6)') thisnode,lat_lims if( masternode ) then !--------------------------------------------------------------------------- ! ... get vertical hybid coordinate coefficients !--------------------------------------------------------------------------- allocate( wrk(nlev+1), stat=astat ) if( astat /= 0 ) then write(*,*) 'inidyn: failed to allocate wrk vector' call endrun end if !--------------------------------------------------------------------------- ! ... check for "pure" pressure hybrid component !--------------------------------------------------------------------------- istat = nf_inq_varid( ncid, 'hyam', vid1 ) if( istat == NF_NOERR ) then call handle_ncerr( nf_get_var_double( ncid, vid1, wrk ), & 'inidyn: getting hyam' ) hyam(:plev) = wrk(nlev-plev+1:nlev) call handle_ncerr( nf_inq_varid( ncid, 'hyai', vid1 ), & 'inidyn: hyai variable not in dyn file' ) call handle_ncerr( nf_get_var_double( ncid, vid1, wrk ), & 'inidyn: getting hyai' ) hyai(:plevp) = wrk(nlev-plev+1:nlev+1) else write(*,*) 'inidyn: hyam not found. Assuming sigma coordinate -- hyam,hyai=0.' hyam(:plev) = 0. hyai(:plevp) = 0. end if call handle_ncerr( nf_inq_varid( ncid, 'hybm', vid1 ), & 'inidyn: hybm variable not in dyn file' ) call handle_ncerr( nf_get_var_double( ncid, vid1, wrk ), & 'inidyn: getting hybm' ) hybm(:plev) = wrk(nlev-plev+1:nlev) call handle_ncerr( nf_inq_varid( ncid, 'hybi', vid1 ), & 'inidyn: hybi variable not in dyn file' ) call handle_ncerr( nf_get_var_double( ncid, vid1, wrk ), & 'inidyn: getting hybi' ) hybi(:plevp) = wrk(nlev-plev+1:nlev+1) deallocate( wrk ) end if #ifdef USE_MPI call mpi_bcast( hyam, plev, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': inidyn: failed to bcast hyam; error = ',istat call endrun endif call mpi_bcast( hybm, plev, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': inidyn: failed to bcast hybm; error = ',istat call endrun endif call mpi_bcast( hyai, plevp, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': inidyn: failed to bcast hyai; error = ',istat call endrun endif call mpi_bcast( hybi, plevp, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': inidyn: failed to bcast hybi; error = ',istat call endrun endif #endif if( masternode ) then !--------------------------------------------------------------------------- ! ... crack the dyn flnm date template !--------------------------------------------------------------------------- call date_cracker( dyn_flnm_date_frmt ) endif !--------------------------------------------------------------------------- ! ... read data into time level 1 !--------------------------------------------------------------------------- ndyn1 = 1 dyndate1 = date(curr_ts) dynsec1 = datesec(curr_ts) call readdata( curr_ts, ndyn1, phis, plonl, platl, pplon ) !--------------------------------------------------------------------------- ! ... if next time sample is not in current file then get next file !--------------------------------------------------------------------------- if( diffdat( date(ntim), datesec(ntim), icdate, icsec ) >= 0. ) then force_synch = .true. last_filename = filename close_file = .true. call nextdyn if( masternode ) then dyn_flsp%nl_filename = active_filename endif force_synch = .false. if( masternode ) then write(*,*) 'inidyn: back from nextdyn with filename,active_filename = ', & trim(filename),', ',trim(active_filename) endif else rmdyn = .false. close_file = .false. call nextdyn if( masternode ) then write(*,*) 'inidyn: back from nextdyn with filename, active_filename = ', & trim(filename), ', ', trim(active_filename) endif close_file = .true. rmdyn = remove end if curr_ts = curr_ts + 1 if( masternode ) then !--------------------------------------------------------------------------- ! ... check that current time sample is an upper bound to the ic time !--------------------------------------------------------------------------- if( diffdat( icdate, icsec, date(curr_ts), datesec(curr_ts)) < 0. ) then write(*,*) 'inidyn: ic time not bounded by specified dynamics ' // & 'input file. check namelist settings.' call endrun end if end if !--------------------------------------------------------------------------- ! ... read data into time level 2 !--------------------------------------------------------------------------- dyndate2 = date(curr_ts) dynsec2 = datesec(curr_ts) ndyn2 = 2 call readdata( curr_ts, ndyn2, phis, plonl, platl, & pplon, oro ) initialization = .false. deallocate( dyn_lats, dyn_lons ) end subroutine inidyn subroutine readdata( itim, idyn, phis, plonl, platl, & pplon, oro ) !--------------------------------------------------------------------------- ! ... read the field data for the specified time sample !--------------------------------------------------------------------------- use netcdf use mo_mpi use mo_regrider, only : regrid_3d, regrid_2d use mo_control, only : xactive_drydep, xactive_emissions implicit none !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: & itim, & ! time index for requested record from netcdf file idyn, & ! index in dynamics input array plonl, & platl, & pplon real, intent(out) :: & phis(plonl,platl,pplon) ! surface geopotential units: m^2/s^2 real, optional, intent(out) :: & oro(plonl,platl,pplon) ! orography, values: 0 over ocean, 1 over land, !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- integer :: & i, j, jl, ju, jll, jul, m, istat, lats integer :: & n2d, n3d integer, save, dimension(3) :: & count2, & start2 integer, save, dimension(4) :: & count, & start real, target :: inp_2d(nlon,nlat) ! global input 2d field real, target :: inp_ice(nlon,nlat) ! global input 2d field real, target :: inp_3d(nlon,nlat,plev) ! global input 3d field real, pointer :: tmp2d(:,:) ! temp 2d field real, pointer :: tmp2d_ice(:,:) ! for ice field real, pointer :: tmp3d(:,:,:) ! temp 3d field at midpoints or interfaces real :: wrk3d(plon,plev,-3:platl+4) ! work 3d field real :: wrk2d(plon,-3:platl+4) ! work 2d field real :: wrk(plon,plat) ! work 2d field real :: wrk2d_ice(plon,platl) ! for ice field logical, save :: has_ice_fld = .false. logical, save :: entered = .false. if( .not. entered ) then if( masternode ) then count2 = (/ nlon, nlat, 1 /) if( plev < nlev ) then count = (/ nlon, nlat, plev, 1 /) else count = (/ nlon, nlat, nlev, 1 /) end if end if entered = .true. end if if( masternode ) then start2 = (/ 1, 1, itim /) if( plev < nlev ) then start = (/ 1, 1, nlev-plev+1, itim /) else start = (/ 1, 1, 1, itim /) end if write(*,*) thisnode,': readdata: diags' write(*,'(i3,'' : readdata : start2 = '',3i6)') thisnode,start2 write(*,'(i3,'' : readdata : count2 = '',3i6)') thisnode,count2 write(*,'(i3,'' : readdata : start = '',4i6)') thisnode,start write(*,'(i3,'' : readdata : count = '',4i6)') thisnode,count end if n2d = nlon*nlat n3d = n2d*plev tmp2d => inp_2d(:,lat_lims(1):lat_lims(2)) tmp3d => inp_3d(:,lat_lims(1):lat_lims(2),:) !----------------------------------------------------------------------- ! ... zonal wind component u (m/s) !----------------------------------------------------------------------- m = 1 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start, count, inp_3d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_3d, n3d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast u; error = ',istat call endrun endif #endif jl = max( 1,base_lat - 3 ) ju = min( plat,base_lat + platl + 4 ) jll = jl - base_lat jul = ju - base_lat lats = ju - jl + 1 write(*,'(i3,'' : readdata: jl,ju,jll,jul = '',4i6)') thisnode,jl,ju,jll,jul call regrid_3d( tmp3d, wrk3d, ndyn_grid, jl, ju, .false., scfac(m) ) uin(:,:,jll:jul,:,idyn) = reshape( wrk3d(:,:,-3:-4+lats), (/plonl,plev,jul-jll+1,pplon/), order=(/1,4,2,3/) ) !----------------------------------------------------------------------- ! ... surface pressure ps (pa) !----------------------------------------------------------------------- m = 3 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast ps; error = ',istat call endrun endif #endif call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) psin(:,jll:jul,:,idyn) = reshape( wrk2d(:,-3:-4+lats), (/plonl,jul-jll+1,pplon/), order=(/1,3,2/) ) #ifdef DEBUG write(*,*) thisnode,': psin from netcdf' write(*,'(5(1p,e20.13))') tmp2d(1,:) write(*,*) thisnode,': psin after regridding' write(*,'(5(1p,e20.13))') wrk2d(:,-3) write(*,*) thisnode,': psin after reshaping' write(*,'(5(1p,e20.13))') psin(:,jll,:,idyn) #endif !----------------------------------------------------------------------- ! ... meridional wind component v (m/s) !----------------------------------------------------------------------- m = 2 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start, count, inp_3d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_3d, n3d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast v; error = ',istat call endrun endif #endif jl = max( 1,base_lat - 2 ) ju = min( plat,base_lat + platl + 3 ) jll = jl - base_lat jul = ju - base_lat lats = ju - jl + 1 #ifdef DEBUG write(*,*) thisnode,': readdata: jl,ju,jll,jul = ',jl,ju,jll,jul #endif call regrid_3d( tmp3d, wrk3d, ndyn_grid, jl, ju, .false., scfac(m) ) vin(:,:,jll:jul,:,idyn) = reshape( wrk3d(:,:,-3:-4+lats), (/plonl,plev,jul-jll+1,pplon/), order=(/1,4,2,3/) ) jl = max( 1,base_lat + 1 ) ju = min( plat,base_lat + platl ) lats = ju - jl + 1 #ifdef DEBUG write(*,'(i3,'' : readdata: jl,ju = '',2i6)') thisnode,jl,ju #endif if( rdoro .and. present( oro ) ) then !----------------------------------------------------------------------- ! ... orography !----------------------------------------------------------------------- m = 4 wrk2d(:,:) = 0. if( masternode ) then write(*,'(i3, '' : readdata: read oro time sample for date = '',i8,'':'',i5)') thisnode,date(itim), datesec(itim) call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast oro; error = ',istat call endrun endif #endif tmp2d_ice => inp_ice(:,lat_lims(1):lat_lims(2)) m = 36 has_ice_fld = vid(m) > 0 .and. nam(m) == 'ICE' if( .not. has_ice_fld ) then tmp2d_ice(:,:) = 0. where( nint(tmp2d(:,:)) == 2 ) tmp2d_ice(:,:) = 1. tmp2d(:,:) = 0. endwhere else if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_ice ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_ice, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast ice; error = ',istat call endrun endif #endif end if call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) wrk2d_ice(:,:) = 0. call regrid_2d( tmp2d_ice, wrk2d_ice, ndyn_grid, jl, ju, .true., scfac(m) ) where( nint(wrk2d_ice(:,:)) == 1 ) wrk2d(:,-3:platl-4) = 2. endwhere oro = reshape( wrk2d(:,-3:platl-4), (/plonl,platl,pplon/), order=(/1,3,2/) ) end if if( rdphis ) then !----------------------------------------------------------------------- ! ... surface geopotential, phis, m2/s2 !----------------------------------------------------------------------- m = 5 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast phis; error = ',istat call endrun endif #endif call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) phis = reshape( wrk2d(:,-3:platl-4), (/plonl,platl,pplon/), order=(/1,3,2/) ) end if if( rdt ) then !----------------------------------------------------------------------- ! ... temperature, t, k !----------------------------------------------------------------------- m = 8 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start, count, inp_3d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_3d, n3d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast t; error = ',istat call endrun endif #endif call regrid_3d( tmp3d, wrk3d, ndyn_grid, jl, ju, .true., scfac(m) ) tin(:,:,:,:,idyn) = reshape( wrk3d(:,:,-3:-4+lats), (/plonl,plev,platl,pplon/), order=(/1,4,2,3/) ) end if !----------------------------------------------------------------------- ! ... specific humidity, q (kg/kg) !----------------------------------------------------------------------- m = 9 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start, count, inp_3d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_3d, n3d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast q; error = ',istat call endrun endif #endif call regrid_3d( tmp3d, wrk3d, ndyn_grid, jl, ju, .true., scfac(m) ) qin(:,:,:,:,idyn) = reshape( wrk3d(:,:,-3:-4+lats), (/plonl,plev,platl,pplon/), order=(/1,4,2,3/) ) !----------------------------------------------------------------------- ! ... limit q to be >= qmin. some of the convection/cloud codes blow ! up when given q=0. !----------------------------------------------------------------------- where( qin(:,:,:,:,idyn) < qmin ) qin(:,:,:,:,idyn) = qmin endwhere if( divdiff ) then !----------------------------------------------------------------------- ! ... surface temperature, ts, k !----------------------------------------------------------------------- if( nam(10) /= nam(8) ) then m = 10 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast ts; error = ',istat call endrun endif #endif call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) tsin(:,:,:,idyn) = reshape( wrk2d(:,-3:platl-4), (/plonl,platl,pplon/), order=(/1,3,2/) ) else tsin(:,:,:,idyn) = tin(:,plev,:,:,idyn) end if !----------------------------------------------------------------------- ! ... surface zonal stress, taux, n/m2 !----------------------------------------------------------------------- m = 11 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast taux; error = ',istat call endrun endif #endif call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) tauxin(:,:,:,idyn) = reshape( wrk2d(:,-3:platl-4), (/plonl,platl,pplon/), order=(/1,3,2/) ) !----------------------------------------------------------------------- ! ... surface meridional stress, tauy, n/m2 !----------------------------------------------------------------------- m = 12 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast tauy; error = ',istat call endrun endif #endif call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) tauyin(:,:,:,idyn) = reshape( wrk2d(:,-3:platl-4), (/plonl,platl,pplon/), order=(/1,3,2/) ) !----------------------------------------------------------------------- ! ... surface heat flux, shflx, w/m2 !----------------------------------------------------------------------- m = 13 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast shflx; error = ',istat call endrun endif #endif call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) hflxin(:,:,:,idyn) = reshape( wrk2d(:,-3:platl-4), (/plonl,platl,pplon/), order=(/1,3,2/) ) !----------------------------------------------------------------------- ! ... surface moisture flux, qflx, kg/m2/s !----------------------------------------------------------------------- m = 14 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast qflx; error = ',istat call endrun endif #endif call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) qflxin(:,:,:,idyn) = reshape( wrk2d(:,-3:platl-4), (/plonl,platl,pplon/), order=(/1,3,2/) ) !----------------------------------------------------------------------- ! ... read snow and downward shortwave radiation flux !----------------------------------------------------------------------- if( mandatory(33) ) then m = 33 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast snow; error = ',istat call endrun endif #endif call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) snowin(:,:,:,idyn) = reshape( wrk2d(:,-3:platl-4), (/plonl,platl,pplon/), order=(/1,3,2/) ) end if if( mandatory(34) ) then m = 34 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast fsds; error = ',istat call endrun endif #endif call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) fsdsin(:,:,:,idyn) = reshape( wrk2d(:,-3:platl-4), (/plonl,platl,pplon/), order=(/1,3,2/) ) end if if( xactive_drydep ) then if( dyn_soilw ) then m = 35 if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast soilw; error = ',istat call endrun endif #endif call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) soilwin(:,:,:,idyn) = reshape( wrk2d(:,-3:platl-4), (/plonl,platl,pplon/), order=(/1,3,2/) ) end if end if end if !----------------------------------------------------------------------- ! ... surface average temperature !----------------------------------------------------------------------- m = 37 if( vid(m) /= 0 ) then if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast ts_avg; error = ',istat call endrun endif #endif call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) ts_avg_in(:,:,:,idyn) = reshape( wrk2d(:,-3:platl-4), (/plonl,platl,pplon/), order=(/1,3,2/) ) end if !----------------------------------------------------------------------- ! ... surface average fsds !----------------------------------------------------------------------- m = 38 if( vid(m) /= 0 ) then if( masternode ) then call handle_ncerr( nf_get_vara_double( ncid, vid(m), start2, count2, inp_2d ), & 'readdata: getting '// trim( nam(m) ) ) endif #ifdef USE_MPI call mpi_bcast( inp_2d, n2d, mpi_double_precision, 0, mpi_comm_world, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': readdata: failed to bcast fsds_avg; error = ',istat call endrun endif #endif call regrid_2d( tmp2d, wrk2d, ndyn_grid, jl, ju, .true., scfac(m) ) fsds_avg_in(:,:,:,idyn) = reshape( wrk2d(:,-3:platl-4), (/plonl,platl,pplon/), order=(/1,3,2/) ) end if write(*,'(i3,'' : readdata: read time sample for date = '',i8,'':'',i5)') & thisnode, date(itim), datesec(itim) end subroutine readdata subroutine chkdyn( ncdate, ncsec, phis, oro, plonl, & platl, pplon ) !----------------------------------------------------------------------- ! ... check that date to interpolate the data to is contained ! in the interval of dynamics data currently in memory. ! if not, then read the next time sample from the dynamics file, ! acquiring the next file if necessary. !----------------------------------------------------------------------- use mo_calendar, only : diffdat use mo_mpi, only : masternode, thisnode use mo_control, only : dyn_flsp, dicldphys use mo_grid, only : nodes use mo_constants,only : phi use mo_cloud, only : inimland implicit none !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: & ncdate, & ! interpolate the input data to (ncdate,ncsec) ncsec, & plonl, & platl, & pplon real, dimension(plonl,platl,pplon), intent(inout) :: & phis, & ! surface geopotential oro ! orography !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- integer :: & istat, lo, hi if( diffdat( ncdate, ncsec, dyndate2, dynsec2 ) >= 0. ) then return end if !----------------------------------------------------------------------- ! (ncdate,ncsec) is outside of current dynamics time interval. ! get the time sample(s) for the new time interval. ! ! make the upper bound of the current time interval the lower bound of the ! new time interval. in the "usual" case we won''t need to read new ! data into this slot. !----------------------------------------------------------------------- ndyn1 = ndyn2 !----------------------------------------------------------------------- ! ... is the new upper bound time sample in the current dynamics file? !----------------------------------------------------------------------- lo = lotim( ncdate, ncsec ) hi = lo + 1 write(*,*) thisnode,': chkdyn: lo,hi,curr_ts = ',lo,hi,curr_ts if( hi <= ntim ) then ! all required data is present !----------------------------------------------------------------------- ! is the current upper bound time sample the same as the lower bound ! time sample for the new time interval? if not we need to acquire ! new lower bound data. !----------------------------------------------------------------------- if( curr_ts == lo ) then dyndate1 = dyndate2 dynsec1 = dynsec2 else dyndate1 = date(lo) dynsec1 = datesec(lo) call readdata( lo, ndyn1, phis, plonl, platl, pplon ) end if else !----------------------------------------------------------------------- ! ... need next file. if last time sample in current file is not the ! current upper bound then read it into the lower bound position ! before getting next file. !----------------------------------------------------------------------- if( curr_ts == ntim ) then dyndate1 = dyndate2 dynsec1 = dynsec2 else dyndate1 = date(ntim) dynsec1 = datesec(ntim) call readdata( ntim, ndyn1, phis, plonl, platl, pplon ) end if if( masternode ) then last_filename = active_filename active_filename = filename dyn_flsp%nl_filename = filename write(*,*) thisnode,': chkdyn: calling nextdyn with filename=',trim(filename) end if call nextdyn if( masternode ) then write(*,*) thisnode,': chkdyn: back from nextdyn with filename=',trim(filename) end if lo = lotim( ncdate, ncsec ) hi = lo + 1 !----------------------------------------------------------------------- ! ... make sure the new upper bound time sample is in this file. !----------------------------------------------------------------------- if( hi > ntim ) then write(*,*) thisnode,': chkdyn: current time not bounded by data on next' // & ' dynamics input file.' call endrun end if !----------------------------------------------------------------------- ! is the current lower bound time sample (i.e., the last time sample ! from the previous file) the correct one? if not we need to acquire ! new lower bound data. !----------------------------------------------------------------------- if( lo > 0 ) then dyndate1 = date(lo) dynsec1 = datesec(lo) call readdata( lo, ndyn1, phis, plonl, platl, pplon ) end if end if !----------------------------------------------------------------------- ! ... read data for new upper bound !----------------------------------------------------------------------- ndyn2 = mod( ndyn1,2 ) + 1 dyndate2 = date(hi) dynsec2 = datesec(hi) call readdata( hi, ndyn2, phis, plonl, platl, & pplon, oro ) curr_ts = hi !---------------------------------------------------------------------- ! ... set land mask !---------------------------------------------------------------------- if( dicldphys ) then write(*,*) 'chkdyn: calling inimland at time - ',ncdate,':',ncsec call inimland( oro, phi, plonl, platl, pplon, nodes ) end if end subroutine chkdyn subroutine interp1( ncdate, ncsec, ps, t, q, ts, & ts_avg, fsds_avg, plonl, platl, pplon ) !--------------------------------------------------------------------------- ! ... interpolate the fields whose values are required at the ! begining of a timestep. !--------------------------------------------------------------------------- use mo_mpi, only : masternode, lastnode use mo_calendar, only : diffdat implicit none !--------------------------------------------------------------------------- ! ... dummy arguments !--------------------------------------------------------------------------- integer, intent(in) :: & ncdate, & ! interpolate the input data to (ncdate,ncsec) ncsec, & plonl, & platl, & pplon real, intent(out) :: ps(plonl,-3:platl+4,pplon) ! surface pressure (time interpolant) real, intent(out) :: t(plonl,plev,platl,pplon) ! temperature (time interpolant) real, intent(out) :: q(plonl,plev,platl,pplon) ! specific humidity (time interpolant) real, intent(out) :: ts(plonl,platl,pplon) ! surface temperature real, intent(out) :: ts_avg(plonl,platl,pplon) ! surface temperature real, intent(out) :: fsds_avg(plonl,platl,pplon) ! surface fsds !--------------------------------------------------------------------------- ! ... local variables !--------------------------------------------------------------------------- integer :: & j, & ! index ip, & ! index jl, & ! index ju ! index real :: & dynday2, & ! (dyndate2,dynsec2) - (dyndate1,dynsec1) in days intpday ! (ncdate,ncsec) - (dyndate1,dynsec1) in days dynday2 = diffdat( dyndate1, dynsec1, dyndate2, dynsec2 ) intpday = diffdat( dyndate1, dynsec1, ncdate, ncsec ) #ifdef DEBUG write(*,*) 'interp1: timing info' write(*,*) 'dynday2 = ',dynday2 write(*,*) 'intpday = ',intpday #endif #ifdef USE_MPI jl = -3 ju = platl + 4 if( masternode ) then jl = 1 end if if( lastnode ) then ju = platl end if #else jl = 1 ju = platl #endif if( pplon > 1 ) then !$omp parallel do private( ip, j ) do ip = 1,pplon do j = jl,ju call intp2d( 1, 0., dynday2, intpday, & psin(1,j,ip,ndyn1), psin(1,j,ip,ndyn2), ps(1,j,ip), plonl ) end do do j = 1,platl call intp2d( plev, 0., dynday2, intpday, & qin(1,1,j,ip,ndyn1), qin(1,1,j,ip,ndyn2), q(1,1,j,ip), plonl ) if( rdt ) then call intp2d( plev, 0., dynday2, intpday, & tin(1,1,j,ip,ndyn1), tin(1,1,j,ip,ndyn2), t(1,1,j,ip), plonl ) end if if( divdiff ) then call intp2d( 1, 0., dynday2, intpday, & tsin(1,j,ip,ndyn1), tsin(1,j,ip,ndyn2), ts(1,j,ip), plonl ) end if if( vid(37) > 0 ) then call intp2d( 1, 0., dynday2, intpday, & ts_avg_in(1,j,ip,ndyn1), ts_avg_in(1,j,ip,ndyn2), ts_avg(1,j,ip), plonl ) end if if( vid(38) > 0 ) then call intp2d( 1, 0., dynday2, intpday, & fsds_avg_in(1,j,ip,ndyn1), fsds_avg_in(1,j,ip,ndyn2), fsds_avg(1,j,ip), plonl ) end if end do end do !$omp end parallel do else !$omp parallel private( j ) !$omp do do j = jl,ju call intp2d( 1, 0., dynday2, intpday, & psin(1,j,1,ndyn1), psin(1,j,1,ndyn2), ps(1,j,1), plonl ) end do !$omp end do !$omp do do j = 1,platl call intp2d( plev, 0., dynday2, intpday, & qin(1,1,j,1,ndyn1), qin(1,1,j,1,ndyn2), q(1,1,j,1), plonl ) if( rdt ) then call intp2d( plev, 0., dynday2, intpday, & tin(1,1,j,1,ndyn1), tin(1,1,j,1,ndyn2), t(1,1,j,1), plonl ) end if if( divdiff ) then call intp2d( 1, 0., dynday2, intpday, & tsin(1,j,1,ndyn1), tsin(1,j,1,ndyn2), ts(1,j,1), plonl ) end if if( vid(37) > 0 ) then call intp2d( 1, 0., dynday2, intpday, & ts_avg_in(1,j,1,ndyn1), ts_avg_in(1,j,1,ndyn2), ts_avg(1,j,1), plonl ) end if if( vid(38) > 0 ) then call intp2d( 1, 0., dynday2, intpday, & fsds_avg_in(1,j,1,ndyn1), fsds_avg_in(1,j,1,ndyn2), fsds_avg(1,j,1), plonl ) end if end do !$omp end do !$omp end parallel end if end subroutine interp1 subroutine interp2( ncdate, ncsec, ps, u, v, & w, cgs, kvh, ts, taux, & tauy, hflx, qflx, qrad, omga, & pdmfuen, pdmfude, pdmfden, pdmfdde, & snow, fsds, soilw, plonl, platl, pplon ) !--------------------------------------------------------------------------- ! ... interpolate the fields whose values are required at ! the midpoint of a timestep. !--------------------------------------------------------------------------- use mo_mpi, only : masternode, lastnode use mo_calendar, only : diffdat use mo_control, only : xactive_drydep, xactive_emissions implicit none !--------------------------------------------------------------------------- ! ... dummy arguments !--------------------------------------------------------------------------- integer, intent(in) :: ncdate ! interpolate the input data to (ncdate,ncsec) integer, intent(in) :: ncsec integer, intent(in) :: plonl integer, intent(in) :: platl integer, intent(in) :: pplon real, intent(out) :: u(plonl,plev,-3:platl+4,pplon) ! u wind (time interpolant) real, intent(out) :: v(plonl,plev,-2:platl+3,pplon) ! v wind (time interpolant) real, intent(out) :: ps(plonl,-3:platl+4,pplon) ! interpolated surface pressure real, intent(out) :: w(plonl,plevp,platl,pplon) ! vertical wind (time interpolant) real, intent(out) :: omga(plonl,plev,platl,pplon) ! vertical velocity( pa/s ) real, intent(out) :: cgs(plonl,plevp,platl,pplon) ! interpolated counter-gradient coefficient real, intent(out) :: kvh(plonl,plevp,platl,pplon) ! interpolated vertical diffusion coefficient real, intent(out) :: ts(plonl,platl,pplon) ! surface temperature real, intent(out) :: taux(plonl,platl,pplon) ! stress in x direction real, intent(out) :: tauy(plonl,platl,pplon) ! stress in y direction real, intent(out) :: hflx(plonl,platl,pplon) ! surface sensible heat flux (w/m2) real, intent(out) :: qflx(plonl,platl,pplon) ! surface water vapor flux (kg/m^2/s) real, intent(out) :: snow(plonl,platl,pplon) ! snow height (m) real, intent(out) :: fsds(plonl,platl,pplon) ! downward shorwave radiation (w/m^2) real, intent(inout) :: soilw(plonl,platl,pplon) ! soil moisture fraction #ifdef DI_CONV_CCM real, dimension(plonl,plev,platl,pplon), intent(out) :: & #else real, dimension(1,1,1,1), intent(out) :: & #endif qrad ! radiative heating tendency (k/s) #if defined DI_TIEDTKE || defined AR_TIEDTKE real, dimension(plonl,plev,platl,pplon), intent(out) :: & #else real, dimension(1,1,1,1), intent(out) :: & #endif pdmfuen, & ! entrainment into updraft pdmfude, & ! detrainment into updraft pdmfden, & ! entrainment into downdraft pdmfdde ! detrainment into downdraft !--------------------------------------------------------------------------- ! ... local variables !--------------------------------------------------------------------------- integer :: & ip, & ! index j ! index integer, dimension(3) :: & jl, & ! index ju ! index real :: & dynday2, & ! (dyndate2,dynsec2) - (dyndate1,dynsec1) in days intpday ! (ncdate,ncsec) - (dyndate1,dynsec1) in days dynday2 = diffdat( dyndate1, dynsec1, dyndate2, dynsec2 ) intpday = diffdat( dyndate1, dynsec1, ncdate, ncsec ) #ifdef DEBUG write(*,*) 'interp2: timing info' write(*,*) 'dynday2 = ',dynday2 write(*,*) 'intpday = ',intpday #endif #ifdef USE_MPI jl(1) = -3 ju(1) = platl + 4 jl(2) = -2 ju(2) = platl + 3 jl(3) = -3 ju(3) = platl + 4 if( masternode ) then jl(:2) = 2 jl(3) = 1 end if if( lastnode ) then ju(:2) = platl - 1 ju(3) = platl end if #else jl(:2) = 2 jl(3) = 1 ju(:2) = platl - 1 ju(3) = platl #endif if( pplon > 1 ) then !$omp parallel do private( j, ip ) do ip = 1,pplon do j = jl(1),ju(1) call intp2d( plev, 0., dynday2, intpday, & uin(1,1,j,ip,ndyn1), uin(1,1,j,ip,ndyn2), u(1,1,j,ip), plonl ) end do do j = jl(2),ju(2) call intp2d( plev, 0., dynday2, intpday, & vin(1,1,j,ip,ndyn1), vin(1,1,j,ip,ndyn2), v(1,1,j,ip), plonl ) end do do j = jl(3),ju(3) call intp2d( 1, 0., dynday2, intpday, & psin(1,j,ip,ndyn1), psin(1,j,ip,ndyn2), ps(1,j,ip), plonl ) end do do j = 1,platl if( divdiff ) then call intp2d( 1, 0., dynday2, intpday, & tsin(1,j,ip,ndyn1), tsin(1,j,ip,ndyn2), ts(1,j,ip), plonl ) call intp2d( 1, 0., dynday2, intpday, & tauxin(1,j,ip,ndyn1), tauxin(1,j,ip,ndyn2), taux(1,j,ip), plonl ) call intp2d( 1, 0., dynday2, intpday, & tauyin(1,j,ip,ndyn1), tauyin(1,j,ip,ndyn2), tauy(1,j,ip), plonl ) call intp2d( 1, 0., dynday2, intpday, & hflxin(1,j,ip,ndyn1), hflxin(1,j,ip,ndyn2), hflx(1,j,ip), plonl ) end if if( divdiff .or. ditiedtke ) then call intp2d( 1, 0., dynday2, intpday, & qflxin(1,j,ip,ndyn1), qflxin(1,j,ip,ndyn2), qflx(1,j,ip), plonl ) end if if( mandatory(33) ) then call intp2d( 1, 0., dynday2, intpday, & snowin(1,j,ip,ndyn1), snowin(1,j,ip,ndyn2), snow(1,j,ip), plonl ) end if if( mandatory(34) ) then call intp2d( 1, 0., dynday2, intpday, & fsdsin(1,j,ip,ndyn1), fsdsin(1,j,ip,ndyn2), fsds(1,j,ip), plonl ) end if if( xactive_drydep .and. dyn_soilw ) then call intp2d( 1, 0., dynday2, intpday, & soilwin(1,j,ip,ndyn1), soilwin(1,j,ip,ndyn2), soilw(1,j,ip), plonl ) end if end do end do !$omp end parallel do else !$omp parallel private( j ) !$omp do do j = jl(1),ju(1) call intp2d( plev, 0., dynday2, intpday, & uin(1,1,j,1,ndyn1), uin(1,1,j,1,ndyn2), u(1,1,j,1), plonl ) end do !$omp end do !$omp do do j = jl(2),ju(2) call intp2d( plev, 0., dynday2, intpday, & vin(1,1,j,1,ndyn1), vin(1,1,j,1,ndyn2), v(1,1,j,1), plonl ) end do !$omp end do !$omp do do j = jl(3),ju(3) call intp2d( 1, 0., dynday2, intpday, & psin(1,j,1,ndyn1), psin(1,j,1,ndyn2), ps(1,j,1), plonl ) end do !$omp end do !$omp do do j = 1,platl if( divdiff ) then call intp2d( 1, 0., dynday2, intpday, & tsin(1,j,1,ndyn1), tsin(1,j,1,ndyn2), ts(1,j,1), plonl ) call intp2d( 1, 0., dynday2, intpday, & tauxin(1,j,1,ndyn1), tauxin(1,j,1,ndyn2), taux(1,j,1), plonl ) call intp2d( 1, 0., dynday2, intpday, & tauyin(1,j,1,ndyn1), tauyin(1,j,1,ndyn2), tauy(1,j,1), plonl ) call intp2d( 1, 0., dynday2, intpday, & hflxin(1,j,1,ndyn1), hflxin(1,j,1,ndyn2), hflx(1,j,1), plonl ) end if if( divdiff .or. ditiedtke ) then call intp2d( 1, 0., dynday2, intpday, & qflxin(1,j,1,ndyn1), qflxin(1,j,1,ndyn2), qflx(1,j,1), plonl ) end if if( mandatory(33) ) then call intp2d( 1, 0., dynday2, intpday, & snowin(1,j,1,ndyn1), snowin(1,j,1,ndyn2), snow(1,j,1), plonl ) end if if( mandatory(34) ) then call intp2d( 1, 0., dynday2, intpday, & fsdsin(1,j,1,ndyn1), fsdsin(1,j,1,ndyn2), fsds(1,j,1), plonl ) end if if( xactive_drydep .and. dyn_soilw ) then call intp2d( 1, 0., dynday2, intpday, & soilwin(1,j,1,ndyn1), soilwin(1,j,1,ndyn2), soilw(1,j,1), plonl ) end if end do !$omp end do !$omp end parallel end if end subroutine interp2 subroutine nextdyn !----------------------------------------------------------------------- ! ... close [and remove] the current dynamics file. open the next ! dynamics file. acquire file from mass store if local version ! doesn''t exist. !----------------------------------------------------------------------- use netcdf use mo_file_utils, only : acquire_file, open_netcdf_file use mo_calendar, only : addsec2dat use mo_calendar, only : diffdat use mo_charutl, only : incstr use mo_control, only : end_time #ifdef USE_MPI use mo_mpi, only : masternode, mpi_comm_comp, mpi_success use mo_mpi, only : thisnode, mpi_integer, mpi_logical #else use mo_mpi, only : masternode, thisnode #endif implicit none !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- integer :: & i, istat, slen, attempts, & recdid, nrec integer, allocatable, dimension(:) :: & next_date, & ! dates in input file next_datesec ! seconds relative to date character(len=128) :: & ctmp character(len=168) :: & loc_file character(len=10) :: & ctime character(len=8) :: & cdate logical :: msread_done, there logical :: get_next_file !----------------------------------------------------------------------- ! ... function declarations !----------------------------------------------------------------------- integer :: fsystem !----------------------------------------------------------------------- ! ... close current dynamics file !----------------------------------------------------------------------- Master_node_only : & if( masternode ) then if( close_file ) then call handle_ncerr( nf_close( ncid ), 'nextdyn: error closing file ' // & trim(lpath) // trim(filename)) end if !----------------------------------------------------------------------- ! ... remove if requested !----------------------------------------------------------------------- if( rmdyn ) then loc_file = last_filename write(*,*) 'nextdyn: removing file = ',trim(lpath) // trim(loc_file) ctmp = 'rm -f ' // trim(lpath) // trim(loc_file) write(*,*) 'nextdyn: fsystem issuing command - ' write(*,*) trim(ctmp) istat = fsystem( trim(ctmp) ) end if #ifdef DEBUG write(*,*) thisnode,': nextdyn: checking for next dyn file, initialization = ',initialization #endif is_initialization : & if( .not. initialization ) then #ifdef NCAR msread_done = .false. do attempts = 1,10 !----------------------------------------------------------------------- ! ... check to see if dynamics file acquired by previous asynch msread is on local disk !----------------------------------------------------------------------- inquire( file = trim( lpath ) // trim( filename ), exist = there, iostat = istat ) #ifdef DEBUG write(*,*) 'nextdyn: checking for next dyn file, istat,there,attempts=',istat,there,attempts #endif if( istat /= 0 ) then write(*,*) 'nextdyn: inquire failed for file ',trim( lpath ) // trim( filename ) write(*,*) 'error code = ',istat call endrun end if if( there ) then msread_done = .true. exit else #ifdef DEBUG write(*,*) 'nextdyn: waiting for dyn file, starting sleep #',attempts #endif call sleep( 30 ) end if end do #ifdef DEBUG write(*,*) thisnode,': nextdyn: msread_done=',msread_done #endif if( .not. msread_done ) then call date_and_time( cdate, ctime ) write(*,*) 'nextdyn: failed to acquire file ',trim( rpath ) // trim(filename),' from mass store' write(*,*) ' at time ',cdate,' ',ctime call endrun end if #else inquire( file = trim( lpath ) // trim( filename ), exist = there, iostat = istat ) if( istat /= 0 ) then write(*,*) 'nextdyn: inquire failed for file ',trim( lpath ) // trim( filename ) write(*,*) 'error code = ',istat call endrun end if if( .not. there ) then call date_and_time( cdate, ctime ) write(*,*) 'nextdyn: file ',trim( rpath ) // trim(filename),' not found' write(*,*) ' at time ',cdate,' ',ctime call endrun end if #endif end if is_initialization end if Master_node_only !--------------------------------------------------------------------------- ! ... all compute nodes wait for master to check last msread !--------------------------------------------------------------------------- #ifdef USE_MPI call mpi_barrier( mpi_comm_comp, istat ) if( istat /= mpi_success ) then write(*,*) 'nextdyn: barrier failed; code = ',istat call endrun end if #endif is_masternode_a : & if( masternode ) then is_initialization_a : & if( .not. initialization ) then call handle_ncerr( nf_open( trim(lpath) // filename, nf_nowrite, ncid ), & 'nextdyn: error opening file ' // trim(lpath) // trim(filename) ) !----------------------------------------------------------------------- ! ... get time info !----------------------------------------------------------------------- call handle_ncerr( nf_inq_unlimdim( ncid, recdid ), & 'nextdyn: inquiring record dimension id' ) call handle_ncerr( nf_inq_dimlen( ncid, recdid, nrec ), & 'nextdyn: inquiring record dimension length' ) if( allocated( next_date ) ) then deallocate( next_date ) end if if( allocated( next_datesec ) ) then deallocate( next_datesec ) end if allocate( next_date(nrec), next_datesec(nrec), stat=istat ) if( istat /= 0 ) then write(*,*) 'nextdyn: failed to allocate next_date ... next_datsec; error = ',istat call endrun end if call handle_ncerr( nf_inq_varid( ncid, 'date', datevid(1) ), & 'nextdyn: date not found in dynamics input file' ) call handle_ncerr( nf_inq_varid( ncid, 'datesec', datevid(2) ), & 'nextdyn: datesec not found in dynamics input file' ) call handle_ncerr( nf_get_var_int( ncid, datevid(1), next_date ), & 'nextdyn: getting date' ) call handle_ncerr( nf_get_var_int( ncid, datevid(2), next_datesec ), & 'nextdyn: getting datesec' ) curr_ts = 0 endif is_initialization_a !----------------------------------------------------------------------- ! ... build filenames for next dynamics file. netcdf files must ! have a .nc extension. !----------------------------------------------------------------------- write(*,*) 'nextdyn: attempting to increment filename ',trim( filename ) if( .not. initialization ) then call incr_flnm( filename, next_date(nrec), next_datesec(nrec), next_datesec(2) - next_datesec(1) ) else call incr_flnm( filename, date(ntim), datesec(ntim), datesec(2) - datesec(1) ) end if write(*,*) 'nextdyn: next filename = ',trim( filename ) is_force_synch : & if( force_synch ) then ncid = open_netcdf_file( filename, lpath, rpath, masteronly=.true. ) active_filename = filename initialization = .false. curr_ts = 0 !----------------------------------------------------------------------- ! ... get time info !----------------------------------------------------------------------- if( allocated( next_date ) ) then deallocate( next_date ) end if if( allocated( next_datesec ) ) then deallocate( next_datesec ) end if call handle_ncerr( nf_inq_unlimdim( ncid, recdid ), & 'nextdyn: inquiring record dimension id' ) call handle_ncerr( nf_inq_dimlen( ncid, recdid, nrec ), & 'nextdyn: inquiring record dimension length' ) allocate( next_date(nrec), next_datesec(nrec), stat=istat ) if( istat /= 0 ) then write(*,*) 'nextdyn: failed to allocate next_date ... next_datsec; error = ',istat call endrun end if call handle_ncerr( nf_inq_varid( ncid, 'date', datevid(1) ), & 'nextdyn: date not found in dynamics input file' ) call handle_ncerr( nf_inq_varid( ncid, 'datesec', datevid(2) ), & 'nextdyn: datesec not found in dynamics input file' ) call handle_ncerr( nf_get_var_int( ncid, datevid(1), next_date ), & 'nextdyn: getting date' ) call handle_ncerr( nf_get_var_int( ncid, datevid(2), next_datesec ), & 'nextdyn: getting datesec' ) !----------------------------------------------------------------------- ! ... build filenames for next dynamics file. netcdf files must ! have a .nc extension. !----------------------------------------------------------------------- call incr_flnm( filename, next_date(nrec), next_datesec(nrec), next_datesec(2) - next_datesec(1) ) endif is_force_synch if( .not. initialization ) then get_next_file = diffdat( next_date(nrec), next_datesec(nrec), end_time%date, end_time%secs ) > 0. else get_next_file = diffdat( date(ntim), datesec(ntim), end_time%date, end_time%secs ) > 0. end if !----------------------------------------------------------------------- ! ... asynch acquire next dyn file !----------------------------------------------------------------------- if( get_next_file ) then write(*,*) thisnode,': nextdyn: acquiring next dyn file' istat = acquire_file( trim(lpath) // filename, trim(rpath) // filename, & lundyn, .false., cosbdyn, .true. ) if( istat /= 0 ) then write(*,*) 'nextdyn: acquire_file: returned ', istat call endrun endif endif is_initialization_b : & if( .not. initialization ) then !--------------------------------------------------------------------------- ! ... get variable ids for fields required to be in the input data !--------------------------------------------------------------------------- do i = 1,mxdynflds if( mandatory(i) ) then call handle_ncerr( nf_inq_varid( ncid, nam(i), vid(i) ), & 'nextdyn: '// trim(nam(i)) // & ' not found in dynamics input file' ) else istat = nf_inq_varid( ncid, nam(i), vid(i) ) if( istat == NF_NOERR ) then write(*,*) 'nextdyn: Found optional variable '// trim(nam(i)) // & ' in dynamics input file' else vid(i) = 0 end if end if end do !----------------------------------------------------------------------- ! ... update date info !----------------------------------------------------------------------- call handle_ncerr( nf_inq_unlimdim( ncid, recdid ), & 'nextdyn: inquiring record dimension id' ) call handle_ncerr( nf_inq_dimlen( ncid, recdid, nrec ), & 'nextdyn: inquiring record dimension length' ) if( nrec /= ntim ) then deallocate( date, datesec ) ntim = nrec allocate( date(ntim), datesec(ntim), stat=istat ) if( istat /= 0 ) then write(*,*) thisnode,': nextdyn: failed to allocate date,datesec; error = ',istat call endrun end if end if call handle_ncerr( nf_inq_varid( ncid, 'date', datevid(1) ), & 'nextdyn: date not found in dynamics input file' ) call handle_ncerr( nf_inq_varid( ncid, 'datesec', datevid(2) ), & 'nextdyn: datesec not found in dynamics input file' ) call handle_ncerr( nf_get_var_int( ncid, datevid(1), date ), & 'nextdyn: getting date' ) call handle_ncerr( nf_get_var_int( ncid, datevid(2), datesec ), & 'nextdyn: getting datesec' ) !----------------------------------------------------------------------- ! ... adjust dates by the specified offset !----------------------------------------------------------------------- do i = 1,ntim call addsec2dat( dynoffset, date(i), datesec(i) ) end do endif is_initialization_b endif is_masternode_a #ifdef USE_MPI call mpi_bcast( initialization, 1, mpi_logical, 0, mpi_comm_comp, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': nextdyn: failed to bcast initialization; error = ',istat call endrun endif is_initialization_c : & if( .not. initialization ) then call mpi_bcast( curr_ts, 1, mpi_integer, 0, mpi_comm_comp, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': nextdyn: failed to bcast curr_ts; error = ',istat call endrun endif call mpi_bcast( nrec, 1, mpi_integer, 0, mpi_comm_comp, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': nextdyn: failed to bcast nrec; error = ',istat call endrun endif if( .not. masternode .and. nrec /= ntim ) then deallocate( date, datesec ) ntim = nrec allocate( date(ntim), datesec(ntim), stat=istat ) if( istat /= 0 ) then write(*,*) thisnode,': nextdyn: failed to allocate date,datesec; error = ',istat call endrun end if end if call mpi_bcast( date, ntim, mpi_integer, 0, mpi_comm_comp, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': nextdyn: failed to bcast date; error = ',istat call endrun endif call mpi_bcast( datesec, ntim, mpi_integer, 0, mpi_comm_comp, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': nextdyn: failed to bcast datesec; error = ',istat call endrun endif call mpi_bcast( vid, mxdynflds, mpi_integer, 0, mpi_comm_comp, istat ) if( istat /= nf_noerr ) then write(*,*) thisnode,': nextdyn: failed to bcast vid; error = ',istat call endrun endif write(*,'(i3,'' : nextdyn: curr_ts,nrec,ntim = '',3i6)') curr_ts,nrec,ntim endif is_initialization_c #endif end subroutine nextdyn subroutine incr_flnm( filename, pdate, psec, delt ) !----------------------------------------------------------------------- ! ... Increment or decrement a date string withing a filename ! the filename date section is assumed to be of the form ! yyyy-dd-mm !----------------------------------------------------------------------- use mo_calendar, only : newdate, caldayr use mo_charutl, only : incstr implicit none !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, optional, intent(in) :: pdate ! present date (yyyymmdd) integer, optional, intent(in) :: psec ! present seconds withing date integer, optional, intent(in) :: delt ! dynamical dataset time sample increment character(len=*), intent(inout) :: filename ! present dynamical dataset filename !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- integer :: beg_date, end_date, pos, pos1, m, slen integer :: year, month, day, date, secs, istat integer :: times(4) real :: calday character(len=128) :: fn_new character(len=6) :: seconds character(len=6) :: num character(len=4) :: frmt pos = index( trim(filename),trim(dyn_flnm_prefix) ) if( pos /= 1 ) then write(*,*) 'incr_flnm: can not find flnm prefix' call endrun end if is_date_frmt : & if( sum( date_template%start(:) ) /= 0 ) then !----------------------------------------------------------------------- ! ... date string in filename !----------------------------------------------------------------------- if( .not. present( pdate ) .or. .not. present( psec ) .or. .not. present( delt ) ) then write(*,*) 'incr_flnm: date, datesec, and delt arguments must be present' call endrun end if fn_new(1:) = trim( dyn_flnm_prefix ) // trim( dyn_flnm_date_frmt ) secs = psec + delt if( secs >= 86400 ) then day = secs/86400 secs = mod( secs,86400 ) date = newdate( pdate, day ) else date = pdate end if year = date/10000 month = mod( date,10000 )/100 day = date - (10000*year + 100*month) times(1) = year times(2) = month times(3) = day times(4) = secs pos = len_trim( dyn_flnm_prefix ) do m = 1,4 slen = date_template%digits(m) if( slen > 0 ) then write(frmt,'(''(i'',i1,'')'')') slen+1 pos1 = pos + date_template%start(m) write(num,frmt) times(m) + int( 10**(date_template%digits(m)) ) fn_new(pos1:pos1+slen-1) = num(2:slen+1) end if end do pos = len_trim( fn_new ) + 1 slen = len_trim( dyn_flnm_suffix ) if( slen > 0 ) then pos = len_trim( fn_new ) + 1 fn_new(pos:) = trim( dyn_flnm_suffix ) end if else is_date_frmt !----------------------------------------------------------------------- ! ... ccm type filename !----------------------------------------------------------------------- pos = len_trim( filename ) fn_new = filename(:pos) write(*,*) 'filename = ',trim(fn_new) if( fn_new(pos-2:) == '.nc' ) then pos = pos - 3 end if istat = incstr( fn_new(:pos), 1 ) if( istat /= 0 ) then write(*,*) 'incr_flnm: incstr returned ', istat write(*,*) ' while trying to decrement ',trim( fn_new ) call endrun end if end if is_date_frmt write(*,*) 'incr_flnm: new filename = ',trim(fn_new) filename = trim(fn_new) end subroutine incr_flnm subroutine qdyn( ncdate, ncsec, lpth, rpth, tim_idx ) !----------------------------------------------------------------------- ! ... query the dynamics module. !----------------------------------------------------------------------- use mo_mpi, only : masternode use mo_calendar, only : diffdat implicit none !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: & ncdate, & ! current simulation date (yyyymmdd) ncsec ! current simulation seconds withing current date (s) character(len=*), intent(out) :: & lpth, & ! local path for dynamics data rpth ! remote path for dynamics data integer, intent(out) :: & tim_idx ! time index for most recently read data if( masternode ) then lpth = trim(lpath) // active_filename rpth = trim(rpath) // active_filename tim_idx = curr_ts end if end subroutine qdyn subroutine intp2d( nlev, t1, t2, tint, f1, & f2, fint, plonl ) !----------------------------------------------------------------------- ! ... linearly interpolate between f1(t1) and f2(t2) to fint(tint). !----------------------------------------------------------------------- implicit none !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: & nlev, & plonl real, intent(in) :: & t1, & ! time level of f1 t2, & ! time level of f2 tint ! interpolant time real, dimension(plonl,nlev), intent(in) :: & f1, & ! field at time t1 f2 ! field at time t2 real, intent(out) :: & fint(plonl,nlev) ! field at time tint !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- integer :: k real :: factor factor = (tint - t1)/(t2 - t1) do k = 1,nlev fint(:,k) = f1(:,k) + (f2(:,k) - f1(:,k))*factor end do end subroutine intp2d integer function lotim( cdate, csec ) !----------------------------------------------------------------------- ! ... return the index of the dynamics time sample that is the lower ! bound of the interval that contains the input date. if ! (cdate,csec) is earlier than the first time sample then 0 is ! returned. if (cdate,csec) is later than the last time sample then ! that index is returned. if (cdate,csec) is equal to the date of a ! dynamics time sample then that index is returned. !----------------------------------------------------------------------- use mo_calendar, only : diffdat implicit none !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: & cdate, & ! date in yyyymmdd csec ! seconds relative to date !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- integer :: i !----------------------------------------------------------------------- ! ... find latest date that is earlier than or equal to (date,sec) !----------------------------------------------------------------------- do i = 1,ntim if( diffdat( cdate, csec, date(i), datesec(i) ) > 0. ) then lotim = i - 1 exit end if if( i == ntim ) then lotim = ntim end if end do end function lotim subroutine bsslzr( bes, n ) !----------------------------------------------------------------------- ! ... return n zeros (or if n>50, approximate zeros), of the ! bessel function j0,in the array bes. the first 50 zeros ! will be given exactly, and the remaining zeros are computed ! by extrapolation, and therefore not exact. !----------------------------------------------------------------------- use mo_constants, only : pi implicit none !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: n ! number of zeros to return real, intent(out) :: bes(n) ! array containing zeros of j0 !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- integer :: j, nn ! loop indices real, save :: bz(50) = & ! table of first 50 zeros (/ 2.4048255577, 5.5200781103, & 8.6537279129, 11.7915344391, 14.9309177086, 18.0710639679, & 21.2116366299, 24.3524715308, 27.4934791320, 30.6346064684, & 33.7758202136, 36.9170983537, 40.0584257646, 43.1997917132, & 46.3411883717, 49.4826098974, 52.6240518411, 55.7655107550, & 58.9069839261, 62.0484691902, 65.1899648002, 68.3314693299, & 71.4729816036, 74.6145006437, 77.7560256304, 80.8975558711, & 84.0390907769, 87.1806298436, 90.3221726372, 93.4637187819, & 96.6052679510, 99.7468198587, 102.8883742542, 106.0299309165, & 109.1714896498, 112.3130502805, 115.4546126537, 118.5961766309, & 121.7377420880, 124.8793089132, 128.0208770059, 131.1624462752, & 134.3040166383, 137.4455880203, 140.5871603528, 143.7287335737, & 146.8703076258, 150.0118824570, 153.1534580192, 156.2950342685 /) nn = n if( n > 50 ) then bes(50) = bz(50) do j = 51,n bes(j) = bes(j-1) + pi end do nn = 49 end if do j = 1,nn bes(j) = bz(j) end do end subroutine bsslzr subroutine gauaw( a, w, k ) !----------------------------------------------------------------------- ! ... calculate sine of latitudes a(k) and weights w(k) for the gaussian ! quadrature. the algorithm is described in davis and rabinowitz, ! journal of research of the nbs, v 56, jan 1956. ! the zeros of the bessel function j0, which are obtained ! from bsslzr, are used as a first guess for the abscissa. !----------------------------------------------------------------------- use mo_constants, only : pi implicit none !----------------------------------------------------------------------- ! ... dummy arguments !----------------------------------------------------------------------- integer, intent(in) :: k ! number of latitudes pole to pole real, intent(out) :: a(k), & ! sine of latitudes w(k) ! weights !----------------------------------------------------------------------- ! ... local variables !----------------------------------------------------------------------- real, parameter :: eps = 1.e-13 real :: & c, & ! constant combination fk, & ! real k xz, & ! abscissa estimate pkm1, & ! | pkm2, & ! |-polynomials pkmrk, & ! | pk, & ! | sp, & ! current iteration latitude increment avsp, & ! |sp| fn ! real n integer :: & kk, & ! k/2 (number of latitudes in hemisphere) is, & ! latitude index iter, & ! iteration counter n, l ! indices !----------------------------------------------------------------------- ! the value eps, used for convergence tests in the iterations, ! can be changed. newton iteration is used to find the abscissas. !----------------------------------------------------------------------- c = (1. - (2./pi)**2)*.25 fk = k kk = k/2 call bsslzr( a, kk ) do is = 1,kk xz = cos( a(is)/sqrt( (fk + .5)**2 + c ) ) !----------------------------------------------------------------------- ! ... this is the first approximation to xz !----------------------------------------------------------------------- do iter = 1,10 pkm2 = 1. pkm1 = xz !----------------------------------------------------------------------- ! ... computation of the legendre polynomial !----------------------------------------------------------------------- do n = 2,k fn = n pk = ((2.*fn - 1.)*xz*pkm1 - (fn - 1.)*pkm2)/fn pkm2 = pkm1 pkm1 = pk end do pkm1 = pkm2 pkmrk = (fk*(pkm1 - xz*pk))/(1. - xz**2) sp = pk/pkmrk xz = xz - sp avsp = abs(sp) if( avsp <= eps ) then exit end if end do if( iter > 10 ) then !----------------------------------------------------------------------- ! ... error exit !----------------------------------------------------------------------- write(*,*) 'gauaw: failed to convergence in 10 iterations' call endrun end if a(is) = xz w(is) = (2.*(1. - xz**2))/(fk*pkm1)**2 end do if( k /= kk*2 ) then !----------------------------------------------------------------------- ! ... for odd k computation of weight at the equator !----------------------------------------------------------------------- a(kk+1) = 0. pk = 2./fk**2 do n = 2,k,2 fn = n pk = pk*fn**2/(fn - 1.)**2 end do w(kk+1) = pk end if !----------------------------------------------------------------------- ! ... complete the sets of abscissas and weights, using the symmetry. !----------------------------------------------------------------------- do n = 1,kk l = k + 1 - n a(l) = -a(n) w(l) = w(n) end do end subroutine gauaw subroutine date_cracker( date_tmp ) !-------------------------------------------------------------------- ! ... crack dyn flnm date template !-------------------------------------------------------------------- implicit none !-------------------------------------------------------------------- ! ... dummy arguments !-------------------------------------------------------------------- character(len=*), intent(in) :: date_tmp !-------------------------------------------------------------------- ! ... local variables !-------------------------------------------------------------------- integer :: m character(len=1) :: date_chars(4) = (/ 'Y', 'M', 'D', 'S' /) !-------------------------------------------------------------------- ! ... look for year, month, and day in template !-------------------------------------------------------------------- do m = 1,4 date_template%start(m) = index( trim(date_tmp),date_chars(m) ) end do write(*,*) ' ' write(*,'('' date_cracker: year,mnth,day,sec start index = '',4i5)') date_template%start(:) write(*,*) ' ' !-------------------------------------------------------------------- ! ... now get the digit count !-------------------------------------------------------------------- do m = 1,4 if( date_template%start(m) > 0 ) then date_template%digits(m) = index( trim(date_tmp), date_chars(m), back=.true. ) date_template%digits(m) = date_template%digits(m) - date_template%start(m) + 1 else date_template%digits(m) = 0 end if end do write(*,*) ' ' write(*,'('' date_cracker: year,mnth,day,sec digit count = '',4i5)') date_template%digits(:) write(*,*) ' ' end subroutine date_cracker end module mo_dyninp