Skip to content

Commit

Permalink
fv3atm issue NOAA-EMC#37: fix the real(8) lat/lon in netcdf file
Browse files Browse the repository at this point in the history
  • Loading branch information
@junwang-noaa
junwang-noaa committed Jan 7, 2020
1 parent 869374c commit 2234939
Show file tree
Hide file tree
Showing 3 changed files with 48 additions and 29 deletions.
2 changes: 1 addition & 1 deletion io/module_write_nemsio.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ subroutine nemsio_first_call(fieldbundle, imo, jmo, &
integer, intent(in) :: wrt_mype, wrt_ntasks, wrt_mpi_comm
integer, intent(in) :: wrt_nbdl, mybdl
integer, intent(in) :: inidate(7)
real, intent(in) :: lat(:), lon(:)
real(8), intent(in) :: lat(:), lon(:)
integer, optional,intent(out) :: rc

!** local vars
Expand Down
11 changes: 5 additions & 6 deletions io/module_write_netcdf.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
real(4), dimension(:,:,:), allocatable :: arrayr4_3d,arrayr4_3d_save
real(8), dimension(:,:,:), allocatable :: arrayr8_3d

real(8) rad2dg,x(im),y(jm)
real(8) x(im),y(jm)
integer :: fieldCount, fieldDimCount, gridDimCount
integer, dimension(:), allocatable :: ungriddedLBound, ungriddedUBound

Expand Down Expand Up @@ -71,7 +71,6 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
!
!!
!
rad2dg = 45./atan(1.0)

call ESMF_FieldBundleGet(fieldbundle, fieldCount=fieldCount, rc=rc); ESMF_ERR_RETURN(rc)

Expand Down Expand Up @@ -247,7 +246,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
if (mype==0) then
if (trim(output_grid) == 'gaussian_grid' .or. &
trim(output_grid) == 'regional_latlon') then
ncerr = nf90_put_var(ncid, im_varid, values=rad2dg*arrayr8(:,1) ); NC_ERR_STOP(ncerr)
ncerr = nf90_put_var(ncid, im_varid, values=arrayr8(:,1) ); NC_ERR_STOP(ncerr)
ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, im_varid, "long_name", "T-cell longitude"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, im_varid, "units", "degrees_E"); NC_ERR_STOP(ncerr)
Expand All @@ -271,15 +270,15 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
ncerr = nf90_put_att(ncid, im_varid, "units", "meters"); NC_ERR_STOP(ncerr)
ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
endif
ncerr = nf90_put_var(ncid, lon_varid, values=rad2dg*arrayr8 ); NC_ERR_STOP(ncerr)
ncerr = nf90_put_var(ncid, lon_varid, values=arrayr8 ); NC_ERR_STOP(ncerr)
endif

