Commit 4475cd3d authored by kbk's avatar kbk
Browse files

flux calc. for point source + combined rot. met. and grid convergence

parent 90c7dc52
!$Id: ncdf_meteo.F90,v 1.15 2005-01-12 19:26:16 kbk Exp $ !$Id: ncdf_meteo.F90,v 1.16 2005-03-31 10:14:20 kbk Exp $
#include "cppdefs.h" #include "cppdefs.h"
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
!BOP !BOP
...@@ -12,9 +12,10 @@ ...@@ -12,9 +12,10 @@
! !
! !USES: ! !USES:
use time, only: string_to_julsecs,time_diff,add_secs,in_interval use time, only: string_to_julsecs,time_diff,add_secs,in_interval
use time, only: jul0,secs0,julianday,secondsofday,timestep use time, only: jul0,secs0,julianday,secondsofday,timestep,simtime
use domain, only: imin,imax,jmin,jmax,az,lonc,latc,conv use domain, only: imin,imax,jmin,jmax,az,lonc,latc,conv
use grid_interpol, only: init_grid_interpol,do_grid_interpol use grid_interpol, only: init_grid_interpol,do_grid_interpol
use grid_interpol, only: to_rotated_lat_lon
use meteo, only: meteo_file,on_grid,calc_met,method,hum_method use meteo, only: meteo_file,on_grid,calc_met,method,hum_method
use meteo, only: airp,u10,v10,t2,hum,tcc use meteo, only: airp,u10,v10,t2,hum,tcc
use meteo, only: tausx,tausy,swr,shf use meteo, only: tausx,tausy,swr,shf
...@@ -36,6 +37,7 @@ ...@@ -36,6 +37,7 @@
integer :: iextr,jextr,textr,tmax=-1 integer :: iextr,jextr,textr,tmax=-1
integer :: grid_scan=1 integer :: grid_scan=1
logical :: point_source=.false. logical :: point_source=.false.
logical :: rotated_meteo_grid=.false.
REALTYPE, allocatable :: met_lon(:),met_lat(:) REALTYPE, allocatable :: met_lon(:),met_lat(:)
REAL_4B, allocatable :: met_times(:) REAL_4B, allocatable :: met_times(:)
...@@ -77,7 +79,10 @@ ...@@ -77,7 +79,10 @@
! Original author(s): Karsten Bolding & Hans Burchard ! Original author(s): Karsten Bolding & Hans Burchard
! !
! $Log: ncdf_meteo.F90,v $ ! $Log: ncdf_meteo.F90,v $
! Revision 1.15 2005-01-12 19:26:16 kbk ! Revision 1.16 2005-03-31 10:14:20 kbk
! flux calc. for point source + combined rot. met. and grid convergence
!
! Revision 1.15 2005/01/12 19:26:16 kbk
! fixed printing of south pole ! fixed printing of south pole
! !
! Revision 1.14 2005/01/12 19:17:47 kbk ! Revision 1.14 2005/01/12 19:17:47 kbk
...@@ -179,7 +184,7 @@ ...@@ -179,7 +184,7 @@
integer :: i,j,n integer :: i,j,n
integer :: err integer :: err
logical :: ok=.true. logical :: ok=.true.
REALTYPE :: x REALTYPE :: olon,olat,rlon,rlat,x
!EOP !EOP
!------------------------------------------------------------------------- !-------------------------------------------------------------------------
include "netcdf.inc" include "netcdf.inc"
...@@ -191,15 +196,6 @@ ...@@ -191,15 +196,6 @@
call open_meteo_file(meteo_file) call open_meteo_file(meteo_file)
if(iextr .eq. 1 .and. jextr .eq. 1) then
point_source = .true.
LEVEL3 'Assuming Point Source meteo forcing'
if (on_grid .eq. .false. ) then
LEVEL3 'Setting on_grid to true'
on_grid=.true.
end if
end if
if (met_lat(1) .gt. met_lat(2)) then if (met_lat(1) .gt. met_lat(2)) then
LEVEL3 'Reverting lat-axis and setting grid_scan to 0' LEVEL3 'Reverting lat-axis and setting grid_scan to 0'
grid_scan = 0 grid_scan = 0
...@@ -219,12 +215,27 @@ ...@@ -219,12 +215,27 @@
if (err /= 0) stop 'ncdf_meteo: Error allocating memory (wrk_dp)' if (err /= 0) stop 'ncdf_meteo: Error allocating memory (wrk_dp)'
wrk_dp = _ZERO_ wrk_dp = _ZERO_
if ( .not. on_grid ) then allocate(beta(E2DFIELD),stat=err)
if (err /= 0) &
stop 'init_meteo_input_ncdf: Error allocating memory (beta)'
beta = _ZERO_
allocate(beta(E2DFIELD),stat=err) if(iextr .eq. 1 .and. jextr .eq. 1) then
if (err /= 0) & point_source = .true.
stop 'init_meteo_input_ncdf: Error allocating memory (beta)' LEVEL3 'Assuming Point Source meteo forcing'
beta = _ZERO_ if ( .not. on_grid) then
LEVEL3 'Setting on_grid to true'
on_grid=.true.
end if
if (rotated_meteo_grid) then
rlon=met_lon(1)
rlat=met_lat(1)
call to_rotated_lat_lon(southpole,olon,olat,rlon,rlat,x)
beta = x
end if
end if
if ( .not. on_grid ) then
allocate(ti(E2DFIELD),stat=err) allocate(ti(E2DFIELD),stat=err)
if (err /= 0) & if (err /= 0) &
...@@ -240,6 +251,7 @@ ...@@ -240,6 +251,7 @@
if (err /= 0) stop & if (err /= 0) stop &
'init_meteo_input_ncdf: Error allocating memory (gridmap)' 'init_meteo_input_ncdf: Error allocating memory (gridmap)'
gridmap(:,:,:) = -999 gridmap(:,:,:) = -999
#if 0
#ifdef MED_15X15MINS_TEST #ifdef MED_15X15MINS_TEST
do i=1,iextr do i=1,iextr
met_lon(i) = -10.125 + (i-1)*1.125 met_lon(i) = -10.125 + (i-1)*1.125
...@@ -258,6 +270,8 @@ ...@@ -258,6 +270,8 @@
end do end do
grid_scan=0 grid_scan=0
#endif #endif
#endif
call init_grid_interpol(imin,imax,jmin,jmax,az, & call init_grid_interpol(imin,imax,jmin,jmax,az, &
lonc,latc,met_lon,met_lat,southpole,gridmap,beta,ti,ui) lonc,latc,met_lon,met_lat,southpole,gridmap,beta,ti,ui)
...@@ -397,8 +411,10 @@ ...@@ -397,8 +411,10 @@
t = loop*timestep t = loop*timestep
do indx=save_n,tmax do indx=save_n,tmax
if (met_times(indx) .gt. real(t + offset)) EXIT if (met_times(indx) .ge. real(t + offset)) EXIT
end do end do
! end of simulation?
if (indx .gt. tmax) then if (indx .gt. tmax) then
LEVEL2 'Need new meteo file' LEVEL2 'Need new meteo file'
call open_meteo_file(meteo_file) call open_meteo_file(meteo_file)
...@@ -406,7 +422,6 @@ ...@@ -406,7 +422,6 @@
do indx=save_n,tmax do indx=save_n,tmax
if (met_times(indx) .gt. real(t + offset)) EXIT if (met_times(indx) .gt. real(t + offset)) EXIT
end do end do
! indx = 1
save_n = 0 save_n = 0
end if end if
start(3) = indx start(3) = indx
...@@ -574,6 +589,7 @@ ...@@ -574,6 +589,7 @@
else else
err = nf_get_var_double(ncid,id,southpole) err = nf_get_var_double(ncid,id,southpole)
if (err .ne. NF_NOERR) go to 10 if (err .ne. NF_NOERR) go to 10
rotated_meteo_grid = .true.
end if end if
LEVEL4 'south pole:' LEVEL4 'south pole:'
LEVEL4 ' lon ',southpole(1) LEVEL4 ' lon ',southpole(1)
...@@ -702,7 +718,7 @@ ...@@ -702,7 +718,7 @@
! !LOCAL VARIABLES: ! !LOCAL VARIABLES:
integer :: i1,i2,istr,j1,j2,jstr integer :: i1,i2,istr,j1,j2,jstr
integer :: i,j,err integer :: i,j,err
REALTYPE :: uu,vv,sinconv,cosconv REALTYPE :: angle,uu,vv,sinconv,cosconv
!EOP !EOP
!----------------------------------------------------------------------- !-----------------------------------------------------------------------
include "netcdf.inc" include "netcdf.inc"
...@@ -712,11 +728,15 @@ ...@@ -712,11 +728,15 @@
err = nf_get_vara_real(ncid,u10_id,start,edges,wrk) err = nf_get_vara_real(ncid,u10_id,start,edges,wrk)
if (err .ne. NF_NOERR) go to 10 if (err .ne. NF_NOERR) go to 10
if (on_grid) then if (on_grid) then
do j=jmin,jmax if (point_source) then
do i=imin,imax u10 = wrk(1,1)
u10(i,j) = wrk(i,j) else
do j=jmin,jmax
do i=imin,imax
u10(i,j) = wrk(i,j)
end do
end do end do
end do end if
else else
!KBKwrk_dp = _ZERO_ !KBKwrk_dp = _ZERO_
call copy_var(grid_scan,wrk,wrk_dp) call copy_var(grid_scan,wrk,wrk_dp)
...@@ -726,43 +746,32 @@ ...@@ -726,43 +746,32 @@
err = nf_get_vara_real(ncid,v10_id,start,edges,wrk) err = nf_get_vara_real(ncid,v10_id,start,edges,wrk)
if (err .ne. NF_NOERR) go to 10 if (err .ne. NF_NOERR) go to 10
if (on_grid) then if (on_grid) then
do j=jmin,jmax if (point_source) then
do i=imin,imax v10 = wrk(1,1)
v10(i,j) = wrk(i,j) else
do j=jmin,jmax
do i=imin,imax
v10(i,j) = wrk(i,j)
end do
end do end do
end do end if
else else
!KBKwrk_dp = _ZERO_ !KBKwrk_dp = _ZERO_
call copy_var(grid_scan,wrk,wrk_dp) call copy_var(grid_scan,wrk,wrk_dp)
call do_grid_interpol(az,wrk_dp,gridmap,ti,ui,v10) call do_grid_interpol(az,wrk_dp,gridmap,ti,ui,v10)
end if end if
#ifdef SPHERICAL ! Rotation of wind due to the combined effect of possible rotation of
! Rotation of wind due to grid convergence ! meteorological grid and possible hydrodynamic grid convergence
do j=jmin,jmax ! (cartesian and curvi-linear grids where conv <> 0.)
do i=imin,imax
if(beta(i,j) .ne. _ZERO_) then
sinconv=sin(beta(i,j))
cosconv=cos(beta(i,j))
uu=u10(i,j)
vv=v10(i,j)
u10(i,j)= uu*cosconv+vv*sinconv
v10(i,j)=-uu*sinconv+vv*cosconv
end if
end do
end do
#else
if (southpole(1) .ne. 0. .or. southpole(2) .ne. -90.) then
LEVEL0 "rotation of wind due to the combined effect of grid convergence"
LEVEL0 "and rotated meteorological grid is not implemented yet."
stop "read_data()"
end if
! Rotation of wind due to grid convergence
do j=jmin,jmax do j=jmin,jmax
do i=imin,imax do i=imin,imax
if(conv(i,j) .ne. _ZERO_) then !KBK angle=-convc(i,j)*deg2rad
sinconv=sin(-conv(i,j)*deg2rad) !KBK angle=beta(i,j)
cosconv=cos(-conv(i,j)*deg2rad) angle=beta(i,j)-conv(i,j)*deg2rad
if(angle .ne. _ZERO_) then
sinconv=sin(angle)
cosconv=cos(angle)
uu=u10(i,j) uu=u10(i,j)
vv=v10(i,j) vv=v10(i,j)
u10(i,j)= uu*cosconv+vv*sinconv u10(i,j)= uu*cosconv+vv*sinconv
...@@ -770,16 +779,19 @@ LEVEL0 "and rotated meteorological grid is not implemented yet." ...@@ -770,16 +779,19 @@ LEVEL0 "and rotated meteorological grid is not implemented yet."
end if end if
end do end do
end do end do
#endif
err = nf_get_vara_real(ncid,airp_id,start,edges,wrk) err = nf_get_vara_real(ncid,airp_id,start,edges,wrk)
if (err .ne. NF_NOERR) go to 10 if (err .ne. NF_NOERR) go to 10
if (on_grid) then if (on_grid) then
do j=jmin,jmax if (point_source) then
do i=imin,imax airp = wrk(1,1)
airp(i,j) = wrk(i,j) else
do j=jmin,jmax
do i=imin,imax
airp(i,j) = wrk(i,j)
end do
end do end do
end do end if
else else
!KBKwrk_dp = _ZERO_ !KBKwrk_dp = _ZERO_
call copy_var(grid_scan,wrk,wrk_dp) call copy_var(grid_scan,wrk,wrk_dp)
...@@ -789,11 +801,15 @@ LEVEL0 "and rotated meteorological grid is not implemented yet." ...@@ -789,11 +801,15 @@ LEVEL0 "and rotated meteorological grid is not implemented yet."
err = nf_get_vara_real(ncid,t2_id,start,edges,wrk) err = nf_get_vara_real(ncid,t2_id,start,edges,wrk)
if (err .ne. NF_NOERR) go to 10 if (err .ne. NF_NOERR) go to 10
if (on_grid) then if (on_grid) then
do j=jmin,jmax if (point_source) then
do i=imin,imax t2 = wrk(1,1)
t2(i,j) = wrk(i,j) else
do j=jmin,jmax
do i=imin,imax
t2(i,j) = wrk(i,j)
end do
end do end do
end do end if
else else
!KBKwrk_dp = _ZERO_ !KBKwrk_dp = _ZERO_
call copy_var(grid_scan,wrk,wrk_dp) call copy_var(grid_scan,wrk,wrk_dp)
...@@ -803,11 +819,15 @@ LEVEL0 "and rotated meteorological grid is not implemented yet." ...@@ -803,11 +819,15 @@ LEVEL0 "and rotated meteorological grid is not implemented yet."
err = nf_get_vara_real(ncid,hum_id,start,edges,wrk) err = nf_get_vara_real(ncid,hum_id,start,edges,wrk)
if (err .ne. NF_NOERR) go to 10 if (err .ne. NF_NOERR) go to 10
if (on_grid) then if (on_grid) then
do j=jmin,jmax if (point_source) then
do i=imin,imax hum = wrk(1,1)
hum(i,j) = wrk(i,j) else
do j=jmin,jmax
do i=imin,imax
hum(i,j) = wrk(i,j)
end do
end do end do
end do end if
else else
!KBKwrk_dp = _ZERO_ !KBKwrk_dp = _ZERO_
call copy_var(grid_scan,wrk,wrk_dp) call copy_var(grid_scan,wrk,wrk_dp)
...@@ -817,11 +837,15 @@ LEVEL0 "and rotated meteorological grid is not implemented yet." ...@@ -817,11 +837,15 @@ LEVEL0 "and rotated meteorological grid is not implemented yet."
err = nf_get_vara_real(ncid,tcc_id,start,edges,wrk) err = nf_get_vara_real(ncid,tcc_id,start,edges,wrk)
if (err .ne. NF_NOERR) go to 10 if (err .ne. NF_NOERR) go to 10
if (on_grid) then if (on_grid) then
do j=jmin,jmax if (point_source) then
do i=imin,imax tcc = wrk(1,1)
tcc(i,j) = wrk(i,j) else
do j=jmin,jmax
do i=imin,imax
tcc(i,j) = wrk(i,j)
end do
end do end do
end do end if
else else
!KBKwrk_dp = _ZERO_ !KBKwrk_dp = _ZERO_
call copy_var(grid_scan,wrk,wrk_dp) call copy_var(grid_scan,wrk,wrk_dp)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment