From 2234939b31a338d3fca55a8c9e992cdfae89c219 Mon Sep 17 00:00:00 2001 From: "Jun.Wang" Date: Tue, 7 Jan 2020 21:31:06 +0000 Subject: [PATCH] fv3atm issue #37: fix the real(8) lat/lon in netcdf file --- io/module_write_nemsio.F90 | 2 +- io/module_write_netcdf.F90 | 11 +++---- io/module_wrt_grid_comp.F90 | 64 ++++++++++++++++++++++++------------- 3 files changed, 48 insertions(+), 29 deletions(-) diff --git a/io/module_write_nemsio.F90 b/io/module_write_nemsio.F90 index 3afd66789..e51f64a52 100644 --- a/io/module_write_nemsio.F90 +++ b/io/module_write_nemsio.F90 @@ -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 diff --git a/io/module_write_netcdf.F90 b/io/module_write_netcdf.F90 index 1fce3d8b9..52cf740e6 100644 --- a/io/module_write_netcdf.F90 +++ b/io/module_write_netcdf.F90 @@ -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 @@ -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) @@ -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) @@ -271,7 +270,7 @@ 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) @@ -279,7 +278,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, 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) @@ -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 diff --git a/io/module_wrt_grid_comp.F90 b/io/module_wrt_grid_comp.F90 index aad3991d0..84769eaea 100644 --- a/io/module_wrt_grid_comp.F90 +++ b/io/module_wrt_grid_comp.F90 @@ -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 @@ -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) @@ -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 @@ -1648,9 +1650,18 @@ 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) @@ -1658,9 +1669,18 @@ subroutine recover_fields(file_bundle,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 ! @@ -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 @@ -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