!!zhx merged the bug fix: IOP Profile bug fixes for the single column model (SCM) simulations #2500 module iop !----------------------------------------------------------------------- !BOP ! ! !MODULE: iop ! ! !DESCRIPTION: ! iop specific routines ! ! !USES: ! use shr_kind_mod, only: r8 => shr_kind_r8 use shr_scam_mod, only: shr_scam_GetCloseLatLon use constituents, only: readtrace, cnst_get_ind, pcnst, cnst_name use string_utils, only: to_lower use pmgrid use prognostics use time_manager, only: timemgr_init, get_curr_date, get_curr_calday,& get_nstep,get_start_date,timemgr_time_inc use cam_abortutils, only: endrun use scamMod use wrap_nf use cam_logfile, only: iulog use phys_control, only: phys_getopts use eul_control_mod,only: eul_nsplit ! ! !PUBLIC TYPES: implicit none private real(r8), allocatable, target :: dqfx3sav(:,:,:,:) real(r8), allocatable, target :: t2sav(:,:,:) real(r8), allocatable, target :: divq3dsav(:,:,:,:) real(r8), allocatable, target :: divt3dsav(:,:,:) real(r8), allocatable, target :: betasav(:) integer :: closelatidx,closelonidx,latid,lonid,levid,timeid real(r8):: closelat,closelon ! ! !PUBLIC MEMBER FUNCTIONS: public :: init_iop_fields public :: readiopdata ! read iop boundary data public :: setiopupdate ! find index in iopboundary data for current time ! public :: scam_use_iop_srf ! !PUBLIC DATA: public betasav, & dqfx3sav, divq3dsav, divt3dsav,t2sav ! ! !REVISION HISTORY: ! Created by John Truesdale ! !EOP ! ! !PRIVATE MEMBER FUNCTIONS: !----------------------------------------------------------------------- contains subroutine init_iop_fields() !------------------------------------------------------------------------------ ! Coupler for converting dynamics output variables into physics input variables ! also writes dynamics variables (on physics grid) to history file !------------------------------------------------------------------------------ implicit none !----------------------------------------------------------------------- if (eul_nsplit>1) then call endrun('iop module cannot be used with eul_nsplit>1') endif if(.not.allocated(betasav)) then allocate (betasav(beglat:endlat)) betasav(:)=0._r8 endif if(.not.allocated(dqfx3sav)) then allocate (dqfx3sav(plon,plev,pcnst,beglat:endlat)) dqfx3sav(:,:,:,:)=0._r8 endif if(.not.allocated(divq3dsav)) then allocate (divq3dsav(plon,plev,pcnst,beglat:endlat)) divq3dsav(:,:,:,:)=0._r8 endif if(.not.allocated(divt3dsav)) then allocate (divt3dsav(plon,plev,beglat:endlat)) divt3dsav(:,:,:)=0._r8 endif if(.not.allocated(t2sav)) then allocate (t2sav(plon,plev,beglat:endlat)) ! temp tendency t2sav(:,:,:)=0._r8 endif end subroutine init_iop_fields subroutine readiopdata( iop_update_surface ) !----------------------------------------------------------------------- ! ! Open and read netCDF file containing initial IOP conditions ! !---------------------------Code history-------------------------------- ! ! Written by J. Truesdale August, 1996, revised January, 1998 ! !----------------------------------------------------------------------- use comsrf use ppgrid, only: begchunk, endchunk use phys_grid, only: clat_p use commap use string_utils, only: to_lower use getinterpnetcdfdata use shr_sys_mod, only: shr_sys_flush use hycoef, only : hyam, hybm use eul_control_mod use error_messages, only : handle_ncerr use netcdf !----------------------------------------------------------------------- implicit none #if ( defined RS6000 ) implicit automatic ( a-z ) #endif !------------------------------Locals----------------------------------- ! logical, intent(in) :: iop_update_surface integer NCID, status integer time_dimID, lev_dimID,lev_varID,mod_dimID,& mod_varID,sps_varID,sps_dimID integer tsec_varID, bdate_varID,varid integer i,j integer nlev, nmod, nsps integer total_levs integer bdate, ntime integer, allocatable :: tsec(:) integer k, m integer icldliq,icldice integer inumliq,inumice logical have_srf ! value at surface is available logical fill_ends ! logical have_cnst(pcnst) real(r8), allocatable :: dplevs( : ) integer, allocatable :: dmods( : ) real(r8) dummy real(r8) lat,xlat real(r8) srf(1) ! value at surface real(r8) pmid(plev) ! pressure at model levels (time n) real(r8) pint(plevp) ! pressure at model interfaces (n ) real(r8) pdel(plev) ! pdel(k) = pint (k+1)-pint (k) real(r8) weight real(r8) tmpdata(1) real(r8) coldata(plev) integer strt4(4),cnt4(4) character(len=16) :: lowername fill_ends= .false. ! ! Open IOP dataset ! call handle_ncerr( nf90_open (iopfile, 0, ncid),& 'readiopdata.F90', __LINE__) ! ! if the dataset is a CAM generated dataset set use_camiop to true ! CAM IOP datasets have a global attribute called CAM_GENERATED_IOP ! if ( nf90_inquire_attribute( ncid, NF90_GLOBAL, 'CAM_GENERATED_FORCING', attnum=i ).EQ. NF90_NOERR ) then use_camiop = .true. else use_camiop = .false. endif !===================================================================== ! ! Read time variables status = nf90_inq_dimid (ncid, 'time', time_dimID ) if (status /= NF90_NOERR) then status = nf90_inq_dimid (ncid, 'tsec', time_dimID ) if (status /= NF90_NOERR) then write(iulog,* )'ERROR - readiopdata.F:Could not find dimension ID for time/tsec' status = NF90_CLOSE ( ncid ) call endrun end if end if call handle_ncerr( nf90_inquire_dimension( ncid, time_dimID, len=ntime ),& 'readiopdata.F90', __LINE__) allocate(tsec(ntime)) status = nf90_inq_varid (ncid, 'tsec', tsec_varID ) call handle_ncerr( nf90_get_var (ncid, tsec_varID, tsec),& 'readiopdata.F90', __LINE__) status = nf90_inq_varid (ncid, 'nbdate', bdate_varID ) if (status /= NF90_NOERR) then status = nf90_inq_varid (ncid, 'bdate', bdate_varID ) if (status /= NF90_NOERR) then write(iulog,* )'ERROR - readiopdata.F:Could not find variable ID for bdate' status = NF90_CLOSE ( ncid ) call endrun end if end if call handle_ncerr( nf90_get_var (ncid, bdate_varID, bdate),& 'readiopdata.F90', __LINE__) ! !====================================================== ! read level data ! status = NF90_INQ_DIMID( ncid, 'lev', lev_dimID ) if ( status .ne. nf90_noerr ) then write(iulog,* )'ERROR - readiopdata.F:Could not find variable dim ID for lev' status = NF90_CLOSE ( ncid ) return end if call handle_ncerr( nf90_inquire_dimension( ncid, lev_dimID, len=nlev ),& 'readiopdata.F90', __LINE__) allocate(dplevs(nlev+1)) status = NF90_INQ_VARID( ncid, 'lev', lev_varID ) if ( status .ne. nf90_noerr ) then write(iulog,* )'ERROR - readiopdata.F:Could not find variable ID for lev' status = NF90_CLOSE ( ncid ) return end if call handle_ncerr( nf90_get_var (ncid, lev_varID, dplevs(:nlev)),& 'readiopdata.F90', __LINE__) ! ===================================================== ! read observed aersol data if(scm_observed_aero .and. .not. iop_update_surface) then status = NF90_INQ_DIMID( ncid, 'mod', mod_dimID ) if ( status .ne. nf90_noerr ) then write(iulog,* )'ERROR - readiopdata.F:Could not find variable dim ID for lev' status = NF90_CLOSE ( ncid ) return end if call handle_ncerr( nf90_inquire_dimension( ncid, mod_dimID, len=nmod ),& 'readiopdata.F90', __LINE__) status = NF90_INQ_DIMID( ncid, 'sps', sps_dimID ) if ( status .ne. nf90_noerr ) then write(iulog,* )'ERROR - readiopdata.F:Could not find variable dim ID for sps' status = NF90_CLOSE ( ncid ) return end if call handle_ncerr( nf90_inquire_dimension( ncid, sps_dimID, len=nsps ),& 'readiopdata.F90', __LINE__) if (.not.allocated(dmods)) then allocate(dmods(nmod)) dmods=-999 end if if (.not.allocated(scm_num)) then allocate(scm_num(nmod)) scm_num= 1.0e30_R8 end if if (.not.allocated(scm_dgnum)) then allocate(scm_dgnum(nmod)) scm_dgnum= 1.0e30_R8 end if if (.not.allocated(scm_std)) then allocate(scm_std(nmod)) scm_std= 1.0e30_R8 end if if (.not.allocated(scm_div)) then allocate(scm_div(nmod,nsps)) scm_div= 1.0e30_R8 end if status = NF90_INQ_VARID( ncid, 'mod', mod_varID ) if ( status .ne. nf90_noerr ) then write(iulog,* )'ERROR - readiopdata.F:Could not find variable ID for mode' status = NF90_CLOSE ( ncid ) return end if call handle_ncerr( nf90_get_var (ncid, mod_varID, dmods(:nmod)),& 'readiopdata.F90', __LINE__) status = NF90_INQ_VARID( ncid, 'scm_num', mod_varID ) if ( status .ne. nf90_noerr ) then write(iulog,* )'ERROR - readiopdata.F:Could not find variable ID for mode' status = NF90_CLOSE ( ncid ) return end if call handle_ncerr( nf90_get_var (ncid, mod_varID, scm_num(:nmod)),& 'readiopdata.F90', __LINE__) status = NF90_INQ_VARID( ncid, 'scm_diam', mod_varID ) if ( status .ne. nf90_noerr ) then write(iulog,* )'ERROR - readiopdata.F:Could not find variable ID for mode' status = NF90_CLOSE ( ncid ) return end if call handle_ncerr( nf90_get_var (ncid, mod_varID, scm_dgnum(:nmod)),& 'readiopdata.F90', __LINE__) status = NF90_INQ_VARID( ncid, 'scm_std', mod_varID ) if ( status .ne. nf90_noerr ) then write(iulog,* )'ERROR - readiopdata.F:Could not find variable ID for mode' status = NF90_CLOSE ( ncid ) return end if call handle_ncerr( nf90_get_var (ncid, mod_varID, scm_std(:nmod)),& 'readiopdata.F90', __LINE__) status = NF90_INQ_VARID( ncid, 'scm_accum_div', sps_varID ) if ( status .ne. nf90_noerr ) then write(iulog,* )'ERROR - readiopdata.F:Could not find variable ID for mode' status = NF90_CLOSE ( ncid ) return end if call handle_ncerr( nf90_get_var (ncid, sps_varID, scm_div(1,:nsps)),& 'readiopdata.F90', __LINE__) status = NF90_INQ_VARID( ncid, 'scm_aitken_div', sps_varID ) if ( status .ne. nf90_noerr ) then write(iulog,* )'ERROR - readiopdata.F:Could not find variable ID for mode' status = NF90_CLOSE ( ncid ) return end if call handle_ncerr( nf90_get_var (ncid, sps_varID, scm_div(2,:nsps)),& 'readiopdata.F90', __LINE__) status = NF90_INQ_VARID( ncid, 'scm_coarse_div', sps_varID ) if ( status .ne. nf90_noerr ) then write(iulog,* )'ERROR - readiopdata.F:Could not find variable ID for mode' status = NF90_CLOSE ( ncid ) return end if call handle_ncerr( nf90_get_var (ncid, sps_varID, scm_div(3,:nsps)),& 'readiopdata.F90', __LINE__) endif !scm_observed_aero !====================================================================== ! !CAM generated forcing already has pressure on millibars ! if (.not. use_camiop) then ! ! convert pressure to millibars ( lev is expressed in pascals in iop datasets ) ! do i=1,nlev dplevs( i ) = dplevs( i )/100._r8 end do endif call shr_scam_GetCloseLatLon(ncid,scmlat,scmlon,closelat,closelon,closelatidx,closelonidx) lonid = 0 latid = 0 levid = 0 timeid = 0 call wrap_inq_dimid(ncid, 'lat', latid) call wrap_inq_dimid(ncid, 'lon', lonid) call wrap_inq_dimid(ncid, 'lev', levid) call wrap_inq_dimid(ncid, 'time', timeid) strt4(1) = closelonidx strt4(2) = closelatidx strt4(3) = iopTimeIdx strt4(4) = 1 cnt4(1) = 1 cnt4(2) = 1 cnt4(3) = 1 cnt4(4) = 1 if (.not. iop_update_surface) then status = nf90_inq_varid( ncid, 'Ps', varid ) if ( status .ne. nf90_noerr ) then have_ps = .false. write(iulog,*)'Could not find variable Ps' if ( .not. use_userdata ) then status = NF90_CLOSE( ncid ) return else if ( get_nstep() .eq. 0 ) write(iulog,*) 'Using value from Analysis Dataset' endif else status = nf90_get_var(ncid, varid, ps(1,1,n3), strt4) have_ps = .true. endif ! for reproducing CAM output don't do interpolation. ! the most expedient way of doing this is to set ! the dataset pressure levels to the current ! scam model levels if ( use_camiop ) then do i = 1, plev dplevs( i ) = 1000.0_r8 * hyam( i ) + ps(1,1,n3) * hybm( i ) / 100.0_r8 end do endif ! add the surface pressure to the pressure level data, so that ! surface boundary condition will be set properly, ! making sure that it is the highest pressure in the array. ! total_levs = nlev+1 dplevs(nlev+1) = ps(1,1,n3)/100 ! ps is expressed in pascals do i= nlev, 1, -1 if ( dplevs(i) .GT. ps(1,1,n3)/100.0_r8) then total_levs = i dplevs(i) = ps(1,1,n3)/100 end if end do if (.not. use_camiop ) then nlev = total_levs endif if ( nlev .eq. 1 ) then write(iulog,*) 'Error - Readiopdata.F: Ps too low!' return endif !===================================================================== status = nf90_inq_varid( ncid, 'Tsair', varid ) if ( status .ne. nf90_noerr ) then have_tsair = .false. else call wrap_get_vara_realx (ncid,varid,strt4,cnt4,tsair) have_tsair = .true. endif ! ! read in Tobs For cam generated iop readin small t to avoid confusion ! with capital T defined in cam ! if ( get_nstep() .le. 1 ) then tobs(:)=t3(1,:,1,1) else tobs(:)= t3(1,:,1,n3) endif !!!!!!!force fill_end to be .true in getinterpncdata () for temperature !!!!!!!!+zhx if ( use_camiop ) then call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx,'t', have_tsair, & !-zhx tsair(1), fill_ends, & tsair(1), .true., & !+zhx dplevs, nlev,ps(1,1,n3),tobs, status ) else call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx,'T', have_tsair, & !-zhx tsair(1), fill_ends, & tsair(1), .true., & !+zhx dplevs, nlev,ps(1,1,n3), tobs, status ) endif if ( status .ne. nf90_noerr ) then have_t = .false. write(iulog,*)'Could not find variable T' if ( .not. use_userdata ) then status = NF90_CLOSE( ncid ) return else write(iulog,*) 'Using value from Analysis Dataset' endif ! ! set T3 to Tobs on first time step ! else have_t = .true. if (.not.use_camiop .and. get_nstep() .eq. 0) then do i=1, PLEV t3(1,i,1,n3)=tobs(i) ! set t to tobs at first time step end do endif endif status = nf90_inq_varid( ncid, 'qsrf', varid ) if ( status .ne. nf90_noerr ) then have_srf = .false. else status = nf90_get_var(ncid, varid, srf(1), strt4) have_srf = .true. endif ! if ( get_nstep() .le. 1 ) then qobs(:)= q3(1,:,1,1,1) else qobs(:)= q3(1,:,1,1,n3) endif !!!!!!!force fill_end to be .true in getinterpncdata () for humidity!!!!!!!!+zhx call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'q', have_srf, & !-zhx srf(1), fill_ends, & srf(1), .true., & !+zhx dplevs, nlev,ps(1,1,n3), qobs, status ) if ( status .ne. nf90_noerr ) then have_q = .false. write(iulog,*)'Could not find variable q' if ( .not. use_userdata ) then status = nf90_close( ncid ) return else write(iulog,*) 'Using values from Analysis Dataset' endif else if (.not. use_camiop .and. get_nstep() .eq. 0) then do i=1, PLEV q3(1,i,1,1,n3)=qobs(i) end do have_q = .true. endif endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'cld', .false., & dummy, fill_ends, dplevs, nlev,ps(1,1,n3), cldobs, status ) if ( status .ne. nf90_noerr ) then have_cld = .false. else have_cld = .true. endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'clwp', .false., & dummy, fill_ends, dplevs, nlev,ps(1,1,n3), clwpobs, status ) if ( status .ne. nf90_noerr ) then have_clwp = .false. else have_clwp = .true. endif ! ! read divq (horizontal advection) ! status = nf90_inq_varid( ncid, 'divqsrf', varid ) if ( status .ne. nf90_noerr ) then have_srf = .false. else status = nf90_get_var(ncid, varid, srf(1), strt4) have_srf = .true. endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, & 'divq', have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), divq(:,1), status ) if ( status .ne. nf90_noerr ) then have_divq = .false. else have_divq = .true. endif ! ! read vertdivq if available ! status = nf90_inq_varid( ncid, 'vertdivqsrf', varid ) if ( status .ne. nf90_noerr ) then have_srf = .false. else status = nf90_get_var(ncid, varid, srf(1), strt4) have_srf = .true. endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'vertdivq', & have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), vertdivq(:,1), status ) if ( status .ne. nf90_noerr ) then have_vertdivq = .false. else have_vertdivq = .true. endif status = nf90_inq_varid( ncid, 'vertdivqsrf', varid ) if ( status .ne. nf90_noerr ) then have_srf = .false. else status = nf90_get_var(ncid, varid, srf(1), strt4) have_srf = .true. endif ! ! add calls to get dynamics tendencies for all prognostic consts ! do m = 1, pcnst call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, trim(cnst_name(m))//'_dten', & have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), divq3d(:,m), status ) if ( status .ne. nf90_noerr ) then have_cnst(m) = .false. divq3d(1:,m)=0._r8 else have_cnst(m) = .true. endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, trim(cnst_name(m))//'_dqfx', & have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), coldata, status ) if ( STATUS .NE. NF90_NOERR ) then dqfxcam=0._r8 else dqfxcam(1,:,m)=coldata(:) endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, trim(cnst_name(m))//'_alph', & have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), tmpdata, status ) if ( status .ne. nf90_noerr ) then ! have_cnst(m) = .false. alphacam(m)=0._r8 else alphacam(m)=tmpdata(1) ! have_cnst(m) = .true. endif end do call cnst_get_ind('NUMLIQ', inumliq, abort=.false.) if ( inumliq > 0 ) then have_srf = .false. call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'NUMLIQ', & have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), numliqobs, status ) if ( status .ne. nf90_noerr ) then have_numliq = .false. else have_numliq = .true. if ( get_nstep() .eq. 0) then do i=1, PLEV q3(1,i,inumliq,1,n3)=numliqobs(i) end do endif endif end if call cnst_get_ind('CLDLIQ', icldliq) have_srf = .false. call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'CLDLIQ', & have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), cldliqobs, status ) if ( status .ne. nf90_noerr ) then have_cldliq = .false. else have_cldliq = .true. if ( get_nstep() .eq. 0) then do i=1, PLEV q3(1,i,icldliq,1,n3)=cldliqobs(i) end do endif endif call cnst_get_ind('CLDICE', icldice) call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'CLDICE', & have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), cldiceobs, status ) if ( status .ne. nf90_noerr ) then have_cldice = .false. else have_cldice = .true. if ( get_nstep() .eq. 0) then do i=1, PLEV q3(1,i,icldice,1,n3)=cldiceobs(i) end do endif endif call cnst_get_ind('NUMICE', inumice, abort=.false.) if ( inumice > 0 ) then have_srf = .false. call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'NUMICE', & have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), numiceobs, status ) if ( status .ne. nf90_noerr ) then have_numice = .false. else have_numice = .true. if ( get_nstep() .eq. 0) then do i=1, PLEV q3(1,i,inumice,1,n3)=numiceobs(i) end do endif endif end if ! ! read divu (optional field) ! status = nf90_inq_varid( ncid, 'divusrf', varid ) if ( status .ne. nf90_noerr ) then have_srf = .false. else status = nf90_get_var(ncid, varid, srf(1), strt4) have_srf = .true. endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divu', & have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), divu, status ) if ( status .ne. nf90_noerr ) then have_divu = .false. else have_divu = .true. endif ! ! read divv (optional field) ! status = nf90_inq_varid( ncid, 'divvsrf', varid ) if ( status .ne. nf90_noerr ) then have_srf = .false. else status = nf90_get_var(ncid, varid, srf(1), strt4) have_srf = .true. endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divv', & have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), divv, status ) if ( status .ne. nf90_noerr ) then have_divv = .false. else have_divv = .true. endif ! ! read divt (optional field) ! status = nf90_inq_varid( ncid, 'divtsrf', varid ) if ( status .ne. nf90_noerr ) then have_srf = .false. else status = nf90_get_var(ncid, varid, srf(1), strt4) have_srf = .true. endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, & 'divT', have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), divt, status ) if ( status .ne. nf90_noerr ) then have_divt = .false. else have_divt = .true. endif ! ! read vertdivt if available ! status = nf90_inq_varid( ncid, 'vertdivTsrf', varid ) if ( status .ne. nf90_noerr ) then have_srf = .false. else status = nf90_get_var(ncid, varid, srf(1), strt4) have_srf = .true. endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'vertdivT', & have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), vertdivt, status ) if ( status .ne. nf90_noerr ) then have_vertdivt = .false. else have_vertdivt = .true. endif ! ! read divt3d (combined vertical/horizontal advection) ! (optional field) status = nf90_inq_varid( ncid, 'divT3dsrf', varid ) if ( status .ne. nf90_noerr ) then have_srf = .false. else status = nf90_get_var(ncid, varid, srf(1), strt4) have_srf = .true. endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'divT3d', & have_srf, srf(1), fill_ends, & dplevs, nlev,ps(1,1,n3), divt3d, status ) if ( status .ne. nf90_noerr ) then have_divt3d = .false. else have_divt3d = .true. endif status = nf90_inq_varid( ncid, 'Ptend', varid ) if ( status .ne. nf90_noerr ) then have_ptend = .false. write(iulog,*)'Could not find variable Ptend. Setting to zero' ptend = 0.0_r8 else status = nf90_get_var(ncid, varid, srf(1), strt4) have_ptend = .true. ptend= srf(1) endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, & 'omega', .true., ptend, fill_ends, & dplevs, nlev,ps(1,1,n3), wfld, status ) if ( status .ne. nf90_noerr ) then have_omega = .false. write(iulog,*)'Could not find variable omega' if ( .not. use_userdata ) then status = nf90_close( ncid ) return else write(iulog,*) 'Using value from Analysis Dataset' endif else have_omega = .true. endif call plevs0(1 ,plon ,plev ,ps(1,1,n3) ,pint,pmid ,pdel) call shr_sys_flush( iulog ) ! ! Build interface vector for the specified omega profile ! (weighted average in pressure of specified level values) ! wfldh(1) = 0.0_r8 do k=2,plev weight = (pint(k) - pmid(k-1))/(pmid(k) - pmid(k-1)) wfldh(k) = (1.0_r8 - weight)*wfld(k-1) + weight*wfld(k) end do wfldh(plevp) = 0.0_r8 status = nf90_inq_varid( ncid, 'usrf', varid ) if ( status .ne. nf90_noerr ) then have_srf = .false. else call wrap_get_vara_realx (ncid,varid,strt4,cnt4,srf) have_srf = .true. endif !!!!!!!force fill_end to be .true in getinterpncdata () for u-wind !!!!!!!!+zhx call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, & !-zhx 'u', have_srf, srf(1), fill_ends, & 'u', have_srf, srf(1), .true., & !+zhx dplevs, nlev,ps(1,1,n3), uobs, status ) if ( status .ne. nf90_noerr ) then have_u = .false. else have_u = .true. if (.not. use_camiop .and. get_nstep() .eq. 0 ) then do i=1, PLEV u3(1,i,1,n3) = uobs(i) ! set u to uobs at first time step end do endif endif status = nf90_inq_varid( ncid, 'vsrf', varid ) if ( status .ne. nf90_noerr ) then have_srf = .false. else call wrap_get_vara_realx (ncid,varid,strt4,cnt4,srf) have_srf = .true. endif !!!!!!!force fill_end to be .true in getinterpncdata () for v-wind !!!!!!!!+zhx call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, & !-zhx 'v', have_srf, srf(1), fill_ends, & 'v', have_srf, srf(1), .true., & !+zhx dplevs, nlev,ps(1,1,n3), vobs, status ) if ( status .ne. nf90_noerr ) then have_v = .false. else have_v = .true. if (.not. use_camiop .and. get_nstep() .eq. 0 ) then do i=1, PLEV v3(1,i,1,n3) = vobs(i) ! set u to uobs at first time step end do endif endif call shr_sys_flush( iulog ) status = nf90_inq_varid( ncid, 'Prec', varid ) if ( status .ne. nf90_noerr ) then have_prec = .false. else call wrap_get_vara_realx (ncid,varid,strt4,cnt4,precobs) have_prec = .true. endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'Q1', & .false., dummy, fill_ends, & ! datasets don't contain Q1 at surface dplevs, nlev,ps(1,1,n3), q1obs, status ) if ( status .ne. nf90_noerr ) then have_q1 = .false. else have_q1 = .true. endif call getinterpncdata( ncid, scmlat, scmlon, ioptimeidx, 'Q2', & .false., dummy, fill_ends, & ! datasets don't contain Q2 at surface dplevs, nlev,ps(1,1,n3), q1obs, status ) if ( status .ne. nf90_noerr ) then have_q2 = .false. else have_q2 = .true. endif ! ! fill in 3d forcing variables if we have both horizontal ! and vertical components, but not the 3d ! if ( .not. have_cnst(1) .and. have_divq .and. have_vertdivq ) then do k=1,plev do m=1,pcnst divq3d(k,m) = divq(k,m) + vertdivq(k,m) enddo enddo have_divq3d = .true. endif if ( .not. have_divt3d .and. have_divt .and. have_vertdivt ) then write(iulog,*) 'Don''t have divt3d - using divt and vertdivt' do k=1,plev divt3d(k) = divt(k) + vertdivt(k) enddo have_divt3d = .true. endif ! ! make sure that use_3dfrc flag is set to true if we only have ! 3d forcing available ! if ( .not. have_divt .or. .not. have_divq ) then use_3dfrc = .true. endif call shr_sys_flush( iulog ) status = nf90_inq_varid( ncid, 'CLAT', varid ) if ( status .eq. nf90_noerr ) then call wrap_get_vara_realx (ncid,varid,strt4,cnt4,clat) clat_p(1)=clat(1) latdeg(1) = clat(1)*45._r8/atan(1._r8) endif status = nf90_inq_varid( ncid, 'beta', varid ) if ( status .ne. nf90_noerr ) then betacam = 0._r8 else status = nf90_get_var(ncid, varid, srf(1), strt4) betacam=srf(1) endif status = nf90_inq_varid( ncid, 'fixmas', varid ) if ( status .ne. nf90_noerr ) then fixmascam=1.0_r8 else status = nf90_get_var(ncid, varid, srf(1), strt4) fixmascam=srf(1) endif else ! if read in surface information status = nf90_inq_varid( ncid, 'Tg', varid ) if (status .ne. nf90_noerr) then write(iulog,*)'Could not find variable Tg' if ( have_tsair ) then write(iulog,*) 'Using Tsair' tground = tsair ! use surface value from T field else write(iulog,*) 'Using T at lowest level' tground = t3(1,plev,1,n3) endif else call wrap_get_vara_realx (ncid,varid,strt4,cnt4,tground) have_tg = .true. endif ! Test for BOTH 'lhflx' and 'lh' without overwriting 'have_lhflx'. ! Analagous changes made for the surface heat flux status = nf90_inq_varid( ncid, 'lhflx', varid ) if ( status .ne. nf90_noerr ) then status = nf90_inq_varid( ncid, 'lh', varid ) if ( status .ne. nf90_noerr ) then have_lhflx = .false. else call wrap_get_vara_realx (ncid,varid,strt4,cnt4,lhflxobs) have_lhflx = .true. endif else call wrap_get_vara_realx (ncid,varid,strt4,cnt4,lhflxobs) have_lhflx = .true. endif status = nf90_inq_varid( ncid, 'shflx', varid ) if ( status .ne. nf90_noerr ) then status = nf90_inq_varid( ncid, 'sh', varid ) if ( status .ne. nf90_noerr ) then have_shflx = .false. else call wrap_get_vara_realx (ncid,varid,strt4,cnt4,shflxobs) have_shflx = .true. endif else call wrap_get_vara_realx (ncid,varid,strt4,cnt4,shflxobs) have_shflx = .true. endif endif call shr_sys_flush( iulog ) status = nf90_close( ncid ) call shr_sys_flush( iulog ) deallocate(dplevs,tsec) return end subroutine readiopdata subroutine setiopupdate !----------------------------------------------------------------------- ! ! Open and read netCDF file to extract time information ! !---------------------------Code history-------------------------------- ! ! Written by John Truesdale August, 1996 ! !----------------------------------------------------------------------- implicit none #if ( defined RS6000 ) implicit automatic (a-z) #endif !------------------------------Locals----------------------------------- integer NCID,i integer tsec_varID, time_dimID integer, allocatable :: tsec(:) integer ntime integer bdate, bdate_varID integer STATUS integer next_date, next_sec, last_date, last_sec integer :: ncsec,ncdate ! current time of day,date integer :: yr, mon, day ! year, month, and day component integer :: start_ymd,start_tod save tsec, ntime, bdate save last_date, last_sec !------------------------------------------------------------------------------ if ( get_nstep() .eq. 0 ) then ! ! Open IOP dataset ! STATUS = NF90_OPEN( iopfile, NF90_NOWRITE, NCID ) ! ! Read time (tsec) variable ! STATUS = NF90_INQ_VARID( NCID, 'tsec', tsec_varID ) if ( STATUS .NE. NF90_NOERR ) write(iulog,*)'ERROR - setiopupdate.F:', & 'Cant get variable ID for tsec' STATUS = NF90_INQ_VARID( NCID, 'bdate', bdate_varID ) if ( STATUS .NE. NF90_NOERR ) then STATUS = NF90_INQ_VARID( NCID, 'basedate', bdate_varID ) if ( STATUS .NE. NF90_NOERR ) & write(iulog,*)'ERROR - setiopupdate.F:Cant get variable ID for bdate' endif STATUS = NF90_INQ_DIMID( NCID, 'time', time_dimID ) if ( STATUS .NE. NF90_NOERR ) then STATUS = NF90_INQ_DIMID( NCID, 'tsec', time_dimID ) if ( STATUS .NE. NF90_NOERR ) then write(iulog,* )'ERROR - setiopupdate.F:Could not find variable dim ID for time' STATUS = NF90_CLOSE ( NCID ) return end if end if if ( STATUS .NE. NF90_NOERR ) & write(iulog,*)'ERROR - setiopupdate.F:Cant get variable dim ID for time' STATUS = NF90_INQUIRE_DIMENSION( NCID, time_dimID, len=ntime ) if ( STATUS .NE. NF90_NOERR )then write(iulog,*)'ERROR - setiopupdate.F:Cant get time dimlen' endif if (.not.allocated(tsec)) allocate(tsec(ntime)) STATUS = NF90_GET_VAR( NCID, tsec_varID, tsec ) if ( STATUS .NE. NF90_NOERR )then write(iulog,*)'ERROR - setiopupdate.F:Cant get variable tsec' endif STATUS = NF90_GET_VAR( NCID, bdate_varID, bdate ) if ( STATUS .NE. NF90_NOERR )then write(iulog,*)'ERROR - setiopupdate.F:Cant get variable bdate' endif ! Close the netCDF file STATUS = NF90_CLOSE( NCID ) ! ! determine the last date in the iop dataset ! call timemgr_time_inc(bdate, 0, last_date, last_sec, inc_s=tsec(ntime)) ! ! set the iop dataset index ! iopTimeIdx=0 do i=1,ntime ! set the first ioptimeidx call timemgr_time_inc(bdate, 0, next_date, next_sec, inc_s=tsec(i)) call get_start_date(yr,mon,day,start_tod) start_ymd = yr*10000 + mon*100 + day if ( start_ymd .gt. next_date .or. (start_ymd .eq. next_date & .and. start_tod .ge. next_sec)) then iopTimeIdx = i endif enddo call get_curr_date(yr,mon,day,ncsec) ncdate=yr*10000 + mon*100 + day if (iopTimeIdx == 0.or.iopTimeIdx .ge. ntime) then call timemgr_time_inc(bdate, 0, next_date, next_sec, inc_s=tsec(1)) write(iulog,*) 'Error::setiopupdate: Current model time does not fall within IOP period' write(iulog,*) ' Current CAM Date is ',ncdate,' and ',ncsec,' seconds' write(iulog,*) ' IOP start is ',next_date,' and ',next_sec,' seconds' write(iulog,*) ' IOP end is ',last_date,' and ',last_sec,' seconds' call endrun endif doiopupdate = .true. !------------------------------------------------------------------------------ ! Check if iop data needs to be updated and set doiopupdate accordingly !------------------------------------------------------------------------------ else ! endstep > 1 call timemgr_time_inc(bdate, 0, next_date, next_sec, inc_s=tsec(iopTimeIdx+1)) call get_curr_date(yr, mon, day, ncsec) ncdate = yr*10000 + mon*100 + day if ( ncdate .gt. next_date .or. (ncdate .eq. next_date & .and. ncsec .ge. next_sec)) then iopTimeIdx = iopTimeIdx + 1 doiopupdate = .true. #if DEBUG > 2 write(iulog,*) 'nstep = ',get_nstep() write(iulog,*) 'ncdate=',ncdate,' ncsec=',ncsec write(iulog,*) 'next_date=',next_date,' next_sec=',next_sec write(iulog,*)'******* do iop update' #endif else doiopupdate = .false. end if endif ! if (endstep .eq. 0 ) ! ! make sure we're ! not going past end of iop data ! if ( ncdate .gt. last_date .or. (ncdate .eq. last_date & .and. ncsec .gt. last_sec)) then if ( .not. use_userdata ) then write(iulog,*)'ERROR - setiopupdate.c:Reached the end of the time varient dataset' stop else doiopupdate = .false. end if endif #if DEBUG > 1 write(iulog,*)'iop time index = ' , ioptimeidx #endif return end subroutine setiopupdate end module iop