Commit 4908e1b3 authored by Espen Sollum's avatar Espen Sollum
Browse files

Fixed a leap year problem in reading the fields

parent 04e2c9e4
......@@ -94,7 +94,7 @@ O_LEV_DBG = g # [0,g]
#LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper -lnetcdff
LIBS = -lgrib_api_f90 -lgrib_api -lm -ljasper $(NCOPT)
FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=large -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) $(NCOPT) $(FUSER) #-Warray-bounds -fcheck=all # -march=native
FFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV) -g -cpp -m64 -mcmodel=large -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -fbacktrace -O$(O_LEV) $(NCOPT) $(FUSER) #-Warray-bounds -fcheck=all # -march=native
DBGFLAGS = -I$(INCPATH1) -I$(INCPATH2) -O$(O_LEV_DBG) -g3 -ggdb3 -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
......
......@@ -280,8 +280,8 @@ module par_mod
! Max/min relative gain/loss of mass for NH3 reactions (1.0=100%)
!*******************************************************
! real, parameter :: max_nh3_gain = 1.e-10, max_nh3_loss=0.99
real, parameter :: max_nh3_gain = 0.99, max_nh3_loss=0.99
real, parameter :: max_nh3_gain = 1.e-10, max_nh3_loss=0.99
! real, parameter :: max_nh3_gain = 0.99, max_nh3_loss=0.99
!************************************
! Unit numbers for input/output files
......
......@@ -21,25 +21,25 @@
subroutine read_hourly_NH3field(itime)
!*****************************************************************************
! *
! Reads the NH3 field into memory *
! *
! AUTHOR: S. Eckhardt, 2019 *
! *
!*****************************************************************************
! *
! Variables: *
! *
! path(numpath) contains the path names *
! lonNH3(nxNH3) longitude of NH3 fields *
! latNH3(nyNH3) latitude of NH3 fields *
! altNH3(nzNH3) altitude of NH3 fields *
! etaNH3(nzNH3) eta-levels of NH3 fields *
! NH3_field(nxNH3,nyNH3,nzNH3) NH3 lifetime
! *
! *
!*****************************************************************************
!*****************************************************************************
! *
! Reads the NH3 field into memory *
! *
! AUTHOR: S. Eckhardt, 2019 *
! *
!*****************************************************************************
! *
! Variables: *
! *
! path(numpath) contains the path names *
! lonNH3(nxNH3) longitude of NH3 fields *
! latNH3(nyNH3) latitude of NH3 fields *
! altNH3(nzNH3) altitude of NH3 fields *
! etaNH3(nzNH3) eta-levels of NH3 fields *
! NH3_field(nxNH3,nyNH3,nzNH3) NH3 lifetime
! *
! *
!*****************************************************************************
use oh_mod
use par_mod
......@@ -60,159 +60,168 @@ subroutine read_hourly_NH3field(itime)
integer :: nid,ierr,xid,yid,zid,tid,vid,timeind
integer :: ix,jy,kz
integer :: jjjjmmdd, hhmmss, jjjj, mm, dd, hh
integer, save :: jjjj_save
real(kind=dp) :: jul,refday,flex_seconds,first_day_of_this_year,hours_this_year
nh3fields_path='/xnilu_wrk/users/ne/AMMONIA/LOSS/ECLv5-GFED-GEIA/NH3LOSS3D_144x143_b/'
! get the date in the simulation
!******************************************************
! get the date in the simulation
!******************************************************
jul=bdate+real(itime,kind=dp)/86400._dp
call caldate(jul,jjjjmmdd,hhmmss)
jul=bdate+real(itime,kind=dp)/86400._dp
call caldate(jul,jjjjmmdd,hhmmss)
! mm=(jjjjmmdd-(jjjjmmdd/10000)*10000)/100
! dd=(jjjjmmdd-(jjjjmmdd/100)*100)
! hh=hhmmss/10000
jjjj=jjjjmmdd/10000
write(jjjj_char,fmt='(i4)') jjjj
! dates in the netcdf files are given in seconds after 2007-1-1
!**************************************************************
refday=juldate(20070101,0)
flex_seconds=(jul-refday)*24.*3600. !convert days to seconds
first_day_of_this_year=juldate(jjjj*10000+101,0)
hours_this_year=(jul-first_day_of_this_year)*24
jjjj=jjjjmmdd/10000
write(jjjj_char,fmt='(i4)') jjjj
if (itime.eq.0) jjjj_save = jjjj
! dates in the netcdf files are given in seconds after 2007-1-1
!**************************************************************
refday=juldate(20070101,0)
flex_seconds=(jul-refday)*24.*3600. !convert days to seconds
first_day_of_this_year=juldate(jjjj*10000+101,0)
hours_this_year=(jul-first_day_of_this_year)*24
! write(*,*) 'first_day_of_this_year,hours_this_year'
! write(*,*) first_day_of_this_year,hours_this_year
! write(*,*) 'itime,jul,jjjjmmdd,hhmmss,jjjj,dd,mm,hh,refday,flex_seconds: '
! write(*,*) itime,jul,jjjjmmdd,hhmmss,jjjj,dd,mm,hh,refday,flex_seconds
! Read NH3 fields and level heights for the actual hour
!******************************************************
! Read NH3 fields and level heights for the actual hour
!******************************************************
! write(mm,fmt='(i2.2)') m
thefile=trim(nh3fields_path)//'NH3LOSS3D_144x143_'//jjjj_char//'_masl.nc'
ierr=nf90_open(trim(thefile),NF_NOWRITE,nid)
if(ierr.ne.0) then
write (*,*) 'The NH3 field at: '//thefile// ' could not be opened'
write (*,*) 'please copy the NH3 fields there, or change the path in the'
write (*,*) 'COMMAND namelist!'
write(*,*) nf_strerror(ierr)
stop
endif
if ((itime.eq.0).and.(lroot)) write(*,*) 'reading NH3 from: ', thefile
! inquire about variables
!ierr=nf_inq_dimid(nid,'longitude',xid)
ierr=nf_inq_dimid(nid,'lon',xid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
!ierr=nf_inq_dimid(nid,'latitude',yid)
ierr=nf_inq_dimid(nid,'lat',yid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimid(nid,'presnivs',zid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimid(nid,'time_counter',tid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
thefile=trim(nh3fields_path)//'NH3LOSS3D_144x143_'//jjjj_char//'_masl.nc'
ierr=nf90_open(trim(thefile),NF_NOWRITE,nid)
if(ierr.ne.0) then
write (*,*) 'The NH3 field at: '//thefile// ' could not be opened'
write (*,*) 'please copy the NH3 fields there, or change the path in the'
write (*,*) 'COMMAND namelist!'
write(*,*) nf_strerror(ierr)
stop
endif
if ((itime.eq.0.or.(jjjj.ne.jjjj_save)).and.(lroot)) write(*,*) 'reading NH3 from: ', thefile
! if (lroot) write(*,*) 'reading NH3 from: ', thefile
! if (lroot) write(*,*) 'jjjj_char: ', jjjj_char
! inquire about variables
!ierr=nf_inq_dimid(nid,'longitude',xid)
ierr=nf_inq_dimid(nid,'lon',xid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
!ierr=nf_inq_dimid(nid,'latitude',yid)
ierr=nf_inq_dimid(nid,'lat',yid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimid(nid,'presnivs',zid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimid(nid,'time_counter',tid)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
! if(m.eq.1) then
! read dimension sizes
ierr=nf_inq_dimlen(nid,xid,nxNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimlen(nid,yid,nyNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimlen(nid,zid,nzNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimlen(nid,tid,ntNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
if (itime.eq.0) then
if (lroot) write(*,*) 'Allocating x,y,z - ',nxNH3,nyNH3,nzNH3
! read dimension sizes
ierr=nf_inq_dimlen(nid,xid,nxNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimlen(nid,yid,nyNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimlen(nid,zid,nzNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_dimlen(nid,tid,ntNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
if (itime.eq.0) then
if (lroot) write(*,*) 'Allocating x,y,z - ',nxNH3,nyNH3,nzNH3
! allocate variables
allocate(lonNH3(nxNH3))
allocate(latNH3(nyNH3))
allocate(altNH3(nzNH3))
allocate(timeNH3(ntNH3))
allocate(NH3LOSS_field(nxNH3,nyNH3,nzNH3))
endif
! read dimension variables
!ierr=nf_inq_varid(nid,'longitude',xid)
ierr=nf_inq_varid(nid,'lon',xid)
ierr=nf_get_var_real(nid,xid,lonNH3)
!rite(*,*) lonNH3
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
allocate(lonNH3(nxNH3))
allocate(latNH3(nyNH3))
allocate(altNH3(nzNH3))
allocate(timeNH3(ntNH3))
allocate(NH3LOSS_field(nxNH3,nyNH3,nzNH3))
! take into account that time dimension may vary for different years
else if (jjjj.ne.jjjj_save) then
deallocate(timeNH3)
allocate(timeNH3(ntNH3))
jjjj_save = jjjj
endif
! read dimension variables
!ierr=nf_inq_varid(nid,'longitude',xid)
ierr=nf_inq_varid(nid,'lon',xid)
ierr=nf_get_var_real(nid,xid,lonNH3)
!rite(*,*) lonNH3
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
! ierr=nf_inq_varid(nid,'latitude',yid)
ierr=nf_inq_varid(nid,'lat',yid)
ierr=nf_get_var_real(nid,yid,latNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_varid(nid,'presnivs',zid)
ierr=nf_get_var_real(nid,zid,altNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_varid(nid,'time_counter',tid)
ierr=nf_get_var_double(nid,tid,timeNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_varid(nid,'vmr_rat_nh3',vid)
! read NH3_field, this reads only a subset from index start the N fields specified in count
! rona hourlyOH ijx=minloc(abs(lonjr-lonOH(ix)),dim=1,mask=abs(lonjr-lonOH(ix)).eq.minval(abs(lonjr-lonOH(ix))))
timeind=minloc(abs(timeNH3-flex_seconds),dim=1,mask=abs(timeNH3-flex_seconds).eq.minval(abs(timeNH3-flex_seconds)))
ierr=nf_inq_varid(nid,'lat',yid)
ierr=nf_get_var_real(nid,yid,latNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_varid(nid,'presnivs',zid)
ierr=nf_get_var_real(nid,zid,altNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_varid(nid,'time_counter',tid)
ierr=nf_get_var_double(nid,tid,timeNH3)
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_inq_varid(nid,'vmr_rat_nh3',vid)
! read NH3_field, this reads only a subset from index start the N fields specified in count
! rona hourlyOH ijx=minloc(abs(lonjr-lonOH(ix)),dim=1,mask=abs(lonjr-lonOH(ix)).eq.minval(abs(lonjr-lonOH(ix))))
timeind=minloc(abs(timeNH3-flex_seconds),dim=1,mask=abs(timeNH3-flex_seconds).eq.minval(abs(timeNH3-flex_seconds)))
! write(*,*) 'NH3 calculated timeind: ',timeind
ierr=nf90_get_var(nid,vid,NH3LOSS_field(:,:,:),start=(/1,1,1,timeind/),count=(/nxNH3,nyNH3,nzNH3,1/))
ierr=nf90_get_var(nid,vid,NH3LOSS_field(:,:,:),start=(/1,1,1,timeind/),count=(/nxNH3,nyNH3,nzNH3,1/))
! write(*,*) 'Lon: ',lonNH3(1),' to ',lonNH3(nxNH3),'lat: ',latNH3(1),' to ',latNH3(nyNH3)
! do jy=1,nyNH3
! write(*,*) (NH3LOSS_field(ix,jy,1),ix=1,nxNH3)
! enddo
! rona original line call check( nf90_get_var(ncid,varid,flx(:,:,1:cnt),start=(/ix,jy,kt/),count=(/nx,ny,cnt/)) )
! rona original line call check( nf90_get_var(ncid,varid,flx(:,:,1:cnt),start=(/ix,jy,kt/),count=(/nx,ny,cnt/)) )
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
if(ierr.ne.0) then
write(*,*) nf_strerror(ierr)
stop
endif
ierr=nf_close(nid)
ierr=nf_close(nid)
return
end subroutine read_hourly_NH3field
......@@ -186,8 +186,10 @@ subroutine timemanager(metdata_format)
call wetdepo(itime,lsynctime,loutnext)
endif
if (NH3LOSS) &
call read_hourly_NH3field(itime)
if (NH3LOSS) then
! write(*,*) 'timemanager> itime: ',itime
call read_hourly_NH3field(itime)
end if
if ((OHREA .or. NH3LOSS) .and. itime .ne. 0 .and. numpart .gt. 0) &
call ohreaction(itime,lsynctime,loutnext)
......
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