Skip to content
Snippets Groups Projects
Select Git revision
  • c9264d9d9ea0ef5d17591be612781273314faae6
  • master default protected
  • est-autem-a-officia-quibusdam-et-dolor
  • dignissimos-libero-alias-distinctio-sequi-mollitia-quia
  • repudiandae-quia-repellat-ipsa-enim-pariatur-quae
  • voluptatibus-ut-earum-fuga-reprehenderit-repudiandae-id
  • doloremque-qui-facere-quo-ea-vel-nostrum
  • quod-expedita-vitae-voluptatum-quo-qui-ipsum
  • aliquam-ratione-assumenda-quos-architecto-tempora-pariatur
  • cupiditate-hic-molestias-facilis-non-qui-praesentium
  • architecto-consequuntur-cupiditate-quo-delectus-similique-sunt
  • eaque-voluptatibus-omnis-labore-aut-qui-possimus
  • dicta-veniam-adipisci-rem-consequatur-ut-delectus
  • beatae-nulla-eum-aliquid-ut-nesciunt-commodi
  • itaque-deserunt-et-quos-non-sit-ut
  • debitis-repellat-tempora-accusantium-quia-ad-nam
  • dignissimos-modi-autem-dolores-fugiat-ipsum-officiis
  • unde-et-enim-aut-aut-dignissimos-atque
  • dolorem-quisquam-mollitia-quia-cum-quam-dolores
  • iure-rem-veritatis-ullam-voluptas-error-ad
  • iste-temporibus-adipisci-error-exercitationem-eaque-omnis
21 results

ADV_wind_stacked.f90