call ESMF_GridGetCoord(wrtGrid, coordDim=2, array=array, rc=rc); ESMF_ERR_RETURN(rc)
call ESMF_ArrayGather(array, arrayr8, rootPet=0, rc=rc); ESMF_ERR_RETURN(rc)
if (mype==0) then
if (trim(output_grid) == 'gaussian_grid' .or. &
trim(output_grid) == 'regional_latlon') then
ncerr = nf90_put_var(ncid, jm_varid, values=rad2dg*arrayr8(1,:) ); NC_ERR_STOP(ncerr)
ncerr = nf90_put_var(ncid, jm_varid, values=arrayr8(1,:) ); NC_ERR_STOP(ncerr)
ncerr = nf90_redef(ncid=ncid); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "long_name", "T-cell latiitude"); NC_ERR_STOP(ncerr)
ncerr = nf90_put_att(ncid, jm_varid, "units", "degrees_N"); NC_ERR_STOP(ncerr)
Expand All @@ -303,7 +302,7 @@ subroutine write_netcdf(fieldbundle, wrtfb, filename, mpi_comm, mype, im, jm, rc
ncerr = nf90_put_att(ncid, jm_varid, "units", "meters"); NC_ERR_STOP(ncerr)
ncerr = nf90_enddef(ncid=ncid); NC_ERR_STOP(ncerr)
endif
ncerr = nf90_put_var(ncid, lat_varid, values=rad2dg*arrayr8 ); NC_ERR_STOP(ncerr)
ncerr = nf90_put_var(ncid, lat_varid, values=arrayr8 ); NC_ERR_STOP(ncerr)
endif

do i=1, fieldCount
Expand Down
64 changes: 42 additions & 22 deletions io/module_wrt_grid_comp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -183,7 +183,7 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc)
character(128) :: FBlist_outfilename(100), outfile_name
character(128),dimension(:,:), allocatable :: outfilename
real(8), dimension(:), allocatable :: slat
real, dimension(:), allocatable :: lat, lon, axesdata
real(8), dimension(:), allocatable :: lat, lon
real(ESMF_KIND_R8), dimension(:,:), pointer :: lonPtr, latPtr
real(ESMF_KIND_R8) :: rot_lon, rot_lat
real(ESMF_KIND_R8) :: geo_lon, geo_lat
Expand Down Expand Up @@ -358,19 +358,20 @@ subroutine wrt_initialize(wrt_comp, imp_state_write, exp_state_write, clock, rc)
wrt_int_state%latstart = lat(1)
wrt_int_state%latlast = lat(jmo)
do j=1,imo
lon(j) = 360./real(imo) *real(j-1)
lon(j) = 360.d0/real(imo,8) *real(j-1,8)
enddo
wrt_int_state%lonstart = lon(1)
wrt_int_state%lonlast = lon(imo)
do j=lbound(latPtr,2),ubound(latPtr,2)
do i=lbound(lonPtr,1),ubound(lonPtr,1)
lonPtr(i,j) = 360./real(imo) * (i-1)
lonPtr(i,j) = 360.d0/real(imo,8) * real(i-1,8)
latPtr(i,j) = lat(j)
enddo
enddo
! print *,'aft wrtgrd, Gaussian, dimi,i=',lbound(lonPtr,1),ubound(lonPtr,1), &
! ' j=',lbound(lonPtr,2),ubound(lonPtr,2),'imo=',imo,'jmo=',jmo
! print *,'aft wrtgrd, lon=',lonPtr(lbound(lonPtr,1),lbound(lonPtr,2)), &
! if(wrt_int_state%mype==0) print *,'aft wrtgrd, lon=',lonPtr(1:5,1), &
! 'lat=',latPtr(1,1:5),'imo,jmo=',imo,jmo
! lonPtr(lbound(lonPtr,1),ubound(lonPtr,2)),'lat=',latPtr(lbound(lonPtr,1),lbound(lonPtr,2)), &
! latPtr(lbound(lonPtr,1),ubound(lonPtr,2))
wrt_int_state%lat_start = lbound(latPtr,2)
Expand Down Expand Up @@ -1622,13 +1623,14 @@ subroutine recover_fields(file_bundle,rc)
character(100) fieldName,uwindname,vwindname
type(ESMF_Field), allocatable :: fcstField(:)
real(ESMF_KIND_R8), dimension(:,:), pointer :: lon, lat
real(ESMF_KIND_R8), dimension(:,:), pointer :: lonloc, latloc
real(ESMF_KIND_R4), dimension(:,:), pointer :: pressfc
real(ESMF_KIND_R4), dimension(:,:), pointer :: uwind2dr4,vwind2dr4
real(ESMF_KIND_R4), dimension(:,:,:), pointer :: uwind3dr4,vwind3dr4
real(ESMF_KIND_R4), dimension(:,:,:), pointer :: cart3dPtr2dr4
real(ESMF_KIND_R4), dimension(:,:,:,:), pointer :: cart3dPtr3dr4
real(ESMF_KIND_R8), dimension(:,:,:,:), pointer :: cart3dPtr3dr8
save lon, lat
save lonloc, latloc
real(ESMF_KIND_R8) :: coslon, sinlon, sinlat
!
! get filed count
Expand All @@ -1648,19 +1650,37 @@ subroutine recover_fields(file_bundle,rc)

if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return

lon = lon * pi/180.
! print *,'in 3DCartesian2wind, lon dim=',lbound(lon,1),ubound(lon,1),lbound(lon,2),ubound(lon,2), &
! 'lon=',lon(lbound(lon,1),lbound(lon,2)), lon(ubound(lon,1),ubound(lon,2))
allocate(lonloc(lbound(lon,1):ubound(lon,1),lbound(lon,2):ubound(lon,2)))
istart = lbound(lon,1)
iend = ubound(lon,1)
jstart = lbound(lon,2)
jend = ubound(lon,2)
!$omp parallel do default(none) shared(lon,lonloc,jstart,jend,istart,iend) &
!$omp private(i,j)
do j=jstart,jend
do i=istart,iend
lonloc(i,j) = lon(i,j) * pi/180.
enddo
enddo

