LCOV - code coverage report
Current view: | top level - physics/cam - boundarydata.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 181 | 366 | 49.5 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 5 | 7 | 71.4 %

              Line data    Source code
           1             : #define _FILE 'physics/cam1/boundarydata.F90 '
           2             : module boundarydata
           3             :   use shr_kind_mod,    only: r8 => shr_kind_r8
           4             :   use spmd_utils,      only: masterproc
           5             :   use ppgrid,          only: pcols, pver, begchunk, endchunk
           6             :   use physics_types,   only: physics_state
           7             :   use cam_abortutils,  only: endrun
           8             : #if ( defined SPMD )
           9             :   use mpishorthand,    only: mpicom, mpir8, mpiint
          10             : #endif
          11             :   use netcdf
          12             :   use error_messages,  only: handle_ncerr
          13             :   use cam_logfile,     only: iulog
          14             :   implicit none
          15             :   private
          16             :   integer, parameter :: ptrtim=12, ptrlon=1
          17             :
          18             :   type boundarydata_type
          19             :      integer           :: ncid
          20             :      integer           :: fieldcnt
          21             :      integer           :: nm
          22             :      integer           :: np
          23             :      integer           :: latsiz
          24             :      integer           :: levsiz
          25             :      integer           :: ncolsiz
          26             :      integer           :: timsiz
          27             :      integer           :: vertextrap
          28             :      logical           :: iszonal, isncol
          29             :      integer           :: ndims
          30             :      integer           :: thistimedim
          31             :      integer           :: psid
          32             :      integer  :: map(4)
          33             :      integer  :: dimids(4)
          34             :      integer,  pointer :: dataid(:) => null()
          35             :      integer,  pointer :: columnmap(:) => null()
          36             :      integer, pointer  :: start(:,:) => null()
          37             :      integer, pointer  :: count(:,:) => null()
          38             :      real(r8), pointer :: lat(:) => null()
          39             :      real(r8), pointer :: zi(:) => null()
          40             :      real(r8), pointer :: pin(:) => null()
          41             :      real(r8), pointer :: cdates(:) => null()
          42             :      real(r8), pointer :: fields(:,:,:,:,:) => null()
          43             :      real(r8), pointer :: datainst(:,:,:,:) => null()
          44             :      real(r8), pointer :: hybi(:) => null()
          45             :      real(r8), pointer :: ps(:,:,:) => null()
          46             :   end type boundarydata_type
          47             :
          48             :   public boundarydata_init
          49             :   public boundarydata_update
          50             :   public boundarydata_vert_interp
          51             :   public boundarydata_type
          52             :
          53             : contains
          54        3072 :   subroutine boundarydata_init(bndyfilename,phys_state,fieldnames,fieldcnt,bndydata,vertextrap)
          55             :     implicit none
          56             :     character(len=*),intent(in) :: bndyfilename
          57             :     type(physics_state), intent(in):: phys_state(begchunk:endchunk)
          58             :     integer,intent(in) :: fieldcnt
          59             :     character(len=*), intent(in) :: fieldnames(fieldcnt)
          60             :     type(boundarydata_type),intent(out) :: bndydata
          61             :     integer,intent(in), optional :: vertextrap ! if 0 set values outside output grid to 0
          62             :                                                ! if 1 set to boundary value
          63             :                                                ! if 2 set to cyclic boundaries
          64             :                                                ! if 3 leave on input data grid and extrapolate later
          65        3072 :     real(r8), pointer :: datain(:,:,:,:,:)
          66             :     integer :: lchnk
          67             :
          68        3072 :     bndydata%fieldcnt=fieldcnt
          69        3072 :     if(present(vertextrap)) then
          70           0 :        bndydata%vertextrap=vertextrap
          71             :     else
          72        3072 :        bndydata%vertextrap=0
          73             :     end if
          74             :     nullify(bndydata%fields)
          75             :
          76        3072 :     call boundarydata_read(phys_state,bndyfilename,fieldcnt,fieldnames,bndydata,datain)
          77             :
          78        3072 :     if(bndydata%iszonal) then
          79        3072 :        call boundarydata_interpolate(phys_state,datain,bndydata)
          80             :
          81             :        allocate(bndydata%datainst(size(bndydata%fields,1),size(bndydata%fields,2), &
          82       18432 :             begchunk:endchunk,bndydata%fieldcnt))
          83             :
          84        3072 :        deallocate(datain)
          85             :     end if
          86        3072 :   end subroutine boundarydata_init
          87             :
          88      741888 :   subroutine boundarydata_update(phys_state, bndydata, update_out)
          89             :     use interpolate_data,only : get_timeinterp_factors
          90             :     type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
          91             :     type(boundarydata_type), intent(inout) :: bndydata
          92             :     logical, intent(out), optional :: update_out
          93             :     real(r8) :: cdate
          94             :     integer :: nm, np, lchnk, j, k, fld, cols, cole, ncol, ndims
          95             :     real(r8) :: fact1, fact2
          96      741888 :     real(r8), allocatable :: datain(:,:,:,:,:)
          97             :     logical :: update
          98             :     integer :: latspan
          99             :     integer :: kmax
         100             :     integer :: count(4), start(4), ierr
         101             :
         102             :
         103      741888 :     call get_data_bounding_date_indices(bndydata%cdates,bndydata%nm,bndydata%np,cdate,update)
         104      741888 :     if(present(update_out)) update_out=update
         105      741888 :     nm= bndydata%nm
         106      741888 :     np= bndydata%np
         107             :
         108      741888 :     call get_timeinterp_factors(.true., np, bndydata%cdates(nm), bndydata%cdates(np), &
         109     1483776 :          cdate, fact1, fact2, _FILE)
         110             :
         111      741888 :     if(size(bndydata%fields,5).eq.2) then
         112      741888 :        nm=1
         113      741888 :        np=2
         114      741888 :        if(update) then ! we need to read in the next month and interpolate
         115           0 :           if(bndydata%isncol) then
         116           0 :              bndydata%fields(:,:,:,:,nm)=bndydata%fields(:,:,:,:,np)
         117           0 :              do lchnk=begchunk,endchunk
         118           0 :                 ncol=phys_state(lchnk)%ncol
         119           0 :                 cols=1
         120           0 :                 cole=cols+bndydata%count(cols,lchnk)-1
         121           0 :                 do while(cole<=ncol)
         122             :
         123           0 :                    if(bndydata%levsiz==1) then
         124           0 :                       ndims=2
         125           0 :                       start=(/bndydata%start(cols,lchnk),bndydata%np,-1,-1/)
         126           0 :                       count=(/bndydata%count(cols,lchnk),1,-1,-1/)
         127             :                    else
         128           0 :                       ndims=3
         129           0 :                       start=(/bndydata%start(cols,lchnk),bndydata%levsiz,bndydata%np,-1/)
         130           0 :                       count=(/bndydata%count(cols,lchnk),1,1,-1/)
         131             :                    end if
         132           0 :                    do fld=1,bndydata%fieldcnt
         133           0 :                       call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(fld) , &
         134           0 :                            bndydata%fields(cols:cole,:,lchnk,fld,np),  &
         135           0 :                            start(1:ndims), count(1:ndims)),&
         136           0 :                            _FILE,__LINE__)
         137             :
         138             :                    end do
         139           0 :                    if(cols==ncol) exit
         140           0 :                    cols=cols+bndydata%count(cols,lchnk)
         141           0 :                    cole=cols+bndydata%count(cols,lchnk)-1
         142             :                 end do
         143             :              end do
         144             :
         145             :           else
         146           0 :              allocate(datain(ptrlon,bndydata%levsiz,bndydata%latsiz,1,bndydata%fieldcnt))
         147           0 :              if(masterproc) then
         148           0 :                 count(1)=ptrlon
         149           0 :                 count(2)=bndydata%levsiz
         150           0 :                 count(3)=bndydata%latsiz
         151           0 :                 count(4)=1
         152           0 :                 start(1)=1
         153           0 :                 start(2)=1
         154           0 :                 start(3)=1
         155           0 :                 start(4)=bndydata%np
         156           0 :                 write(iulog,*) 'boundarydata reading data for month: ',bndydata%np
         157           0 :                 do fld=1,bndydata%fieldcnt
         158           0 :                    call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(fld), &
         159           0 :                         datain(:,:,:,:,fld), start, count),_FILE,__LINE__)
         160             :                 end do
         161             :              end if
         162             : #ifdef SPMD
         163           0 :              call mpibcast (datain, bndydata%levsiz*bndydata%latsiz*1*bndydata%fieldcnt, mpir8, 0, mpicom, ierr)
         164             : #endif
         165           0 :              bndydata%fields(:,:,:,:,nm) = bndydata%fields(:,:,:,:,np)
         166           0 :              call boundarydata_interpolate(phys_state,datain,bndydata)
         167           0 :              deallocate(datain)
         168             :           end if
         169             :        end if
         170             :     end if
         171      741888 :     kmax = size(bndydata%fields,2)
         172             :
         173     2596608 :     do fld=1,bndydata%fieldcnt
         174    10073448 :        do lchnk=begchunk,endchunk
         175     7476840 :           if(bndydata%isncol) then
         176           0 :              latspan = phys_state(lchnk)%ncol
         177             :           else
         178     7476840 :              latspan=phys_state(lchnk)%ulatcnt
         179             :           end if
         180   704677680 :           do k=1,kmax
         181 11409491205 :              do j=1,latspan
         182           0 :                 bndydata%datainst(j,k,lchnk,fld)=bndydata%fields(j,k,lchnk,fld,nm)*fact1 + &
         183 11402014365 :                      bndydata%fields(j,k,lchnk,fld,np)*fact2
         184             :              end do
         185             :           end do
         186             :        end do
         187             :     end do
         188      741888 :   end subroutine boundarydata_update
         189             :
         190             :
         191        3072 :   subroutine boundarydata_read(phys_state,bndyfilename,fieldcnt,fieldnames,bndydata,datain)
         192             :     !-----------------------------------------------------------------------
         193             :     !
         194             :     ! Purpose:
         195             :     ! Do initial read of time-variant boundary dataset, containing
         196             :     ! 12 monthly fields as a function of latitude and pressure.  Determine the two
         197             :     ! consecutive months between which the current date lies.
         198             :     !
         199             :     ! Method:
         200             :     !
         201             :     ! Author: NCAR CMS
         202             :     !-----------------------------------------------------------------------
         203             :     use ioFileMod, only : getfil
         204             :     use bnddyi_mod, only: bnddyi
         205             :
         206             :     implicit none
         207             :     type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
         208             :     character(len=*),intent(in) :: bndyfilename
         209             :     integer,intent(in) :: fieldcnt
         210             :     character(len=*), intent(in) :: fieldnames(fieldcnt)
         211             :     type(boundarydata_type), intent(inout) :: bndydata
         212             :     real(r8), pointer :: datain(:,:,:,:,:)  !
         213             :     !
         214             :     ! Local variables
         215             :     !
         216             :     integer :: londimid
         217             :     integer :: latdimid
         218             :     integer :: ncoldimid
         219             :     integer :: levdimid
         220             :     integer :: ilevdimid
         221             :     integer :: timdimid
         222             :     integer :: ndims
         223             :     integer :: dimlen
         224             :     integer :: ilevsiz
         225             :     integer :: ncolsiz
         226             :     character(len=nf90_max_name) :: dimname
         227             :
         228             :
         229             : !    integer ::  ncid              ! netcdf id for file
         230             :     integer ::  dateid            ! netcdf id for date variable
         231             :     integer ::  secid             ! netcdf id for seconds variable
         232             :     integer ::  lonid             ! netcdf id for longitude variable
         233             :     integer ::  ncolid             ! netcdf id for longitude variable
         234             :     integer ::  latid             ! netcdf id for latitude variable
         235             :     integer ::  levid             ! netcdf id for level variable
         236             :     integer ::  timid             ! netcdf id for time variable
         237             :     integer :: hybid
         238             :
         239             :     integer ::  dataid  ! netcdf id for data fields
         240             :
         241             :     integer ::  lonsiz            ! size of longitude dimension on tracer dataset
         242             :     integer ::  levsiz            ! size of level dimension on tracer dataset
         243             :     integer ::  latsiz            ! size of latitude dimension on tracer dataset
         244             :
         245             :     integer ::  j,n,k,nt,id          ! indices
         246             :     integer ::  ki,ko,ji,jo       ! indices
         247             :     integer :: date_tr(ptrtim), sec_tr(ptrtim)
         248             :
         249             :     integer :: dimids(4), start(4), count(4)
         250        3072 :     integer, pointer :: columnmap(:)
         251             :     real(r8) :: calday        ! current calendar day
         252        3072 :     real(r8), pointer :: pin(:)
         253        3072 :     real(r8), allocatable :: tmp_ps(:,:), tmp_fld(:,:,:)
         254             :     integer :: mincid,maxcid
         255             :     real(r8), allocatable, target:: lati(:)
         256             :     integer :: cols, cole
         257             :     integer ::  ierr, dimcnt
         258             :     integer ::  i, ncol, lchnk
         259             :     character(len=256) :: locfn    ! netcdf local filename to open
         260             :
         261             :     !
         262             :     !-----------------------------------------------------------------------
         263             :     !
         264             :     ! SPMD: Master reads dataset and does broadcast.  All subsequent interpolation is
         265             :     !       done in every process.  This is not required, one could remove this conditional
         266             :     !       and read the dataset independently on each task.
         267             :     !
         268        3072 :     if(masterproc) then
         269           4 :        write(iulog,*)'boundarydata_read: Reading from: ', trim(bndyfilename)
         270             : #ifndef USE_MASTERPROC
         271             :     end if
         272             : #endif
         273        3072 :     call getfil(bndyfilename, locfn)
         274             :     call handle_ncerr( nf90_open(locfn, 0, bndydata%ncid),&
         275        3072 :          _FILE,__LINE__)
         276             :
         277             :     !       write(iulog,*)'boundarydata_read: NCOPN returns id ',bndydata%ncid,' for file ',trim(locfn)
         278             :     !
         279             :     !------------------------------------------------------------------------
         280             :     ! Read tracer data
         281             :     !------------------------------------------------------------------------
         282             :     !
         283             :     ! Get dimension info
         284             :     !
         285        3072 :     nullify(columnmap)
         286        3072 :     nullify(pin)
         287             :     call handle_ncerr( nf90_inquire(bndydata%ncid, bndydata%ndims, unlimiteddimid=timdimid), &
         288        3072 :          _FILE,__LINE__)
         289        3072 :     ncolsiz=-1
         290        3072 :     levsiz=-1
         291        3072 :     lonsiz=-1
         292        3072 :     latsiz=-1
         293        3072 :     bndydata%isncol = .false.
         294       15360 :     do i=1,bndydata%ndims
         295             :        call handle_ncerr( nf90_inquire_dimension(bndydata%ncid, i, dimname, dimlen),&
         296       12288 :             _FILE,__LINE__)
         297       15360 :        if (dimname(1:3).eq.'lat') then
         298        3072 :           latdimid=i
         299        3072 :           latsiz=dimlen
         300        9216 :        else if (dimname(1:3) .eq.'lon') then
         301        3072 :           londimid=i
         302        3072 :           lonsiz=dimlen
         303        6144 :        else if (dimname(1:4) .eq. 'ncol') then
         304           0 :           ncoldimid=i
         305           0 :           ncolsiz=dimlen
         306           0 :           bndydata%isncol=.true.
         307        6144 :        else if (dimname(1:3) .eq. 'lev') then
         308        3072 :           levdimid=i
         309        3072 :           levsiz=dimlen
         310        3072 :        else if (dimname(1:4) .eq. 'ilev') then
         311       12288 :           ilevdimid=i
         312       12288 :           ilevsiz=dimlen
         313        3072 :        else if (dimname(1:4) .eq. 'time') then
         314        3072 :           if(timdimid/=i) then
         315           0 :              timdimid=i
         316             :           end if
         317        3072 :           bndydata%timsiz=dimlen
         318             :        else
         319           0 :           write(iulog,*) 'Warning: do not know how to handle dimension ',&
         320           0 :                trim(dimname), ' in boundarydata.F90:313'
         321             :        end if
         322             :     end do
         323             :
         324        3072 :     bndydata%iszonal = (latsiz>0 .and. lonsiz<=1)
         325        3072 :     if (bndydata%iszonal) then
         326        9216 :        allocate(bndydata%lat(latsiz))
         327             :     end if
         328        3072 :     if(bndydata%isncol) then
         329             : !       allocate (columnmap(ncolsiz))
         330             : !       call handle_ncerr( nf90_inq_varid(bndydata%ncid, 'ncol'    , ncolid),&
         331             : !            _FILE,__LINE__)
         332             : !       call handle_ncerr( nf90_get_var(bndydata%ncid,ncolid,columnmap), &
         333             : !            _FILE,__LINE__)
         334           0 :        if(levsiz>0) then
         335           0 :           allocate(bndydata%fields(pcols,levsiz,begchunk:endchunk,fieldcnt,2))
         336             :
         337           0 :           ierr = nf90_inq_varid(bndydata%ncid, 'PS', bndydata%psid)
         338           0 :           if(ierr.eq.NF90_NOERR) then
         339           0 :              allocate(bndydata%ps(pcols,begchunk:endchunk,2))
         340           0 :              allocate(bndydata%hybi(levsiz+1))
         341             :              call handle_ncerr(nf90_inq_varid(bndydata%ncid,'hybi',hybid),&
         342           0 :                   _FILE,__LINE__)
         343             :              call handle_ncerr( nf90_get_var(bndydata%ncid, hybid, bndydata%hybi ),&
         344           0 :                   _FILE,__LINE__)
         345             :           else
         346           0 :              call endrun('BOUNDARYDATA_READ: Did not recognize a vertical coordinate variable')
         347             :           end if
         348             :        else
         349           0 :           levsiz=1
         350           0 :           allocate(bndydata%fields(pcols,1,begchunk:endchunk,fieldcnt,2))
         351             :        end if
         352             :     else
         353       21504 :        allocate(datain(lonsiz,levsiz,latsiz,2,fieldcnt))
         354             :        !
         355             :        ! Check dimension info
         356             :        !
         357        3072 :        if (lonsiz/=ptrlon) then
         358           0 :           call endrun ('BOUNDARYDATA_READ: longitude dependence not implemented')
         359             :        endif
         360             :
         361        3072 :        if (bndydata%timsiz /= ptrtim) then
         362           0 :           write(iulog,*)'BOUNDARYDATA_READ: timsiz=',bndydata%timsiz,' must = ptrtim=',ptrtim
         363           0 :           call endrun
         364             :        end if
         365        3072 :        if( bndydata%vertextrap.lt.3) then
         366        9216 :           allocate(pin(levsiz))
         367             :        else
         368           0 :           allocate(bndydata%pin(levsiz))
         369           0 :           pin => bndydata%pin
         370             :        end if
         371        9216 :        allocate(bndydata%lat(latsiz))
         372             :
         373             :
         374       18432 :        allocate(datain(ptrlon,levsiz,latsiz,2,fieldcnt))
         375             :
         376             :        call handle_ncerr( nf90_inq_varid(bndydata%ncid, 'lat'    , latid),&
         377        3072 :             _FILE,__LINE__)
         378             :     end if
         379             :     !
         380             :     ! Determine necessary dimension and variable id's
         381             :     !
         382        9216 :     allocate(bndydata%cdates(bndydata%timsiz))
         383             :
         384             :     call handle_ncerr(nf90_inq_varid(bndydata%ncid, 'date'   , dateid), &
         385        3072 :          _FILE,__LINE__)
         386             :     call handle_ncerr( nf90_get_var(bndydata%ncid, dateid, date_tr),&
         387        3072 :          _FILE,__LINE__)
         388        3072 :     ierr = nf90_inq_varid(bndydata%ncid, 'datesec', secid)
         389        3072 :     if(ierr==NF90_NOERR) then
         390             :        call handle_ncerr( nf90_get_var(bndydata%ncid, secid , sec_tr),&
         391        3072 :             _FILE,__LINE__)
         392             :     else
         393           0 :        sec_tr=0
         394             :     end if
         395             :
         396        3072 :     if (mod(date_tr(1),10000)/100 /= 1) then
         397           0 :        call endrun ('(boundarydata_read): error when cycling data: 1st month must be 1')
         398             :     end if
         399        3072 :     if (mod(date_tr(ptrtim),10000)/100 /= 12) then
         400           0 :        call endrun ('(boundarydata_read): error when cycling data: last month must be 12')
         401             :     end if
         402             :     !
         403             :     !    return the calander dates of the file data
         404             :     !
         405       39936 :     do n=1,ptrtim
         406       39936 :        call bnddyi(date_tr(n), sec_tr(n), bndydata%cdates(n))
         407             :     end do
         408             : !    else
         409             : !       call handle_ncerr( nf90_inq_varid(bndydata%ncid, 'time', dateid),&
         410             : !            _FILE,__LINE__)
         411             : !
         412             : !       call handle_ncerr( nf90_get_var(bndydata%ncid, dateid, bndydata%cdates),&
         413             : !            _FILE,__LINE__)
         414             : !
         415             : !    end if
         416             : #ifdef USE_MASTERPROC
         417             :  else
         418             :     allocate(bndydata%cdates(ptrtim))
         419             :  end if
         420             : #ifdef SPMD
         421             :  call mpibcast (bndydata%cdates, ptrtim, mpir8, 0, mpicom, ierr)
         422             : #endif
         423             : #endif
         424        3072 :  bndydata%nm=12
         425        3072 :  bndydata%np=1
         426        3072 :  call get_data_bounding_date_indices(bndydata%cdates,bndydata%nm,bndydata%np)
         427             : #ifdef USE_MASTERPROC
         428             :  if(masterproc) then
         429             : #endif
         430             :     !
         431             :     ! Obtain entire date and sec variables. Assume that will always
         432             :     ! cycle over 12 month data.
         433             :     !
         434             :     !
         435             :     ! Obtain input data latitude and level arrays.
         436             :     !
         437        3072 :     if(bndydata%iszonal) then
         438             :        call handle_ncerr( nf90_get_var(bndydata%ncid, latid, bndydata%lat),&
         439        3072 :             _FILE,__LINE__)
         440        3072 :        ierr = nf90_inq_varid(bndydata%ncid, 'lev'    , levid)
         441             :        call handle_ncerr( nf90_get_var(bndydata%ncid, levid, pin ),&
         442        3072 :             _FILE,__LINE__)
         443             :     end if
         444             :
         445        9216 :     allocate(bndydata%dataid(fieldcnt))
         446        3072 :     if(masterproc) then
         447           4 :        write(iulog,*) 'boundarydata reading data for months: ',bndydata%nm,bndydata%np
         448             :     end if
         449       10752 :     do i=1,fieldcnt
         450        7680 :        call handle_ncerr( nf90_inq_varid(bndydata%ncid, fieldnames(i)   , bndydata%dataid(i)),&
         451       18432 :             _FILE,__LINE__)
         452             :     end do
         453        3072 :     if(bndydata%isncol) then
         454             :        allocate(bndydata%start(pcols,begchunk:endchunk), &
         455           0 :             bndydata%count(pcols,begchunk:endchunk))
         456             :
         457             : !
         458             : !  For i/o efficiency we read in a block of data which includes the data needed on this
         459             : !  processor but which may in fact include data not needed here.  physics cids are just the
         460             : !  offset into the file.
         461             : !
         462             :
         463           0 :        bndydata%start=-1
         464           0 :        bndydata%count=1
         465           0 :        mincid=2147483647
         466           0 :        maxcid=-1
         467           0 :        do lchnk=begchunk,endchunk
         468           0 :           ncol=phys_state(lchnk)%ncol
         469           0 :           i=minval(phys_state(lchnk)%cid(1:ncol))
         470           0 :           if(i < mincid) mincid = i
         471           0 :           i=maxval(phys_state(lchnk)%cid(1:ncol))
         472           0 :           if(i > maxcid) maxcid = i
         473             :        end do
         474             :
         475           0 :        allocate(tmp_ps(mincid:maxcid,2))
         476           0 :        start=(/mincid,bndydata%nm,1,-1/)
         477           0 :        if(bndydata%np>bndydata%nm) then
         478           0 :           count=(/maxcid-mincid+1,2,-1,-1/)
         479             :        else
         480           0 :           count=(/maxcid-mincid+1,1,-1,-1/)
         481             :        end if
         482           0 :        if(associated(bndydata%ps) ) then
         483             :           call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%psid , &
         484           0 :                tmp_ps(:,1:count(2)), start(1:2), &
         485             :                count(1:2)),&
         486           0 :                _FILE,__LINE__)
         487           0 :           if(bndydata%np<bndydata%nm) then
         488           0 :              start(2)=bndydata%np
         489             :              call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%psid , &
         490             :                   tmp_ps(:,2:2), start(1:2), &
         491             :                   count(1:2)),&
         492           0 :                   _FILE,__LINE__)
         493             :           end if
         494             :
         495           0 :           do lchnk=begchunk,endchunk
         496           0 :              do n=1,phys_state(lchnk)%ncol
         497           0 :                 bndydata%ps(n ,lchnk,:) = tmp_ps(phys_state(lchnk)%cid(n),:)
         498             :              end do
         499             :           end do
         500           0 :           deallocate(tmp_ps)
         501             :
         502             :        end if
         503             :
         504           0 :        if(levsiz>1) then
         505             :           dimcnt=3
         506             :        else
         507           0 :           dimcnt=2
         508             :        end if
         509           0 :        start(2)=1
         510           0 :        count(2)=levsiz
         511             :
         512           0 :        if(bndydata%np>bndydata%nm) then
         513           0 :           count(dimcnt)=2
         514             :        else
         515           0 :           count(dimcnt)=1
         516             :        end if
         517           0 :        start(dimcnt)=bndydata%nm
         518             :
         519           0 :        allocate(tmp_fld(mincid:maxcid,count(2),2))
         520             :
         521           0 :        do i=1,fieldcnt
         522           0 :           call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(i) , &
         523           0 :                tmp_fld(:,:,1:count(dimcnt)), &
         524             :                start(1:dimcnt), count(1:dimcnt)),&
         525           0 :                _FILE,__LINE__)
         526             :
         527           0 :           do lchnk=begchunk,endchunk
         528           0 :              do n=1,phys_state(lchnk)%ncol
         529           0 :                 bndydata%fields(n,:,lchnk,i,1:count(dimcnt)) = tmp_fld(phys_state(lchnk)%cid(n),:,:)
         530             :              end do
         531             :           end do
         532             :        end do
         533           0 :        if(bndydata%np<bndydata%nm) then
         534           0 :           start(dimcnt)=bndydata%np
         535           0 :           do i=1,fieldcnt
         536           0 :              call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(i) , &
         537             :                   tmp_fld(:,:,2:2), &
         538             :                   start(1:dimcnt), count(1:dimcnt)),&
         539           0 :                   _FILE,__LINE__)
         540           0 :              do lchnk=begchunk,endchunk
         541           0 :                 do n=1,phys_state(lchnk)%ncol
         542           0 :                    bndydata%fields(n,:,lchnk,i,1:count(dimcnt)) = tmp_fld(phys_state(lchnk)%cid(n),:,:)
         543             :                 end do
         544             :              end do
         545             :           end do
         546             :        end if
         547           0 :        deallocate(tmp_fld)
         548             : !       deallocate(columnmap)
         549             :     else
         550             :        !
         551             :        ! get the dimension orientation info from the first variable assume but verify that
         552             :        ! all variables requested have the same orientation
         553             :        !
         554        3072 :        allocate(bndydata%start(4,1),bndydata%count(4,1))
         555           0 :        call handle_ncerr( nf90_inquire_variable(bndydata%ncid,bndydata%dataid(1), &
         556        3072 :             ndims=bndydata%ndims,dimids=bndydata%dimids),_FILE,__LINE__)
         557             :
         558       18432 :        bndydata%start=1
         559       15360 :        do id=1,bndydata%ndims
         560       15360 :           if(bndydata%dimids(id).eq.londimid) then
         561        3072 :              bndydata%map(id)=1
         562        3072 :              bndydata%count(id,1)=lonsiz
         563        9216 :           else if(bndydata%dimids(id).eq.levdimid) then
         564        3072 :              bndydata%map(id)=lonsiz
         565        3072 :              bndydata%count(id,1)=levsiz
         566        6144 :           else if(bndydata%dimids(id).eq.latdimid) then
         567        3072 :              bndydata%map(id)=lonsiz
         568        6144 :              if(any(bndydata%dimids.eq.levdimid)) bndydata%map(id)=bndydata%map(id)*levsiz
         569        3072 :              bndydata%count(id,1)=latsiz
         570        3072 :           else if(bndydata%dimids(id).eq.timdimid) then
         571        3072 :              bndydata%map(id)=lonsiz*latsiz
         572        6144 :              if(any(bndydata%dimids.eq.levdimid)) bndydata%map(id)=bndydata%map(id)*levsiz
         573        3072 :              bndydata%count(id,1)=2
         574        3072 :              bndydata%start(id,1)=bndydata%nm
         575        3072 :              bndydata%thistimedim=id
         576             :           else
         577           0 :              write(iulog,*) __LINE__,fieldnames(i),id,bndydata%dimids(id),londimid, &
         578           0 :                   levdimid,latdimid,timdimid
         579           0 :              call endrun(_FILE)
         580             :           end if
         581             :        end do
         582             :
         583       10752 :        do i=1,fieldcnt
         584        7680 :           call handle_ncerr( nf90_inq_varid(bndydata%ncid, fieldnames(i), bndydata%dataid(i)),&
         585       15360 :                _FILE,__LINE__)
         586             :
         587           0 :           call handle_ncerr( nf90_inquire_variable(bndydata%ncid,bndydata%dataid(i), &
         588        7680 :                ndims=ndims,dimids=dimids),_FILE,__LINE__)
         589             :           if(ndims/=bndydata%ndims .or. dimids(1)/=bndydata%dimids(1).or.&
         590        7680 :                dimids(2)/=bndydata%dimids(2) .or. dimids(3)/=bndydata%dimids(3)) then
         591           0 :              call endrun('BOUNDARYDATA_READ: Variable dims or order does not match')
         592             :           end if
         593             :
         594       10752 :           if(bndydata%np .gt. bndydata%nm) then
         595           0 :              call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(i), &
         596             :                   datain(:,:,:,:,i),bndydata%start(:,1), bndydata%count(:,1), &
         597           0 :                   map=bndydata%map),_FILE,__LINE__)
         598             :           else
         599        7680 :              bndydata%count(bndydata%thistimedim,1)=1
         600        7680 :              bndydata%start(bndydata%thistimedim,1)=bndydata%nm
         601           0 :              call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(i), &
         602           0 :                   datain(:,:,:,1:1,i), bndydata%start(:,1), bndydata%count(:,1), &
         603        7680 :                   map=bndydata%map), _FILE,__LINE__)
         604             :
         605        7680 :              bndydata%start(bndydata%thistimedim,1)=bndydata%np
         606           0 :              call handle_ncerr( nf90_get_var(bndydata%ncid, bndydata%dataid(i), &
         607           0 :                   datain(:,:,:,2:2,i), bndydata%start(:,1), bndydata%count(:,1), &
         608        7680 :                   map=bndydata%map), _FILE,__LINE__)
         609             :
         610             :           endif
         611             :        end do
         612             :
         613             :     end if
         614             : #ifdef USE_MASTERPROC
         615             :  end if
         616             : #ifdef SPMD
         617             :  call mpibcast (levsiz, 1, mpiint, 0, mpicom, ierr)
         618             :  call mpibcast (latsiz, 1, mpiint, 0, mpicom, ierr)
         619             : #endif
         620             : #endif
         621        3072 :  bndydata%latsiz=latsiz
         622        3072 :  bndydata%levsiz=levsiz
         623             : #ifdef USE_MASTERPROC
         624             :  if(.not. masterproc) then
         625             :     if( bndydata%vertextrap.lt.3) then
         626             :        allocate(pin(levsiz))
         627             :     else
         628             :        allocate(bndydata%pin(levsiz))
         629             :        pin => bndydata%pin
         630             :     end if
         631             :     allocate(bndydata%lat(latsiz))
         632             :     allocate(datain(ptrlon,levsiz,latsiz,2,fieldcnt))
         633             :  endif
         634             : #ifdef SPMD
         635             :  call mpibcast (bndydata%lat, latsiz, mpir8, 0, mpicom, ierr)
         636             :  call mpibcast (pin, levsiz, mpir8, 0, mpicom, ierr)
         637             :  call mpibcast (datain, levsiz*latsiz*2*fieldcnt, mpir8, 0, mpicom, ierr)
         638             :
         639             : #endif
         640             : #endif
         641             :  ! Convert input pressure from millibars to pascals.
         642        3072 :  if(associated(pin)) then
         643      175104 :     pin=pin*100._r8
         644        3072 :     if(bndydata%vertextrap.lt.3) then
         645        9216 :        allocate(bndydata%zi(levsiz))
         646             :     !
         647             :     !
         648             :     ! Convert input pressure levels to height (m).
         649             :     !
         650      175104 :        do k=1,levsiz
         651      175104 :           bndydata%zi(k) = 7.0e3_r8 * log (1.0e5_r8 / pin(k))
         652             :        end do
         653        3072 :        deallocate(pin)
         654             :     end if
         655             :  end if
         656        6144 : end subroutine boundarydata_read
         657             :
         658        3072 :   subroutine boundarydata_interpolate(phys_state, datain, bndydata)
         659             :     use ref_pres, only : pref_mid
         660             :     use interpolate_data,only : interp_type, lininterp_init, &
         661             :          lininterp_finish, lininterp
         662             :     use physconst, only: pi
         663             :
         664             :     type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
         665             :     real(r8),intent(in)                  :: datain(:,:,:,:,:)
         666             :     type(boundarydata_type), intent(inout) :: bndydata
         667             :     type(interp_type) :: interp_wgts, lev_wgts
         668             :
         669             :     integer :: k, lchnk, nt, j, fcnt
         670             :     real(r8) :: zo(pver)
         671             :     real(r8) :: lato(pcols)
         672             :     integer :: ulatcnt
         673             :     integer :: maxlatcnt
         674             :     integer :: timesize, tvalout
         675             :
         676             :     !------------------------------------------------------------------------
         677             :     ! Interpolate tracer data to model grid
         678             :     !------------------------------------------------------------------------
         679             :     !
         680             :     ! Loop over all input times.
         681             :     !
         682             :
         683        3072 :     timesize=2
         684             :
         685        3072 :     maxlatcnt=1
         686       15456 :     do lchnk=begchunk,endchunk
         687       15456 :        maxlatcnt=max(maxlatcnt,phys_state(lchnk)%ulatcnt)
         688             :     end do
         689        3072 :     if(bndydata%vertextrap.lt.3) then
         690             :        !
         691             :        ! Convert approximate cam pressure levels to height (m).
         692             :        !
         693      288768 :        do k=1,pver
         694      288768 :           zo (k) = 7.0e3_r8 * log (1.0e5_r8 / pref_mid(k))
         695             :        end do
         696             :
         697        3072 :        call lininterp_init(bndydata%zi,size(bndydata%zi),zo,pver,bndydata%vertextrap,lev_wgts)
         698        3072 :        if(.not. associated(bndydata%fields)) then
         699       18432 :           allocate(bndydata%fields(maxlatcnt,pver,begchunk:endchunk,bndydata%fieldcnt,timesize))
         700    97890876 :           bndydata%fields=0_r8
         701             :        end if
         702             :     else
         703           0 :        if(.not. associated(bndydata%fields)) then
         704           0 :           allocate(bndydata%fields(maxlatcnt,bndydata%levsiz,begchunk:endchunk,bndydata%fieldcnt,timesize))
         705           0 :           bndydata%fields=0_r8
         706             :        end if
         707             :     endif
         708       15456 :     do lchnk=begchunk,endchunk
         709       12384 :        ulatcnt=phys_state(lchnk)%ulatcnt
         710             :
         711             :        !
         712             :        ! Convert cam model latitudes to degrees.
         713             :        ! Input model latitudes already in degrees.
         714             :        !
         715      203068 :        do j=1,ulatcnt
         716      203068 :           lato(j) = phys_state(lchnk)%ulat(j)*180._r8/pi
         717             :        end do
         718             :
         719       12384 :        call lininterp_init(bndydata%lat,size(bndydata%lat),lato(1:ulatcnt),ulatcnt,1,interp_wgts)
         720       12384 :        timesize =  size(datain,4)
         721       43344 :        do fcnt=1,bndydata%fieldcnt
         722      105264 :           do nt = 1,timesize
         723       61920 :              if(timesize.gt.1) then
         724       61920 :                 tvalout=nt
         725             :              else
         726             :                 tvalout=2
         727             :              end if
         728       92880 :             if(bndydata%vertextrap.lt.3) then
         729             :                 call lininterp(transpose(datain(1,:,:,nt,fcnt)),bndydata%latsiz,bndydata%levsiz, &
         730    26916600 :                      bndydata%fields(1:ulatcnt,:,lchnk,fcnt,tvalout), ulatcnt, pver, interp_wgts, lev_wgts)
         731             :              else
         732           0 :                 do k=1,bndydata%levsiz
         733           0 :                    call lininterp(datain(1,k,:,nt,fcnt),bndydata%latsiz, &
         734           0 :                         bndydata%fields(1:ulatcnt,k,lchnk,fcnt,tvalout), ulatcnt, interp_wgts)
         735             :                 end do
         736             :              end if
         737             :           end do
         738             :        end do                 ! end loop over time samples
         739       15456 :        call lininterp_finish(interp_wgts)
         740             :     end do
         741        3072 :     if(bndydata%vertextrap.lt.3) &
         742        3072 :          call lininterp_finish(lev_wgts)
         743             :
         744        3072 :     return
         745        3072 :   end subroutine boundarydata_interpolate
         746             :
         747      744960 :   subroutine get_data_bounding_date_indices(cdates,nm,np, cdayout, update)
         748        3072 :     use time_manager, only: get_curr_date, get_perp_date, get_curr_calday, &
         749             :          is_perpetual
         750             :     real(r8), intent(in) :: cdates(ptrtim)
         751             :     real(r8), intent(out),optional :: cdayout
         752             :     logical, intent(out) ,optional :: update
         753             :     integer, intent(inout) :: nm, np
         754             :     integer :: n, np1
         755             :     real(r8) :: calday
         756             :     integer :: yr, mon, day   ! components of a date
         757             :     integer :: ncdate         ! current date in integer format [yyyymmdd]
         758             :     integer :: ncsec          ! current time of day [seconds]
         759             :
         760      744960 :     calday = get_curr_calday()
         761      744960 :     if(present(cdayout)) cdayout=calday
         762      744960 :     if(present(update)) update=.false. ! initialize output variable
         763             :
         764             :     if(min(nm,np) .ge. 1 .and. max(nm,np) .le. 12 .and.  &
         765      744960 :          calday>cdates(nm) .and. calday<=cdates(np)) return
         766      744960 :     if((nm==12 .and. np==1) .and. (calday <= cdates(np) .or. &
         767             :          calday > cdates(nm))) return
         768             :
         769           0 :     if(present(update)) update=.true.
         770             :
         771           0 :     if(calday <= cdates(1) .or. calday > cdates(12)) then
         772           0 :        nm=12
         773           0 :        np=1
         774             :     else
         775           0 :        nm=0
         776           0 :        do n=1,ptrtim-1
         777           0 :           if(calday>cdates(n) .and. calday<=cdates(n+1)) then
         778           0 :              nm=n
         779           0 :              np=n+1
         780             :           end if
         781             :        end do
         782           0 :        if(nm .eq. 0) then
         783           0 :           if ( is_perpetual() ) then
         784           0 :              call get_perp_date(yr, mon, day, ncsec)
         785             :           else
         786           0 :              call get_curr_date(yr, mon, day, ncsec)
         787             :           end if
         788           0 :           ncdate = yr*10000 + mon*100 + day
         789             :
         790           0 :           write(iulog,*)'model date:', ncdate, ncsec,'boundary data dates:', cdates
         791           0 :           call endrun('get_data_bounding_date_indices: Failed to find dates bracketing dates')
         792             :        end if
         793             :     end if
         794             :
         795      744960 :   end subroutine get_data_bounding_date_indices
         796             :
         797             :
         798             :   !================================================================================================
         799           0 :   subroutine boundarydata_vert_interp(lchnk, ncol, levsiz, fldcnt, pin, pmid, datain, dataout)
         800             :     !-----------------------------------------------------------------------
         801             :     !
         802             :     ! Purpose: Interpolate ozone from current time-interpolated values to model levels
         803             :     !
         804             :     ! Method: Use pressure values to determine interpolation levels
         805             :     !
         806             :     ! Author: Bruce Briegleb
         807             :     !
         808             :     !--------------------------------------------------------------------------
         809             :     implicit none
         810             :     ! Arguments
         811             :     !
         812             :     integer, intent(in) :: lchnk               ! chunk identifier
         813             :     integer, intent(in) :: ncol                ! number of atmospheric columns
         814             :     integer, intent(in) :: levsiz
         815             :     integer, intent(in) :: fldcnt
         816             :     real(r8), intent(in) :: pin(levsiz)
         817             :     real(r8), intent(in) :: pmid(pcols,pver)   ! level pressures (mks)
         818             :     real(r8), intent(in) :: datain(pcols,levsiz,fldcnt)
         819             :     real(r8), intent(out) :: dataout(pcols,pver,fldcnt)    ! ozone mass mixing ratio
         820             :     !
         821             :     ! local storage
         822             :     !
         823             :
         824             :     integer ::  i                   ! longitude index
         825             :     integer ::  k, kk, kkstart      ! level indices
         826             :     integer ::  kupper(pcols)       ! Level indices for interpolation
         827             :     integer ::  kount               ! Counter
         828             :     integer ::  fld
         829             :     real(r8) dpu                ! upper level pressure difference
         830             :     real(r8) dpl                ! lower level pressure difference
         831             :     !--------------------------------------------------------------------------
         832             :     !
         833             :     ! Initialize index array
         834             :     !
         835           0 :     do i=1,ncol
         836           0 :        kupper(i) = 1
         837             :     end do
         838             :
         839           0 :     do k=1,pver
         840             :        !
         841             :        ! Top level we need to start looking is the top level for the previous k
         842             :        ! for all longitude points
         843             :        !
         844             :        kkstart = levsiz
         845           0 :        do i=1,ncol
         846           0 :           kkstart = min0(kkstart,kupper(i))
         847             :        end do
         848             :        kount = 0
         849             :        !
         850             :        ! Store level indices for interpolation
         851             :        !
         852           0 :        do kk=kkstart,levsiz-1
         853           0 :           do i=1,ncol
         854           0 :              if (pin(kk).lt.pmid(i,k) .and. pmid(i,k).le.pin(kk+1)) then
         855           0 :                 kupper(i) = kk
         856           0 :                 kount = kount + 1
         857             :              end if
         858             :           end do
         859             :           !
         860             :           ! If all indices for this level have been found, do the interpolation and
         861             :           ! go to the next level
         862             :           !
         863           0 :           if (kount.eq.ncol) then
         864           0 :              do fld=1,fldcnt
         865           0 :                 do i=1,ncol
         866           0 :                    dpu = pmid(i,k) - pin(kupper(i))
         867           0 :                    dpl = pin(kupper(i)+1) - pmid(i,k)
         868           0 :                    dataout(i,k,fld) = (datain(i,kupper(i),fld )*dpl + &
         869           0 :                         datain(i,kupper(i)+1,fld)*dpu)/(dpl + dpu)
         870             :
         871             :                 end do
         872             :              end do
         873             :              goto 35
         874             :           end if
         875             :        end do
         876             :        !
         877             :        ! If we've fallen through the kk=1,levsiz-1 loop, we cannot interpolate and
         878             :        ! must extrapolate from the bottom or top data level for at least some
         879             :        ! of the longitude points.
         880             :        !
         881           0 :        do fld=1,fldcnt
         882           0 :           do i=1,ncol
         883           0 :              if (pmid(i,k) .lt. pin(1)) then
         884           0 :                 dataout(i,k,fld) = datain(i,1,fld)*pmid(i,k)/pin(1)
         885           0 :              else if (pmid(i,k) .gt. pin(levsiz)) then
         886           0 :                 dataout(i,k,fld) = datain(i,levsiz,fld)
         887             :              else
         888           0 :                 dpu = pmid(i,k) - pin(kupper(i))
         889           0 :                 dpl = pin(kupper(i)+1) - pmid(i,k)
         890           0 :                 dataout(i,k,fld) = (datain(i,kupper(i),fld )*dpl + &
         891           0 :                      datain(i,kupper(i)+1,fld)*dpu)/(dpl + dpu)
         892             :              end if
         893             :           end do
         894             :        end do
         895           0 :        if (kount.gt.ncol) then
         896           0 :           call endrun ('ozone_data_vert_interp: Bad ozone data: non-monotonicity suspected')
         897             :        end if
         898           0 : 35     continue
         899             :     end do
         900      744960 :   end subroutine boundarydata_vert_interp
         901             : #if 0
         902             :   subroutine ncol_read_bracket(cid,columnmap,start,count,ncol)
         903             :     integer, intent(in) :: cid(:), columnmap(:), ncol
         904             :     integer, intent(out) :: start(:), count(:)
         905             :
         906             :     integer :: i, j, tcol
         907             :
         908             :     tcol = size(columnmap)
         909             :     count=1
         910             :     do i=1,ncol
         911             : #if 1
         912             :         start(i)=cid(i)
         913             :         count(i)=1
         914             : #else
         915             :        do j=1,tcol
         916             :           if(columnmap(j).eq.cid(i)) then
         917             :              start(i)=j
         918             :              count(i)=1
         919             :              exit
         920             :           end if
         921             :        end do
         922             : #endif
         923             :     end do
         924             :     do i=1,ncol-1
         925             :        do j=1,ncol-i
         926             :           if(columnmap(start(i+j)).eq.columnmap(start(i)+j)) then
         927             :              count(i)=count(i)+1
         928             :           else
         929             :              exit
         930             :           end if
         931             :        end do
         932             :     end do
         933             :     write(iulog,*) __LINE__,cid(1),cid(ncol),minval(cid(1:ncol)),maxval(cid(1:ncol))
         934             :
         935             :   end subroutine ncol_read_bracket
         936             : #endif
         937           0 : end module boundarydata

Generated by: LCOV version 1.14