Skip to content

GitLab

Commit 7ff942a7 authored by Sabine's avatar Sabine
Browse files

implementation for nested metfields

parent 9c52f03c
Hide whitespace changes
Inline Side-by-side
src/calculate_watercycle.f90
View file @ 7ff942a7
......@@ -34,11 +34,11 @@ subroutine calculate_watercycle(partnumber,itime)
implicit none
real(kind=dp) :: jul
integer :: itime,i,j,jjjjmmdd,ihmmss,numshortout,numshortall
integer :: itime,i,j,k,jjjjmmdd,ihmmss,numshortout,numshortall
integer :: ix,jy,ixp,jyp,ind,indz,indzp,il,m,indexh,itage,nage
integer :: interp_time
integer :: interp_time, ngrid
integer :: partnumber,forparticle ! -1 inititalize q, partnumber real do deposition
real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4
real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,xtn,ytn
real centerposx_real, centerposy_real, centerposx, centerposy, diff, xl(2), yl(2)
real qv1(2),qvprof(2),qvi,dz,dz1,dz2
real e_minus_p1(2), e_minus_pi
......@@ -91,13 +91,39 @@ subroutine calculate_watercycle(partnumber,itime)
end do
33 continue
ix=xtra1(i)
jy=ytra1(i)
! Determine which nesting level to be used
!*****************************************
ngrid=0
do k=numbnests,1,-1
if ((xtra1(i).gt.xln(k)).and.(xtra1(i).lt.xrn(k)).and. &
(ytra1(i).gt.yln(k)).and.(ytra1(i).lt.yrn(k))) then
ngrid=k
goto 23
endif
end do
23 continue
! Determine nested grid coordinates
!**********************************
if (ngrid.gt.0) then
xtn=(xtra1(i)-xln(ngrid))*xresoln(ngrid)
ytn=(ytra1(i)-yln(ngrid))*yresoln(ngrid)
ix=int(xtn)
jy=int(ytn)
ddy=ytn-real(jy)
ddx=xtn-real(ix)
else
jy=int(ytra1(i))
ddx=xtra1(i)-real(ix)
ddy=ytra1(i)-real(jy)
endif
ixp=ix+1
jyp=jy+1
if (jy.eq.numygrid) jyp=numygrid
ddx=xtra1(i)-real(ix)
ddy=ytra1(i)-real(jy)
rddx=1.-ddx
rddy=1.-ddy
p1=rddx*rddy
......@@ -126,11 +152,17 @@ subroutine calculate_watercycle(partnumber,itime)
do m=1,2
indexh=memind(m)
qv1(m)=p1*qv(ix ,jy ,ind,indexh) &
if (ngrid.eq.0) then
qv1(m)=p1*qv(ix ,jy ,ind,indexh) &
+p2*qv(ixp,jy ,ind,indexh) &
+p3*qv(ix ,jyp,ind,indexh) &
+p4*qv(ixp,jyp,ind,indexh)
else
qv1(m)=p1*qvn(ix ,jy ,ind,indexh,ngrid) &
+p2*qvn(ixp,jy ,ind,indexh,ngrid) &
+p3*qvn(ix ,jyp,ind,indexh,ngrid) &
+p4*qvn(ixp,jyp,ind,indexh,ngrid)
endif
end do
qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt
......@@ -140,13 +172,20 @@ subroutine calculate_watercycle(partnumber,itime)
do m=1,2
indexh=memind(m)
e_minus_p1(m)=p1*e_minus_p(ix ,jy ,indexh) &
if (ngrid.eq.0) then
e_minus_p1(m)=p1*e_minus_p(ix ,jy ,indexh) &
+p2*e_minus_p(ixp,jy ,indexh) &
+p3*e_minus_p(ix ,jyp,indexh) &
+p4*e_minus_p(ixp,jyp,indexh)
else
e_minus_p1(m)=p1*e_minus_p_nest(ix ,jy ,indexh,ngrid) &
+p2*e_minus_p_nest(ixp,jy ,indexh,ngrid) &
+p3*e_minus_p_nest(ix ,jyp,indexh,ngrid) &
+p4*e_minus_p_nest(ixp,jyp,indexh,ngrid)
endif
end do
e_minus_pi=(e_minus_p1(1)*dt2+dt1*e_minus_p1(2))*dtt
! e_minus_pi=e_minus_p1(memtime(1)) !test how the results change without interpolation - does not change a lot
! write(*,*) 'all: ',i,diff,e_minus_p(ix,jy) &
......@@ -173,7 +212,7 @@ subroutine calculate_watercycle(partnumber,itime)
! we cannot be sure 3600 is the second timestep, could also be 1800!
if ( ((diff.ge.0).or.((e_minus_pi/24).gt.-2.0)) .and. (itage.le.3600)) then
if (((diff.ge.0).or.((e_minus_pi)/24.gt.-2.0)) .and. (itage.le.3600)) then
! write(*,*) 'terminated: ',i,diff,e_minus_p(ix,jy) &
! ,partnumber,val_q(i),qvi,itime,itra1(i),itage
itra1(i)=-999999999
......@@ -185,8 +224,12 @@ subroutine calculate_watercycle(partnumber,itime)
centerposy=(centerposy_real-ylat0)/dy
if (diff.gt.0) then
call drydepokernel(1,diff,centerposx,centerposy,nage,1)
if (nested_output.eq.1) &
call drydepokernel_nest(1,diff,centerposx,centerposy,nage,1)
else
call wetdepokernel(1,abs(diff),centerposx,centerposy,nage,1)
if (nested_output.eq.1) &
call wetdepokernel_nest(1,abs(diff),centerposx,centerposy,nage,1)
endif
endif
endif
......
src/com_mod.f90
View file @ 7ff942a7
......@@ -496,6 +496,7 @@ module com_mod
& rhon, drhodzn, tthn, qvhn, clwcn, ciwcn, clwn, clwchn, ciwchn
real,allocatable,dimension(:,:,:,:) :: ctwcn
real,allocatable, dimension (:,:,:,:) :: uul,vvl
real,allocatable, dimension (:,:,:,:,:) :: uuln,vvln
integer,allocatable,dimension(:,:,:,:) :: cloudshn
integer(kind=1),allocatable,dimension(:,:,:,:,:) :: cloudsn
......@@ -679,6 +680,7 @@ module com_mod
real, allocatable, dimension(:,:) :: xscav_frac1
real, allocatable, dimension(:,:,:) :: e_minus_p ! E-P field for the watercycle
real, allocatable, dimension(:,:,:,:) :: e_minus_p_nest ! E-P field for the watercycle
real, allocatable, dimension(:) :: xtra1_q, ytra1_q, val_q
! eso: Moved from timemanager
......
src/euler_rain.f90
View file @ 7ff942a7
......@@ -26,10 +26,8 @@
real watercolumn(2),watermass(nzmax,2)
real dwatercolumndt(0:nxmax-1,0:nymax-1)
real divwatercolumn(0:nxmax-1,0:nymax-1,3)
integer nmax
parameter(nmax=nxmax*nymax)
write(89,*) xlon0+dx,ylat0+dy,nx-2,ny-2,dx,dy,1,10000.,3
! write(89,*) xlon0+dx,ylat0+dy,nx-2,ny-2,dx,dy,1,10000.,3
! Loop over the whole grid
!*************************
......@@ -157,7 +155,7 @@
! Set values at the poles equal to values 1 deg equatorward
!**********************************************************
do n=1,2
n=memind(1)
do 80 ix=1,nx-1
e_minus_p(ix,ny-1,n)=e_minus_p(ix,ny-2,n)
......@@ -169,5 +167,4 @@
write(89,*) (e_minus_p(ix,jy,n),ix=1,nx-2)
85 continue
end do
end
src/getfields.f90
View file @ 7ff942a7
......@@ -148,6 +148,7 @@ subroutine getfields(itime,nstop,metdata_format)
memtime(2)=wftime(indj+1)
if (WATERCYCLE) &
call euler_rain(itime) ! calculate the integrated E-P
call euler_rain_nest(itime) ! calculate the integrated E-P
nstop = 1
goto 40
endif
......@@ -184,8 +185,10 @@ subroutine getfields(itime,nstop,metdata_format)
call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
memtime(1)=wftime(indj)
memind(2)=2
if (WATERCYCLE) &
if (WATERCYCLE) then
call euler_rain(itime) ! calculate the integrated E-P
call euler_rain_nest(itime) ! calculate the integrated E-P
endif
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
else
......@@ -201,8 +204,10 @@ subroutine getfields(itime,nstop,metdata_format)
end if
call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
memtime(2)=wftime(indj+1)
if (WATERCYCLE) &
if (WATERCYCLE) then
call euler_rain(itime) ! calculate the integrated E-P
call euler_rain_nest(itime) ! calculate the integrated E-P
endif
nstop = 1
goto 60
endif
......
src/makefile
View file @ 7ff942a7
......@@ -90,11 +90,11 @@ O_LEV_DBG = g # [0,g]
## LIBRARIES
LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff # -fopenmp
FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) # -Warray-bounds -fcheck=all # -march=native
#FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) -Warray-bounds -fcheck=all # -march=native
#FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) # -Warray-bounds -fcheck=all # -march=native
FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(FUSER) -Warray-bounds -fcheck=all # -march=native
DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace -Wall -fdump-core $(FUSER) # -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all
#DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace -Wall -fdump-core $(FUSER) -Warray-bounds -fcheck=all
#DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace -Wall -fdump-core $(FUSER) # -ffpe-trap=invalid,overflow,denormal,underflow,zero -Warray-bounds -fcheck=all
DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV_DBG) -fbacktrace -Wall -fdump-core $(FUSER) -Warray-bounds -fcheck=all
LDFLAGS = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
LDDEBUG = $(DBGFLAGS) -L$(LIBPATH1) $(LIBS) #-L$(LIBPATH2)
......@@ -126,6 +126,7 @@ OBJECTS_SERIAL = \
concoutput_surf.o concoutput_surf_nest.o \
getfields.o \
euler_rain.o \
euler_rain_nests.o \
readwind_ecmwf.o
## For MPI version
......@@ -291,6 +292,7 @@ fluxoutput.o: com_mod.o flux_mod.o outg_mod.o par_mod.o
get_settling.o: com_mod.o par_mod.o
getfields.o: com_mod.o par_mod.o class_gribfile_mod.o
euler_rain.o: com_mod.o par_mod.o
euler_rain_nests.o: com_mod.o par_mod.o
getfields_mpi.o: com_mod.o par_mod.o mpi_mod.o class_gribfile_mod.o
gethourlyOH.o: com_mod.o oh_mod.o par_mod.o
getrb.o: par_mod.o
......
src/outgrid_init.f90
View file @ 7ff942a7
......@@ -221,6 +221,8 @@ subroutine outgrid_init
if (WATERCYCLE) then
allocate(uul(0:nxmax-1,0:nymax-1,nuvzmax,2))
allocate(vvl(0:nxmax-1,0:nymax-1,nuvzmax,2))
allocate(uuln(0:nxmaxn-1,0:nymaxn-1,nuvzmax,2,maxnests))
allocate(vvln(0:nxmaxn-1,0:nymaxn-1,nuvzmax,2,maxnests))
endif
! Extra field for totals at MPI root process
......@@ -300,6 +302,7 @@ subroutine outgrid_init
if (stat.ne.0) write(*,*)'ERROR: could not allocate gridunc'
if (WATERCYCLE) then
allocate(e_minus_p(0:nxmax-1,0:nymax-1,2))
allocate(e_minus_p_nest(0:nxmaxn-1,0:nymaxn-1,2,maxnests))
endif
endif
......
src/par_mod.f90
View file @ 7ff942a7
......@@ -159,7 +159,7 @@ module par_mod
! Maximum dimensions of the nested input grids
!*********************************************
integer,parameter :: maxnests=0,nxmaxn=451,nymaxn=226
integer,parameter :: maxnests=1,nxmaxn=351,nymaxn=151
! nxmax,nymax maximum dimension of wind fields in x and y
! direction, respectively
......
src/readwind_ecmwf.f90
View file @ 7ff942a7
......@@ -487,16 +487,16 @@ subroutine readwind_ecmwf(indj,n,uuh,vvh,wwh)
do i=0,nxmin1
do j=0,nymin1
uuh(i,j,1)=u10(i,j,1,n)
vvh(i,j,1)=v10(i,j,1,n)
qvh(i,j,1,n)=qvh(i,j,2,n)
tth(i,j,1,n)=tt2(i,j,1,n)
if (WATERCYCLE) then
do ki=1,nuvz
uul(i,j,ki,n)=uuh(i,j,ki)
vvl(i,j,ki,n)=vvh(i,j,ki)
end do
endif
uuh(i,j,1)=u10(i,j,1,n)
vvh(i,j,1)=v10(i,j,1,n)
qvh(i,j,1,n)=qvh(i,j,2,n)
tth(i,j,1,n)=tt2(i,j,1,n)
end do
end do
......
src/readwind_nests.f90
View file @ 7ff942a7
......@@ -57,7 +57,7 @@ subroutine readwind_nests(indj,n,uuhn,vvhn,wwhn)
real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
integer :: indj,i,j,k,n,levdiff2,ifield,iumax,iwmax,l
integer :: indj,i,j,k,ki,n,levdiff2,ifield,iumax,iwmax,l
! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
......@@ -436,6 +436,12 @@ subroutine readwind_nests(indj,n,uuhn,vvhn,wwhn)
vvhn(i,j,1,l)=v10n(i,j,1,n,l)
qvhn(i,j,1,n,l)=qvhn(i,j,2,n,l)
tthn(i,j,1,n,l)=tt2n(i,j,1,n,l)
if (WATERCYCLE) then
do ki=1,nuvz
uuln(i,j,ki,n,l)=uuhn(i,j,ki,l)
vvln(i,j,ki,n,l)=vvhn(i,j,ki,l)
end do
endif
end do
end do
......
src/timemanager.f90
View file @ 7ff942a7
......@@ -127,6 +127,7 @@ subroutine timemanager(metdata_format)
open(89,file=path(2)(1:length(2))//'budget.dat')
open(90,file=path(2)(1:length(2))//'budget_nest.dat')
! First output for time 0
!************************
......@@ -745,6 +746,7 @@ subroutine timemanager(metdata_format)
if (WATERCYCLE) then
close(89) ! Euler rain, water budget file
close(90) ! Euler rain, water budget file
endif
! De-allocate memory and end
......@@ -772,7 +774,7 @@ subroutine timemanager(metdata_format)
deallocate(outheight,outheighthalf)
deallocate(oroout, area, volume)
if (WATERCYCLE) then
deallocate(uul,vvl)
deallocate(uul,vvl,uuln,vvln)
endif
end subroutine timemanager
Supports Markdown
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