Select Git revision
ADV_wind_stacked.f90
Behrens, Prof. Dr. Jörn authored
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