Blame
  • Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    ADV_wind_stacked.f90 33.06 KiB
    !*****************************************************************
    !
    ! MODULE NAME:
    !    ADV_wind
    ! FUNCTION:
    !    calculate the windfield for the advection problem
    ! CONTAINS:
    !-----------------------------------------------------------------
    !
    ! NAME:
    !    slm_windfield
    ! FUNCTION:
    !    calculate the advecting force for simple advection
    ! SYNTAX:
    !    real.arr= slm_windfield(real.arr, real)
    ! ON INPUT:
    !    r_coord: coordinates of point        real
    !    r_time:  time coordinate (optional)    real
    ! ON OUTPUT:
    !    r_field: windfield            real
    ! CALLS:
    !
    ! COMMENTS:
    !
    !-----------------------------------------------------------------
    !
    ! PUBLIC:
    !
    ! COMMENTS:
    !
    ! USES:
    !    MISC_globalparam, MISC_error
    ! LIBRARIES:
    !
    ! REFERENCES:
    !
    ! VERSION(S):
    !    1. original version        j. behrens    12/97
    !    2. bug fix concerning interp.    j. behrens    2/98
    !    3. compliant to amatos 1.0    j. behrens    12/2000
    !    4. compliant to amatos 1.2    j. behrens    3/2002
    !    5. new version for course tracer transport    j. behrens 01/2012
    !
    !*****************************************************************
        MODULE ADV_wind
          USE GRID_api
          USE FLASH_parameters
          PRIVATE
          INTEGER, PARAMETER        :: i_ioerr=0
          REAL (KIND = GRID_SR)     :: r_intervallen, &
                                       r_scalfacx, r_scalfacy, r_scalfacz, &
                                       r_readlast, r_scalinvx, r_scalinvy, r_scalinvz
    !       REAL (KIND = GRID_SR), DIMENSION(:,:,:,:), ALLOCATABLE :: r_flowx
    !       REAL (KIND = GRID_SR), DIMENSION(:,:,:,:), ALLOCATABLE :: r_flowy
    !       REAL (KIND = GRID_SR), DIMENSION(:,:,:,:), ALLOCATABLE :: r_flowz
    !       REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE       :: r_lat
    !       REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE       :: r_lon
    !       REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE       :: r_alt
    !       REAL (KIND = GRID_SR), DIMENSION(:), ALLOCATABLE       :: r_times
          REAL*4, DIMENSION(:,:,:,:), ALLOCATABLE :: r_flowx
          REAL*4, DIMENSION(:,:,:,:), ALLOCATABLE :: r_flowy
          REAL*4, DIMENSION(:,:,:,:), ALLOCATABLE :: r_flowz
          REAL*4, DIMENSION(:), ALLOCATABLE       :: r_lat
          REAL*4, DIMENSION(:), ALLOCATABLE       :: r_lon
          REAL*4, DIMENSION(:), ALLOCATABLE       :: r_alt
          REAL*8, DIMENSION(:), ALLOCATABLE       :: r_times
          INTEGER (KIND = GRID_SI)  :: i_lon, i_lat, i_layer, i_timesteps
          INTEGER (KIND = GRID_SI)  :: i_pointsx, i_pointsy, i_pointsz, &
                                       i_timeinterval, i_intervalnum
    
          PUBLIC                    :: slm_windfield, slm_windinit, slm_windquit
          PUBLIC                    :: slm_vertwind
    
          CONTAINS
    !*****************************************************************
          FUNCTION slm_windfield(r_coord, i_layerid, r_time) RESULT (r_field)
    
    !---------- local declarations
    
          IMPLICIT NONE
    
          REAL (KIND = GRID_SR), DIMENSION(GRID_dimension), INTENT(in) :: r_coord
          INTEGER (KIND = GRID_SI), INTENT(in)                         :: i_layerid
          REAL (KIND = GRID_SR), INTENT(in), OPTIONAL                  :: r_time
          REAL (KIND = GRID_SR), DIMENSION(GRID_dimension)             :: r_field
    
          REAL (KIND = GRID_SR)                                        :: r_tim
    
    !---------- set time
    
          IF(present(r_time)) THEN
            r_tim= r_time
          ELSE
            r_tim= 0.0_GRID_SR
          END IF
    
    !---------- decide, if data has to be read
    
          data_read: IF(r_readlast <= r_tim) THEN
    
    !---------- update values for next open
    
            r_readlast    = r_readlast+ r_intervallen
            i_timeinterval= i_timeinterval+ 1
          END IF data_read
    
    !---------- interpolate to coordinate
    
          r_field= data_interpol(r_coord, i_layerid, i_timeinterval)
    
          RETURN
          END FUNCTION slm_windfield
    
    !*****************************************************************
          FUNCTION slm_vertwind(r_coord, i_layerid) RESULT (r_wwind)
    
    !---------- declarations
    
          IMPLICIT NONE
    
    !---------- i/o parameters
    
          REAL (KIND = GRID_SR), DIMENSION(GRID_dimension), INTENT(in) :: r_coord
          INTEGER (KIND = GRID_SI), INTENT(in)                         :: i_layerid
          REAL (KIND = GRID_SR)                                        :: r_wwind
    
    !---------- interpolate to coordinate
    ! Note that the time interval will be globally adjusted in the 
    ! horizontal wind field computation!
    
          r_wwind= vertical_interpol(r_coord, i_layerid, i_timeinterval)
    
          RETURN
          END FUNCTION slm_vertwind
    !*****************************************************************
          SUBROUTINE slm_windinit(p_control)
    
    !---------- local declarations
    
          IMPLICIT NONE
    
          TYPE (control_struct)     :: p_control
          INTEGER (KIND = GRID_SI)  :: i_iofil= 19
          CHARACTER (len=80)        :: a_filrow
          INTEGER (KIND = GRID_SI)  :: i_iost, i_ioend, i_alct, i_cnt
          INTEGER (KIND = GRID_SI)  :: i_mon, i_day, i_year
          REAL (KIND = GRID_SR)     :: r_tim, r_start, r_t
    
    
    !---------- read current data from NetCDF file
    
          CALL read_netcdf_currents(p_control)
    
    !---------- initialize some values
    
          r_readlast    = 0.0_GRID_SR
          i_timeinterval= 0_GRID_SI
          r_intervallen = 21600._GRID_SR   ! we get a new wind value every 6 hours
          r_start       = 1937256._GRID_SR ! this is the hours from 1800-01-01 for 2021...
    
    !---------- compute starting time frame
    
          i_year= p_control%tst%tst_int(1,2)
          i_mon = p_control%tst%tst_int(1,3)
          i_day = p_control%tst%tst_int(1,4)
          r_tim = time_from_datecomponents(i_year,i_mon,i_day,0,0,0,r_start)
          time_loop: DO i_cnt=1,i_timesteps-1
            r_t = r_times(i_cnt)
            IF(r_times(i_cnt) .LT. r_tim) THEN
              IF(r_times(i_cnt+1) .GE. r_tim) THEN
                i_timeinterval = i_cnt
                exit time_loop
              END IF
            END IF
          END DO time_loop 
    
    !      RETURN
          END SUBROUTINE slm_windinit
    
    !*****************************************************************
          SUBROUTINE slm_windquit
    
    !---------- local declarations
    
          IMPLICIT NONE
    
    !---------- deallocate wind data arrays
    
          DEALLOCATE(r_flowx, r_flowy, r_flowz)
          DEALLOCATE(r_lon, r_lat, r_alt, r_times)
    
    !      RETURN
          END SUBROUTINE slm_windquit
    
    !*****************************************************************
          SUBROUTINE read_netcdf_currents(p_control)
    
    !---------- local declarations
    
          IMPLICIT NONE
          INCLUDE "netcdf.inc"
          
    !---------- input parameters      
          TYPE (control_struct), INTENT(in)      :: p_control
          
    !---------- local variables
          CHARACTER (len=io_fillen)              :: c_uname, c_vname, c_wname
          INTEGER (KIND = GRID_SI)               :: i_alct, i_ncstat
          INTEGER (KIND = GRID_SI)               :: i_filexid, i_fileyid, i_filezid
          INTEGER (KIND = GRID_SI)               :: i_dimid, i_varid
          INTEGER (KIND = GRID_SI)               :: i_tm, i_ln, i_lt, i_ly
          REAL (KIND = GRID_SR), DIMENSION(:,:,:,:), ALLOCATABLE :: r_auxx
          REAL (KIND = GRID_SR), DIMENSION(:,:,:,:), ALLOCATABLE :: r_auxy
          REAL (KIND = GRID_SR), DIMENSION(:,:,:,:), ALLOCATABLE :: r_auxz
    
    !---------- initialize file names
          WRITE(GRID_parameters%ioout,*) &
            '***** ***** ***** ***** INITIALIZE  ***** ***** ***** *****'
          WRITE(GRID_parameters%ioout,*) &
            '***** reading NetCDF file header ...                  *****'
    
          WRITE(c_uname,*) trim(p_control%tst%tst_char(1))
          c_uname= adjustl(c_uname)
          WRITE(c_vname,*) trim(p_control%tst%tst_char(2))
          c_vname= adjustl(c_vname)
          WRITE(c_wname,*) trim(p_control%tst%tst_char(3))
          c_wname= adjustl(c_wname)
    
    !---------- open current files
    
          i_ncstat= nf_open(c_uname,NF_NOWRITE,i_filexid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not open data file for u')
          i_ncstat= nf_open(c_vname,NF_NOWRITE,i_fileyid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not open data file for v')
          i_ncstat= nf_open(c_wname,NF_NOWRITE,i_filezid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not open data file for w')
        
    !---------- determine lon/lat/time/levels dimension sizes (grid size of currents field)
    
          i_ncstat= nf_inq_dimid(i_filexid, 'lon', i_dimid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not identify lon dimension')
          i_ncstat= nf_inq_dimlen(i_filexid, i_dimid, i_lon)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read lon dimension')
          
          i_ncstat= nf_inq_dimid(i_filexid, 'lat', i_dimid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not identify lat dimension')
          i_ncstat= nf_inq_dimlen(i_filexid, i_dimid, i_lat)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read lat dimension')
    
          i_ncstat= nf_inq_dimid(i_filexid, 'time', i_dimid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not identify time dimension')
          i_ncstat= nf_inq_dimlen(i_filexid, i_dimid, i_timesteps)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read time dimension')
    
          i_ncstat= nf_inq_dimid(i_filexid, 'level', i_dimid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not identify time dimension')
          i_ncstat= nf_inq_dimlen(i_filexid, i_dimid, i_layer)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read time dimension')
    
    !---------- make sure layer number in file is larger than required from parameters
    
          IF(i_layer .LT. p_control%tst%tst_int(1,1)) &
            CALL grid_error(c_error='[read_netcdf_currents]: too few layers in NetCDF file')
    
    !---------- check dimension match in v- and w-component files
    
          i_ncstat= nf_inq_dimid(i_fileyid, 'lon', i_dimid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not identify lon dimension')
          i_ncstat= nf_inq_dimlen(i_fileyid, i_dimid, i_ln)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read lon dimension')
          IF(i_lon /= i_ln) &
            CALL grid_error(c_error='[read_netcdf_currents]: lon dimension mismatch')
    
          i_ncstat= nf_inq_dimid(i_fileyid, 'lat', i_dimid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not identify lat dimension')
          i_ncstat= nf_inq_dimlen(i_fileyid, i_dimid, i_lt)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read lat dimension')
          IF(i_lat /= i_lt) &
            CALL grid_error(c_error='[read_netcdf_currents]: lat dimension mismatch')
    
          i_ncstat= nf_inq_dimid(i_fileyid, 'time', i_dimid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not identify time dimension')
          i_ncstat= nf_inq_dimlen(i_fileyid, i_dimid, i_tm)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read time dimension')
          IF(i_timesteps /= i_tm) &
            CALL grid_error(c_error='[read_netcdf_currents]: time dimension mismatch')
    
          i_ncstat= nf_inq_dimid(i_fileyid, 'level', i_dimid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not identify time dimension')
          i_ncstat= nf_inq_dimlen(i_fileyid, i_dimid, i_ly)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read time dimension')
          IF(i_layer /= i_ly) &
            CALL grid_error(c_error='[read_netcdf_currents]: layer dimension mismatch')
    
          i_ncstat= nf_inq_dimid(i_filezid, 'lon', i_dimid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not identify lon dimension')
          i_ncstat= nf_inq_dimlen(i_filezid, i_dimid, i_ln)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read lon dimension')
          IF(i_lon /= i_ln) &
            CALL grid_error(c_error='[read_netcdf_currents]: lon dimension mismatch')
    
          i_ncstat= nf_inq_dimid(i_filezid, 'lat', i_dimid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not identify lat dimension')
          i_ncstat= nf_inq_dimlen(i_filezid, i_dimid, i_lt)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read lat dimension')
          IF(i_lat /= i_lt) &
            CALL grid_error(c_error='[read_netcdf_currents]: lat dimension mismatch')
    
          i_ncstat= nf_inq_dimid(i_filezid, 'time', i_dimid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not identify time dimension')
          i_ncstat= nf_inq_dimlen(i_filezid, i_dimid, i_tm)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read time dimension')
          IF(i_timesteps /= i_tm) &
            CALL grid_error(c_error='[read_netcdf_currents]: time dimension mismatch')
    
          i_ncstat= nf_inq_dimid(i_filezid, 'level', i_dimid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not identify time dimension')
          i_ncstat= nf_inq_dimlen(i_filezid, i_dimid, i_ly)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read time dimension')
          IF(i_ly .LT. p_control%tst%tst_int(1,1)) &
            CALL grid_error(c_error='[read_netcdf_currents]: w-layer dimension mismatch')
    
    !---------- allocate latitude and longitude coordinate arrays
    
          ALLOCATE(r_lat(i_lat), r_lon(i_lon), r_alt(i_layer), r_times(i_timesteps), stat= i_alct)
          IF(i_alct /= 0) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not allocate lat/lon/alt/time field')
    
    !---------- read latitude and longitude coordinate values
    
          i_ncstat= nf_inq_varid(i_filexid, 'lon', i_varid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not determine lon varid')
          i_ncstat= nf_get_var_real(i_filexid, i_varid, r_lon)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read lon data')
    
          i_ncstat= nf_inq_varid(i_filexid, 'lat', i_varid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not determine lat varid')
          i_ncstat= nf_get_var_real(i_filexid, i_varid, r_lat)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read lat data')
    
          i_ncstat= nf_inq_varid(i_filexid, 'level', i_varid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not determine level varid')
          i_ncstat= nf_get_var_real(i_filexid, i_varid, r_alt)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read level data')
    
          i_ncstat= nf_inq_varid(i_filexid, 'time', i_varid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not determine time varid')
          i_ncstat= nf_get_var_double(i_filexid, i_varid, r_times)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read time data')
    
    !---------- copy globally available pressure level information
    
          r_stclayerpressure(1:i_stcnumlayers) = r_alt(1:i_stcnumlayers)
    
    
    !---------- allocate current data arrays
          WRITE(GRID_parameters%ioout,*) &
            '***** reading NetCDF file core data ...               *****'
    
          ALLOCATE(r_flowx(i_lon, i_lat, i_layer, i_timesteps), &
                   r_flowy(i_lon, i_lat, i_layer, i_timesteps), &
                   r_flowz(i_lon, i_lat, i_layer, i_timesteps), &
                   r_auxx(i_timesteps, i_layer, i_lat, i_lon), &
                   r_auxy(i_timesteps, i_layer, i_lat, i_lon), &
                   r_auxz(i_timesteps, i_layer, i_lat, i_lon), stat= i_alct)
          IF(i_alct /= 0) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not allocate currents fields')
    
    !---------- read x-/y-direction data of currents
    
          i_ncstat= nf_inq_varid(i_filexid, 'uwnd', i_varid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not determine varid of uwnd')
          i_ncstat= nf_get_var_real(i_filexid, i_varid, r_auxx)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read uwnd data')
    
          i_ncstat= nf_inq_varid(i_fileyid, 'vwnd', i_varid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not determine varid of vwnd')
          i_ncstat= nf_get_var_real(i_fileyid, i_varid, r_auxy)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read vwnd data')
    
          i_ncstat= nf_inq_varid(i_filezid, 'omega', i_varid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not determine varid of omega')
          i_ncstat= nf_get_var_real(i_filezid, i_varid, r_auxz)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not read omega data')
    
    !---------- reshape fields for better accessibility
          WRITE(GRID_parameters%ioout,*) &
            '***** reshaping core data arrays ...                  *****'
    
          r_flowx = reshape(r_auxx, shape(r_flowx), order=(/4,3,2,1/))
          r_flowy = reshape(r_auxy, shape(r_flowy), order=(/4,3,2,1/))
          r_flowz = reshape(r_auxz, shape(r_flowz), order=(/4,3,2,1/))
          DEALLOCATE(r_auxx, r_auxy, r_auxz)
    
    !---------- Fix mask values
          WRITE(GRID_parameters%ioout,*) &
            '***** cleaning core data arrays ...                   *****'
    
          DO i_tm=1,i_timesteps
            DO i_ly= 1,i_layer
              DO i_lt= 1,i_lat
                DO i_ln= 1,i_lon
                  IF(r_flowx(i_ln,i_lt,i_ly,i_tm) .LE. -8.99E+30) &
                    r_flowx(i_ln,i_lt,i_ly,i_tm)= 0.0
                  IF(r_flowy(i_ln,i_lt,i_ly,i_tm) .LE. -8.99E+30) &
                    r_flowy(i_ln,i_lt,i_ly,i_tm)= 0.0
                  IF(r_flowz(i_ln,i_lt,i_ly,i_tm) .LE. -8.99E+30) &
                    r_flowz(i_ln,i_lt,i_ly,i_tm)= 0.0
                END DO
              END DO
            END DO
          END DO
    
    !---------- close currents files
    
          i_ncstat= nf_close(i_filexid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not close uwnd data file')
          i_ncstat= nf_close(i_fileyid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not close vwnd data file')
          i_ncstat= nf_close(i_filezid)
          IF(i_ncstat /= NF_NOERR) &
            CALL grid_error(c_error='[read_netcdf_currents]: could not close omega data file')
    
          WRITE(GRID_parameters%ioout,*) &
            '***** finished ingesting wind data                    *****'
          WRITE(GRID_parameters%ioout,*) &
            '***** ***** ***** ***** ***** ***** ***** ***** ***** *****'
    
    !      RETURN
          END SUBROUTINE read_netcdf_currents
    
    !*****************************************************************
    
          FUNCTION data_interpol(r_cart, i_layer, i_tim) RESULT (r_field)
    
    !---------- local declarations
    
          IMPLICIT NONE
          REAL (KIND = GRID_SR), DIMENSION(GRID_dimension)    :: r_cart
          INTEGER (KIND = GRID_SI)                            :: i_layer
          INTEGER (KIND = GRID_SI)                            :: i_tim
          REAL (KIND = GRID_SR), DIMENSION(GRID_dimension)    :: r_field
    
          REAL (KIND = GRID_SR), DIMENSION(GRID_dimspherical) :: r_inter
          REAL (KIND = GRID_SR), DIMENSION(GRID_dimspherical) :: r_coord
          INTEGER (KIND = GRID_SI)                            :: i_lox, i_hix, i_loy, i_hiy, i_cnt
          REAL (KIND = GRID_SR)                               :: r_dx, r_dy, r_lx, r_hx, &
            r_ly, r_hy, r_l1, r_h1, r_l2, r_h2, r_dxi, r_dyi
          REAL (KIND = GRID_SR)                               :: r_scalx, r_scaly
          REAL (KIND = GRID_SR)                               :: r_deg
          REAL (KIND = GRID_SR)                               :: r_earth
         
    !---------- initialize radians
    
          r_deg= 360._GRID_SR/(r_earthrad*2.*GRID_PI)
         
    !---------- compute lat/lon coordinates from Cartesian (given) coordinates
    
          r_coord = grid_kartgeo(r_cart)
         
    !---------- correct lon coordinates to 0..360 degrees
    
          IF(r_coord(1) .LT. 0_GRID_SR) r_coord(1)= r_coord(1)+ 360._GRID_SR
    
    !---------- find wind box corresponding to coordinate
    
          i_lox=1_GRID_SI
          i_hix=2_GRID_SI
          i_loy=1_GRID_SI
          i_hiy=2_GRID_SI
    
          determine_lon: DO i_cnt=1, i_lon-1
            IF(REAL(r_lon(i_cnt),GRID_SR) .LT. r_coord(1)) THEN
              IF(REAL(r_lon(i_cnt +1),GRID_SR) .GE. r_coord(1)) THEN
                i_lox= i_cnt
                i_hix= i_cnt+1
                exit determine_lon
              END IF
            END IF
          END DO determine_lon
          IF(REAL(r_lon(i_lon),GRID_SR) .LT. r_coord(1)) THEN
            i_lox = i_lon-1
            i_hix = i_lon
          END IF
    
          determine_lat: DO i_cnt=1, i_lat-1
            IF(REAL(r_lat(i_cnt),GRID_SR) .GE. r_coord(2)) THEN
              IF(REAL(r_lat(i_cnt +1),GRID_SR) .LT. r_coord(2)) THEN
                i_loy= i_cnt
                i_hiy= i_cnt+1
                exit determine_lat
              END IF
            END IF
          END DO determine_lat
          IF(REAL(r_lat(i_lat),GRID_SR) .LT. r_coord(2)) THEN
            i_loy = i_lat-1
            i_hiy = i_lat
          END IF
          
    
    !---------- calculate weights for bilinear interpolation
    !           NOTE: We do not consider the distance stretching with height!
    
          r_dx= REAL(r_lon(i_hix),GRID_SR)- REAL(r_lon(i_lox),GRID_SR)
          r_dy= REAL(r_lat(i_hiy),GRID_SR)- REAL(r_lat(i_loy),GRID_SR)
          r_lx= r_coord(1) - REAL(r_lon(i_lox),GRID_SR)
          r_ly= r_coord(2) - REAL(r_lat(i_loy),GRID_SR)
          r_hx= REAL(r_lon(i_hix),GRID_SR) - r_coord(1)
          r_hy= REAL(r_lat(i_hiy),GRID_SR) - r_coord(2)
          r_scalx=r_deg
          r_scaly=r_deg*cos(r_coord(2))
    
    !---------- linear interpolation in x-direction
    
          IF(r_dx /= 0.0_GRID_SR) THEN
            r_dxi= 1._GRID_SR/r_dx
            r_l1= (r_hx* REAL(r_flowx(i_lox, i_loy, i_layer, i_tim),GRID_SR)+ &
                   r_lx* REAL(r_flowx(i_hix, i_loy, i_layer, i_tim),GRID_SR))* r_dxi
            r_h1= (r_hx* REAL(r_flowx(i_lox, i_hiy, i_layer, i_tim),GRID_SR)+ &
                   r_lx* REAL(r_flowx(i_hix, i_hiy, i_layer, i_tim),GRID_SR))* r_dxi
            r_l2= (r_hx* REAL(r_flowy(i_lox, i_loy, i_layer, i_tim),GRID_SR)+ &
                   r_lx* REAL(r_flowy(i_hix, i_loy, i_layer, i_tim),GRID_SR))* r_dxi
            r_h2= (r_hx* REAL(r_flowy(i_lox, i_hiy, i_layer, i_tim),GRID_SR)+ &
                   r_lx* REAL(r_flowy(i_hix, i_hiy, i_layer, i_tim),GRID_SR))* r_dxi
          ELSE
            r_l1= REAL(r_flowx(i_lox, i_loy, i_layer, i_tim),GRID_SR)
            r_h1= REAL(r_flowx(i_lox, i_hiy, i_layer, i_tim),GRID_SR)
            r_l2= REAL(r_flowy(i_lox, i_loy, i_layer, i_tim),GRID_SR)
            r_h2= REAL(r_flowy(i_lox, i_hiy, i_layer, i_tim),GRID_SR)
          END IF
    
    !---------- linear interpolation in y-direction
    
          IF(r_dy /= 0.0_GRID_SR) THEN
            r_dyi= 1._GRID_SR/r_dy
            r_inter(1)= (r_hy* r_l1+ r_ly* r_h1)* r_dyi* r_scalx
            r_inter(2)= (r_hy* r_l2+ r_ly* r_h2)* r_dyi* r_scaly
          ELSE
            r_inter(1)= r_l1* r_scalx
            r_inter(2)= r_l2* r_scaly
          END IF
    
    !---------- calculate wind vector in kartesian coordinates
    
    	  r_field= kartwind(r_inter,r_coord)
    
          RETURN
          END FUNCTION data_interpol
    
    !*****************************************************************
          FUNCTION vertical_interpol(r_cart, i_layer, i_tim) RESULT (r_inter)
    
    !---------- local declarations
    
          IMPLICIT NONE
          REAL (KIND = GRID_SR), DIMENSION(GRID_dimension) :: r_cart
          INTEGER (KIND = GRID_SI)                         :: i_layer
          INTEGER (KIND = GRID_SI)                         :: i_tim
          REAL (KIND = GRID_SR)                            :: r_inter
    
          REAL (KIND = GRID_SR), DIMENSION(GRID_dimspherical) :: r_coord
          INTEGER (KIND = GRID_SI)                         :: i_lox, i_hix, i_loy, i_hiy, i_cnt
          REAL (KIND = GRID_SR)                            :: r_dx, r_dy, r_lx, r_hx, &
            r_ly, r_hy, r_l1, r_h1, r_l2, r_h2, r_dxi, r_dyi
          REAL (KIND = GRID_SR)                            :: r_scalx, r_scaly
          REAL (KIND = GRID_SR)                            :: r_deg
          REAL (KIND = GRID_SR)                            :: r_earth
              
    !---------- compute lat/lon coordinates from Cartesian (given) coordinates
    
          r_coord = grid_kartgeo(r_cart)
    
    !---------- find wind box corresponding to coordinate
    
          i_lox=1_GRID_SI
          i_hix=2_GRID_SI
          i_loy=1_GRID_SI
          i_hiy=2_GRID_SI
    
          determine_lon: DO i_cnt=1, i_lon-1
            IF(REAL(r_lon(i_cnt),GRID_SR) .LT. r_coord(1)) THEN
              IF(REAL(r_lon(i_cnt +1),GRID_SR) .GE. r_coord(1)) THEN
                i_lox= i_cnt
                i_hix= i_cnt+1
                exit determine_lon
              END IF
            END IF
          END DO determine_lon
          IF(REAL(r_lon(i_lon),GRID_SR) .LT. r_coord(1)) THEN
            i_lox = i_lon-1
            i_hix = i_lon
          END IF
    
          determine_lat: DO i_cnt=1, i_lat-1
            IF(REAL(r_lat(i_cnt),GRID_SR) .LT. r_coord(2)) THEN
              IF(REAL(r_lat(i_cnt +1),GRID_SR) .GE. r_coord(2)) THEN
                i_loy= i_cnt
                i_hiy= i_cnt+1
                exit determine_lat
              END IF
            END IF
          END DO determine_lat
          IF(REAL(r_lat(i_lat),GRID_SR) .LT. r_coord(2)) THEN
            i_loy = i_lat-1
            i_hiy = i_lat
          END IF
          
    
    !---------- calculate weights for bilinear interpolation
    !           NOTE: We do not consider the distance stretching with height!
    
          r_dx= REAL(r_lon(i_hix),GRID_SR)- REAL(r_lon(i_lox),GRID_SR)
          r_dy= REAL(r_lat(i_hiy),GRID_SR)- REAL(r_lat(i_loy),GRID_SR)
          r_lx= r_coord(1) - REAL(r_lon(i_lox),GRID_SR)
          r_ly= r_coord(2) - REAL(r_lat(i_loy),GRID_SR)
          r_hx= REAL(r_lon(i_hix),GRID_SR) - r_coord(1)
          r_hy= REAL(r_lat(i_hiy),GRID_SR) - r_coord(2)
    
    !---------- linear interpolation in x-direction
    
          IF(r_dx /= 0.0_GRID_SR) THEN
            r_dxi= 1._GRID_SR/r_dx
            r_l1= (r_hx* REAL(r_flowz(i_lox, i_loy, i_layer, i_tim),GRID_SR)+ &
                   r_lx* REAL(r_flowz(i_hix, i_loy, i_layer, i_tim),GRID_SR))* r_dxi
            r_h1= (r_hx* REAL(r_flowz(i_lox, i_hiy, i_layer, i_tim),GRID_SR)+ &
                   r_lx* REAL(r_flowz(i_hix, i_hiy, i_layer, i_tim),GRID_SR))* r_dxi
          ELSE
            r_l1= REAL(r_flowz(i_lox, i_loy, i_layer, i_tim),GRID_SR)
            r_h1= REAL(r_flowz(i_lox, i_hiy, i_layer, i_tim),GRID_SR)
          END IF
    
    !---------- linear interpolation in y-direction
    
          IF(r_dy /= 0.0_GRID_SR) THEN
            r_dyi= 1._GRID_SR/r_dy
            r_inter= (r_hy* r_l1+ r_ly* r_h1)* r_dyi
          ELSE
            r_inter= r_l1
          END IF
    
          RETURN
          END FUNCTION vertical_interpol
    
    !*****************************************************************
          FUNCTION time_from_datecomponents(i_year, i_month, i_day, i_hour, i_minute, &
                                            i_second, r_start, l_leap) RESULT (r_longtime)
    
    !---------- local declarations
    
          IMPLICIT NONE
          INTEGER (KIND = GRID_SI), INTENT(in)  :: i_year
          INTEGER (KIND = GRID_SI), INTENT(in)  :: i_month
          INTEGER (KIND = GRID_SI), INTENT(in)  :: i_day
          INTEGER (KIND = GRID_SI), INTENT(in)  :: i_hour
          INTEGER (KIND = GRID_SI), INTENT(in)  :: i_minute
          INTEGER (KIND = GRID_SI), INTENT(in)  :: i_second
          REAL (KIND = GRID_SR), INTENT(in)     :: r_start
          LOGICAL, INTENT(in), OPTIONAL         :: l_leap
          REAL (KIND = GRID_SR)                 :: r_longtime
          
          INTEGER (KIND = GRID_SI), DIMENSION(12) :: i_montharray
          LOGICAL                                 :: l_leapyear
          REAL (KIND = GRID_SR)                   :: r_accumhours
          INTEGER (KIND = GRID_SI)                :: i_cnt, i_current
         
    !---------- conversion function -- DOES NOT WORK (2022-01-25, JB)
    ! reference: https://wiki.usask.ca/display/MESH/How+to+manually+calculate+NetCDF+times
    
    !       r_longtime = real((1461*(i_year + 4800 + (i_month - 14)/12))/4 + &
    !                        (367*(i_month - 2 - 12*((i_month - 14)/12)))/12 - &
    !                        (3*((i_year + 4900 + (i_month - 14)/12)/100))/4 + i_day - 32075) + &
    !                        i_hour/24.0 + i_minute/24.0/60.0 + i_second/24.0/60.0/60.0
    
    !---------- handle optional leap year variable
    
          IF(present(l_leap)) THEN
            l_leapyear= l_leap
          ELSE
            l_leapyear= .FALSE.
          END IF
    
    !---------- define month's days
    
          IF(l_leapyear) THEN
            i_montharray = (/ 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
          ELSE
            i_montharray = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
          END IF
    
    !---------- compute accumulated hours from start for the months
    
          r_accumhours = 0.
          i_current    = 0
          DO i_cnt=1,i_month
            r_accumhours = r_accumhours+ 24. * i_current
            i_current = i_montharray(i_cnt)
          END DO
    
    !---------- compute accumulated hours from days, hours, minutes and seconds
          
          r_accumhours = r_accumhours+ 24. * (i_day-1) + i_hour + i_minute/60 + i_second/3600
    
    !---------- add offset
    
          r_longtime = r_accumhours + r_start
    
          RETURN
          END FUNCTION time_from_datecomponents
    
    !*****************************************************************
          SUBROUTINE datecomponents_from_time(r_longtime, r_start, i_year, i_month, i_day, i_hour, &
                                              i_minute, i_second, l_leap)
    
    !---------- local declarations
    
          IMPLICIT NONE
          REAL (KIND = GRID_SR), INTENT(in)       :: r_longtime
          REAL (KIND = GRID_SR), INTENT(in)       :: r_start
          INTEGER (KIND = GRID_SI), INTENT(inout) :: i_year
          INTEGER (KIND = GRID_SI), INTENT(out)   :: i_month
          INTEGER (KIND = GRID_SI), INTENT(out)   :: i_day
          INTEGER (KIND = GRID_SI), INTENT(out)   :: i_hour
          INTEGER (KIND = GRID_SI), INTENT(out)   :: i_minute
          INTEGER (KIND = GRID_SI), INTENT(out)   :: i_second
          LOGICAL, OPTIONAL                       :: l_leap
        
          INTEGER (KIND = GRID_SI)               :: i_auxa, i_auxb, i_auxc
          
          INTEGER (KIND = GRID_SI), DIMENSION(12) :: i_montharray
          LOGICAL                                 :: l_leapyear
          REAL (KIND = GRID_SR)                   :: r_accumhours, r_time, r_remhours, r_hours
          INTEGER (KIND = GRID_SI)                :: i_cnt, i_current
    
    !---------- conversion function -- DOES NOT WORK (2022-01-25, JB)
    ! reference: https://wiki.usask.ca/display/MESH/How+to+manually+calculate+NetCDF+times
    
    !       i_auxa = 4*(floor(r_longtime) + 1401 + (((4*floor(r_longtime) + &
    !                274277)/146097)*3)/4 - 38) + 3
    !       i_auxb = mod(i_auxa, 1461)/4
    !       i_auxc = 5*i_auxb + 2
    ! 
    !       i_day = (mod(i_auxc, 153))/5 + 1
    !       i_month = mod((i_auxc/153 + 2), 12) + 1
    !       i_year = i_auxa/1461 - 4716 + (12 + 2 - i_month)/12
    !       i_second = mod(nint((r_longtime - floor(r_longtime))*24.0*60.0*60.0), 60)
    !       i_minute = mod(nint((r_longtime - floor(r_longtime))*24.0*60.0 - i_second/60.0), 60)
    !       i_hour = nint((r_longtime - floor(r_longtime))*24.0 - i_minute/60.0 - i_second/60.0/60.0)
    
    !---------- handle optional leap year variable
    
          IF(present(l_leap)) THEN
            l_leapyear= l_leap
          ELSE
            l_leapyear= .FALSE.
          END IF
    
    !---------- define month's days
    
          IF(l_leapyear) THEN
            i_montharray = (/ 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
          ELSE
            i_montharray = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
          END IF
    
    !---------- compute month
    
          r_time       = r_longtime - r_start
          r_accumhours = 0.
          i_current    = 0
          i_cnt        = 1
          DO WHILE ((r_accumhours .LT. r_time) .AND. (i_cnt .LE. 12))
            i_current    = i_montharray(i_cnt)
            r_accumhours = r_accumhours+ 24. * i_current
            i_cnt        = i_cnt + 1
          END DO
          i_month = i_cnt
    
    !---------- from remaining hours compute days
    
          r_remhours = r_time - (24. * sum(i_montharray(1:i_month)))
          i_day      = floor(r_remhours/24) + 1
          r_hours    = modulo(r_remhours,24._GRID_SR)
    
    !---------- now hours, minutes, seconds
    
          i_hour     = floor(r_hours)
          r_hours    = r_hours - REAL(i_hour,GRID_SR)
          
          i_minute   = floor(r_hours*60.)
          r_hours    = r_hours*60. - REAL(i_minute,GRID_SR)
          
          i_second   = floor(r_hours*60.)
    
          END SUBROUTINE datecomponents_from_time
    
    !*****************************************************************
    	  FUNCTION kartwind(r_lamphi, r_lpcoor) RESULT (r_kart)
    
    !---------- local declarations
    
    	  IMPLICIT NONE
    	  REAL (KIND = GRID_SR), DIMENSION(GRID_dimspherical) :: r_lamphi, r_lpcoor
    	  REAL (KIND = GRID_SR), DIMENSION(GRID_dimension)    :: r_kart
    	  REAL (KIND = GRID_SR)                               :: r_sl, r_cl, r_sp, r_cp
    
    !---------- calculate constants
    
    	  r_sl= sin(r_lpcoor(1))
    	  r_cl= cos(r_lpcoor(1))
    	  r_sp= sin(r_lpcoor(2))
    	  r_cp= cos(r_lpcoor(2))
    
    !---------- calculate 3D kartesian components
    
    	  r_kart(1)= -r_sl* r_lamphi(1)- r_sp* r_cl* r_lamphi(2)
    	  r_kart(2)=  r_cl* r_lamphi(1)- r_sp* r_sl* r_lamphi(2)
    	  r_kart(3)=                     r_cp*       r_lamphi(2)
    
    	  RETURN
    	  END FUNCTION kartwind
    
    !*****************************************************************
        END MODULE ADV_wind