!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
module current_meters,1
!BOP
! !MODULE: current_meters
!
! !DESCRIPTION:
! This module collects data to compare with current meter data.
! NOTE: this module currently does not work. old CM-5 routines
! are appended but must be re-done.
!
! !REVISION HISTORY:
! SVN:$Id: current_meters.F90 808 2006-04-28 17:06:38Z njn01 $
! !USES:
use kinds_mod
implicit none
private
save
!EOP
!BOC
!EOC
!***********************************************************************
contains
!***********************************************************************
!BOP
! !IROUTINE: init_current_meters
! !INTERFACE:
subroutine init_current_meters
! !DESCRIPTION:
! Initializes all necessary variables for current meter diagnostics.
!
! !REVISION HISTORY:
! same as module
!EOP
!BOC
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!EOC
end subroutine init_current_meters
!***********************************************************************
end module current_meters
!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#ifdef old_current_meters
integer num_cmeters_max, num_buoys_max, num_cmeters
parameter ( num_cmeters_max = 1, num_buoys_max = 1100)
integer nnbr_buoys
parameter ( nnbr_buoys = 4 )
integer num_buoys(num_cmeters_max)
TYPE buoy_xy(num_buoys_max,num_cmeters_max,2)
character*80 cmeter_file(num_cmeters_max)
common/hydro_buoys_scalar/ num_cmeters , num_buoys, buoy_xy
common/hydro_buoys_char/ cmeter_file
integer, dimension(num_buoys_max, num_cmeters_max
& , nnbr_buoys, 2) :: ADDR_BUOYS
TYPE, dimension(num_buoys_max, num_cmeters_max
& , nnbr_buoys) :: DIST_BUOYS
common/hydro_buoys_array/ ADDR_BUOYS, DIST_BUOYS
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c begin file cmeters.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
subroutine init_cmeters,8
c-----------------------------------------------------------------------
c Initialize everything necessary for writing data for comparison
c with hydrographic sections.
c-----------------------------------------------------------------------
use grid
use io
use time_management
implicit none
integer (kind=int_kind) ::
& nu ! i/o unit attached to input file
integer syear, smonth, sday, start_year, cmeter, buoy, nlen
integer lenc, i, num_temp
external lenc
real lat,lon,lat1,lon1
character*80 temp_file
double precision pos_buoys(num_buoys_max,2)
integer, dimension(num_buoys_max,nnbr_buoys,2) :: ADDR_TEMP
double precision, dimension(num_buoys_max,nnbr_buoys) :: DIST_TEMP
c****************************************************************
c-----------------------------------------------------------------------
c Read in names of the current meter files used to get times and positions.
c-----------------------------------------------------------------------
call get_unit(nu)
open(nu,file = in_cmeters, status = 'old')
num_cmeters = 0
10 continue
if(num_cmeters .gt. num_cmeters_max) then
write(stdout,*)' '
write(stdout,*)' Too many current meter files listed in file '
& ,in_cmeters
write(stdout,*)' Increase the value of num_cmeters_max'
stop
endif
c-----------------------------------------------------------------------
c Read in_cmeters file until an error (or end-of-file) occurs
c so the user does not have to specify the number of files to
c be read in. The file name goes into the cmeter_file array.
c Also note that if we encounter 'none' anywhere,
c then no files are read in.
c-----------------------------------------------------------------------
read(nu,'(a)',err = 20) temp_file
if(temp_file .eq. 'none') then
write(stdout,*)' '
write(stdout,*)' No current meter files will be read'
num_cmeters = 0
goto 99
endif
num_cmeters = num_cmeters + 1
nlen = lenc(temp_file)
cmeter_file(num_cmeters) = temp_file(1:nlen)
goto 10
20 continue
write(stdout,*)' '
write(stdout,"(' Attempting to read in ',i4,
& ' current meter files')") num_cmeters
close(nu)
call release_unit(nu)
c-----------------------------------------------------------------------
c Now read in data from each current meter file. Give warning if there
c are more buoys on a given current meter than allocated space.
c
c Current meter files generally have several instruments at the
c same location (different depths) so we are actually interested
c in only one horizontal location since well dump velocity at all
c depths.
c-----------------------------------------------------------------------
write(stdout,*)' '
do cmeter = 1,num_cmeters
temp_file = '/home/ocean/data/cmeters/'/ /cmeter_file(cmeter)
call get_unit(nu)
open(nu,file = temp_file, status = 'old')
write(stdout,*)' Reading file: ',cmeter_file(cmeter)
read(nu,*)num_temp
lat = -999.
lon = -999.
num_buoys(cmeter) = 0
do buoy = 1,num_temp
read(nu,11)lat1,lon1
if(lat1.ne.lat .or. lon1.ne.lon) then
if(num_buoys(cmeter).eq.num_buoys_max) then
write(stdout,*)' '
write(stdout,*)' *** Warning ***'
write(stdout,*)
& ' More buoys on this current meter file than allowed'
write(stdout,"(' Proceeding with first ',i4,' buoys')")
& num_buoys_max
num_buoys(cmeter) = num_buoys_max
goto 22
endif
num_buoys(cmeter) = num_buoys(cmeter) + 1
lat = lat1
lon = lon1
buoy_xy(num_buoys(cmeter),cmeter,1) = lon/radian
buoy_xy(num_buoys(cmeter),cmeter,2) = lat/radian
endif
enddo
22 continue
close(nu)
call release_unit(nu)
enddo
11 format(f7.3,1x,f8.3,3(1x,i2),1x,3(1x,i2),1x,a1,1x,a)
c-----------------------------------------------------------------------
c Compute the neighbor addresses for each buoy.
c-----------------------------------------------------------------------
do cmeter = 1,num_cmeters
do buoy = 1,num_buoys(cmeter)
pos_buoys(buoy,1) = buoy_xy(buoy,cmeter,1)
pos_buoys(buoy,2) = buoy_xy(buoy,cmeter,2)
enddo
call gather_set(ADDR_TEMP, nnbr_buoys, num_buoys(cmeter),
& num_buoys_max, ULAT, ULONG, CALCU, pos_buoys, DIST_TEMP)
ADDR_BUOYS(:,cmeter,:,:) = ADDR_TEMP
DIST_BUOYS(:,cmeter,: ) = DIST_TEMP
enddo
99 continue
return
end subroutine init_cmeters
c**************************************************
subroutine data_cmeters,8
c-----------------------------------------------------------------------
c Gather neighboring data and dump to file for all buoys.
c-----------------------------------------------------------------------
use grid
use time_management
use prognostic
implicit none
integer cmeter, buoy, nlen, lenc, k, n, nf, nu, iostat, iml, imr
external lenc
character*80 temp_file
character*10 cday
integer, dimension(num_buoys_max,nnbr_buoys,2) :: ADDR_TEMP
double precision, dimension(num_buoys_max,nnbr_buoys,km)
& :: WORK1, WORK2, WORK3
real, dimension(num_buoys_max,nnbr_buoys,km,3)
& :: FIELDS_BUOYS
double precision, dimension(imt,jmt,km) :: FIELD
integer, dimension(512) :: ITEMP
real, dimension(num_buoys_max,nnbr_buoys) :: DIST_TEMP
real, dimension(num_buoys_max,2) :: XY_TEMP
c*********************************************
c-----------------------------------------------------------------------
c Loop over all buoys and check to see where data is needed.
c Note that U and V are in local coordinates, so we must rotate
c them to true zonal and meridional components for non-polar grids.
c-----------------------------------------------------------------------
do cmeter = 1,num_cmeters
ADDR_TEMP(:,:,:) = ADDR_BUOYS(:,cmeter,:,:)
WORK1 = c0 ! initialize
FIELD = U
call gather(WORK1, FIELD, ADDR_TEMP, nnbr_buoys
& , num_buoys(cmeter))
FIELD = V
call gather(WORK2, FIELD, ADDR_TEMP, nnbr_buoys
& , num_buoys(cmeter))
#if polar_grid
FIELDS_BUOYS(:,:,:,1) = WORK1
FIELDS_BUOYS(:,:,:,2) = WORK2
#else
c-----------------------------------------------------------------------
c Rotate velocity vector to lat-long grid for non-polar grids.
c-----------------------------------------------------------------------
do k = 1,km
FIELD(:,:,k) = ANGLE ! redundant but easy
enddo
call gather(WORK3, FIELD, ADDR_TEMP, nnbr_buoys
& , num_buoys(cmeter))
FIELDS_BUOYS(:,:,:,1) = WORK1*cos(WORK3) + WORK2*sin(-WORK3)
FIELDS_BUOYS(:,:,:,2) = WORK2*cos(WORK3) - WORK1*sin(-WORK3)
#endif
call wcalc(FIELD,U,V,this_block) ! FIELD contains W
call gather(WORK1, FIELD, ADDR_TEMP, nnbr_buoys
& , num_buoys(cmeter))
FIELDS_BUOYS(:,:,:,3) = WORK1
write (cday,'(i10)') iday + 1000000000 ! allows 2.7 million years
nlen = lenc(cmeter_file(cmeter))
iml = 10 - movie_digits + 1
imr = lenc(out_cmeters)
temp_file = out_cmeters(1:imr)/ /cmeter_file(cmeter)(1:nlen)
& / /'.'/ /cday(iml:10)
nu = 31
call CMF_file_open(nu,temp_file,iostat)
write(stdout,*)' '
write(stdout,*)' Writing file: ',temp_file
ITEMP(1) = num_buoys_max
ITEMP(2) = num_buoys(cmeter)
ITEMP(3) = nnbr_buoys
ITEMP(4) = km
DIST_TEMP = -999.
do k = 1,nnbr_buoys
do n = 1,num_buoys(cmeter)
DIST_TEMP(n,k) = DIST_BUOYS(n,cmeter,k)
enddo
enddo
XY_TEMP = -999.
do n = 1,num_buoys(cmeter)
XY_TEMP(n,1) = buoy_xy(n,cmeter,1)
XY_TEMP(n,2) = buoy_xy(n,cmeter,2)
enddo
call CMF_cm_array_to_file(nu,ITEMP,iostat)
if (iostat.eq.-1) then
if(my_task.eq.master_task)
& write(stdout,*) ' i/o error writing ',ITEMP
call write_status(iostat)
stop
endif
call CMF_cm_array_to_file(nu,DIST_TEMP,iostat)
if (iostat.eq.-1) then
if(my_task.eq.master_task)
& write(stdout,*) ' i/o error writing ',DIST_TEMP
call write_status(iostat)
stop
endif
call CMF_cm_array_to_file(nu,XY_TEMP,iostat)
if (iostat.eq.-1) then
if(my_task.eq.master_task)
& write(stdout,*) ' i/o error writing ',XY_TEMP
call write_status(iostat)
stop
endif
call CMF_cm_array_to_file(nu,FIELDS_BUOYS,iostat)
if (iostat.eq.-1) then
if(my_task.eq.master_task)
& write(stdout,*) ' i/o error writing ',FIELDS_BUOYS
call write_status(iostat)
stop
endif
call CMF_file_close(nu,iostat)
enddo
end subroutine data_cmeters
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c end file cmeters.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c begin file gather.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c***********************************************************************
c this set of routines searches for and gathers the nnbr closest
c points on an ocean grid to a given set of target points.
c
c written by: Phil Jones, T-3, Los Alamos National Laboratory
c date last revised: 6 December 1994
c***********************************************************************
subroutine gather(DSTARR, SRCARR, IADDR, nnbr, nlist) 16,4
c-----------------------------------------------------------------------
c this routine gathers nnbr values from srcarr and places the
c result in dstarr given a list of addresses computed in a
c previous setup routine (gather_set).
c-----------------------------------------------------------------------
implicit none
c-----------------------------------------------------------------------
c intent (in):
c
c nlist = number of search points
c
c nnbr = number of neighboring points to gather for each
c target point
c
c SRCARR = the ocean variable to gather
c
c IADDR = the list of gather addresses computed in
c a previous setup routine
c-----------------------------------------------------------------------
integer nlist, nnbr, k
integer, dimension(nlist, nnbr, 2) :: IADDR
real, dimension(imt, jmt, km) :: SRCARR
c-----------------------------------------------------------------------
c intent(out):
c
c DSTARR = values at all the neighboring ocean points
c-----------------------------------------------------------------------
real, dimension(nlist, nnbr, km) :: DSTARR
c-----------------------------------------------------------------------
c local variables:
c-----------------------------------------------------------------------
integer i, n
c-----------------------------------------------------------------------
c gather up the nnbr neighbors for all model levels.
c-----------------------------------------------------------------------
DSTARR = 0.0
do k = 1,km
forall (n=1:nlist, i=1:nnbr)
& DSTARR(n,i,k) = SRCARR(IADDR(n,i,1),IADDR(n,i,2),k)
enddo
c-----------------------------------------------------------------------
return
end subroutine gather
c***********************************************************************
c***********************************************************************
c***********************************************************************
subroutine gather_set(IADDR, nnbr, nlist, nlmax, OLAT, OLON, OMSK, 3
& plist, NBR_DIST)
c-----------------------------------------------------------------------
c this subroutine sets up address arrays for the nnbr closests
c points to a set of target points.
c-----------------------------------------------------------------------
implicit none
c-----------------------------------------------------------------------
c intent(in):
c
c nlist = number of points in list of search points
c
c nnbr = number of neighboring points to gather for
c each target point
c
c OLAT, OLON = array of latitudes and longitudes on the
c ocean grid
c
c OMSK = ocean mask
c
c plist = list of points for which neighboring ocean
c points are desired
c-----------------------------------------------------------------------
integer nlist, nnbr, nlmax
logical, dimension(imt, jmt) :: OMSK
real, dimension(imt, jmt) :: OLAT, OLON
real, dimension(nlmax, 2) :: plist
c-----------------------------------------------------------------------
c intent(out):
c
c IADDR = the addresses of the nnbr closest ocean points to the
c list of search points
c-----------------------------------------------------------------------
integer, dimension(nlmax, nnbr, 2) :: IADDR
real, dimension(nlmax, nnbr) :: NBR_DIST
c-----------------------------------------------------------------------
c local variables:
c-----------------------------------------------------------------------
integer i, j, n
real *8 rlat, rlon, fac1, fac2, fac3, dmin, wttmp, dsum
real, dimension(imt, jmt) :: OARC1, OARC2, OARC3, DIST
c-----------------------------------------------------------------------
c set up some arc-length constants and initialize some stuff.
c-----------------------------------------------------------------------
OARC3 = COS(OLAT)
OARC1 = COS(OLON)*OARC3
OARC2 = SIN(OLON)*OARC3
OARC3 = SIN(OLAT)
IADDR = 0
c-----------------------------------------------------------------------
c loop through the list of destination points. compute the
c arclength from this point to all the ocean grid points.
c-----------------------------------------------------------------------
do n=1,nlist
rlon = plist(n, 1)
rlat = plist(n, 2)
fac3 = COS(rlat)
fac1 = COS(rlon)*fac3
fac2 = SIN(rlon)*fac3
fac3 = SIN(rlat)
DIST = ACOS(OARC1*fac1 + OARC2*fac2 + OARC3*fac3)
c-----------------------------------------------------------------------
c find the closest nnbr points. eliminate land points.
c-----------------------------------------------------------------------
where (.NOT. OMSK) DIST = 10.0
do i=1,nnbr
IADDR(n,i,:) = MINLOC(DIST)
NBR_DIST(n,i) = DIST(IADDR(n,i,1),IADDR(n,i,2))
DIST(IADDR(n,i,1),IADDR(n,i,2)) = 10.0
end do
end do
c-----------------------------------------------------------------------
return
end subroutine gather_set
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
c end file gather.F
c|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
#endif