CALL ESMF_LogWrite("call recover field get coord 2",ESMF_LOGMSG_INFO,rc=RC)

call ESMF_GridGetCoord(fieldgrid, coordDim=2, farrayPtr=lat, rc=rc)

if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, line=__LINE__, file=__FILE__)) return

lat = lat * pi/180.
! print *,'in 3DCartesian2wind, lat dim=',lbound(lat,1),ubound(lat,1),lbound(lat,2),ubound(lat,2), &
! 'lat=',lat(lbound(lon,1),lbound(lon,2)), lat(ubound(lon,1),ubound(lon,2))
allocate(latloc(lbound(lat,1):ubound(lat,1),lbound(lat,2):ubound(lat,2)))
istart = lbound(lat,1)
iend = ubound(lat,1)
jstart = lbound(lat,2)
jend = ubound(lat,2)
!$omp parallel do default(none) shared(lat,latloc,jstart,jend,istart,iend) &
!$omp private(i,j)
do j=jstart,jend
do i=istart,iend
latloc(i,j) = lat(i,j) * pi/180.d0
enddo
enddo
first_getlatlon = .false.
endif
!
Expand Down Expand Up @@ -1718,18 +1738,18 @@ subroutine recover_fields(file_bundle,rc)
! update u , v wind
!$omp parallel do default(shared) private(i,j,k,coslon,sinlon,sinlat)
do k=kstart,kend
!!$omp parallel do default(none) shared(uwind3dr4,vwind3dr4,lon,lat,cart3dPtr3dr4,jstart,jend,istart,iend,k) &
!!$omp private(i,j,coslon,sinlon,sinlat)
!$omp parallel do default(none) shared(uwind3dr4,vwind3dr4,lonloc,latloc,cart3dPtr3dr4,jstart,jend,istart,iend,k) &
!$omp private(i,j,coslon,sinlon,sinlat)
do j=jstart, jend
do i=istart, iend
coslon = cos(lon(i,j))
sinlon = sin(lon(i,j))
sinlat = sin(lat(i,j))
coslon = cos(lonloc(i,j))
sinlon = sin(lonloc(i,j))
sinlat = sin(latloc(i,j))
uwind3dr4(i,j,k) = cart3dPtr3dr4(1,i,j,k) * coslon &
+ cart3dPtr3dr4(2,i,j,k) * sinlon
vwind3dr4(i,j,k) =-cart3dPtr3dr4(1,i,j,k) * sinlat*sinlon &
+ cart3dPtr3dr4(2,i,j,k) * sinlat*coslon &
+ cart3dPtr3dr4(3,i,j,k) * cos(lat(i,j))
+ cart3dPtr3dr4(3,i,j,k) * cos(latloc(i,j))
enddo
enddo
enddo
Expand All @@ -1749,18 +1769,18 @@ subroutine recover_fields(file_bundle,rc)
call ESMF_FieldGet(ufield, localDe=0, farrayPtr=uwind2dr4,rc=rc)
call ESMF_FieldGet(vfield, localDe=0, farrayPtr=vwind2dr4,rc=rc)
! update u , v wind
!$omp parallel do default(none) shared(uwind2dr4,vwind2dr4,lon,lat,cart3dPtr2dr4,jstart,jend,istart,iend) &
!$omp parallel do default(none) shared(uwind2dr4,vwind2dr4,lonloc,latloc,cart3dPtr2dr4,jstart,jend,istart,iend) &
!$omp private(i,j,k,coslon,sinlon,sinlat)
do j=jstart, jend
do i=istart, iend
coslon = cos(lon(i,j))
sinlon = sin(lon(i,j))
sinlat = sin(lat(i,j))
coslon = cos(lonloc(i,j))
sinlon = sin(lonloc(i,j))
sinlat = sin(latloc(i,j))
uwind2dr4(i,j) = cart3dPtr2dr4(1,i,j) * coslon &
+ cart3dPtr2dr4(2,i,j) * sinlon
vwind2dr4(i,j) =-cart3dPtr2dr4(1,i,j) * sinlat*sinlon &
+ cart3dPtr2dr4(2,i,j) * sinlat*coslon &
+ cart3dPtr2dr4(3,i,j) * cos(lat(i,j))
+ cart3dPtr2dr4(3,i,j) * cos(latloc(i,j))
enddo
enddo
endif
Expand Down

0 comments on commit 2234939

Please sign in to comment.