Skip to content

Commit

Permalink
Merge pull request #1152 from LIS-navari/feature/lambert_4_lvt_DAobs
Browse files Browse the repository at this point in the history 
Support for Lambert projection added to LVT LISDAobs reader.
  • Loading branch information
@cmclaug2
cmclaug2 authored Oct 11, 2022
2 parents a7b4416 + 76eda29 commit 9d9c5fe
Showing 1 changed file with 159 additions and 84 deletions.
243 changes: 159 additions & 84 deletions lvt/datastreams/LISDAobs/LISda_obsMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@ module LISda_obsMod
!
! !REVISION HISTORY:
! 02 Oct 2008 Sujay Kumar Initial Specification
! July 18 2022 Madi Navari Added support for Lambert projection.
!
!EOP

Expand Down Expand Up @@ -113,6 +114,9 @@ subroutine LISda_obsInit(i)
logical :: file_exists
character*20 :: gridtype
integer :: rc
integer :: truelat1, truelat2
integer :: domain_resolution
integer :: standard_lon

if(.not.allocated(lisdaobs)) then
allocate(lisdaobs(LVT_rc%nDataStreams))
Expand Down Expand Up @@ -168,105 +172,176 @@ subroutine LISda_obsInit(i)
if(file_exists) then
ios = nf90_open(path=domFile,mode=NF90_NOWRITE,ncid=ftn)
call LVT_verify(ios,'Error in nf90_open in readObsDomainInput')

ios = nf90_get_att(ftn, NF90_GLOBAL, 'MAP_PROJECTION',map_proj)
call LVT_verify(ios, 'Error in nf90_get_att: MAP_PROJECTION')

if(map_proj.eq."EQUIDISTANT CYLINDRICAL" .or. map_proj.eq."EASE V2") then

ios = nf90_inq_dimid(ftn,"east_west",ncId)
call LVT_verify(ios,&
'Error in nf90_inq_dimid in readObsDomainInput:east_west')
ios = nf90_inq_dimid(ftn,"east_west",ncId)
call LVT_verify(ios,&
'Error in nf90_inq_dimid in readObsDomainInput:east_west')

ios = nf90_inq_dimid(ftn,"north_south",nrId)
call LVT_verify(ios,&
'Error in nf90_inq_dimid in readObsDomainInput:north_south')
ios = nf90_inq_dimid(ftn,"north_south",nrId)
call LVT_verify(ios,&
'Error in nf90_inq_dimid in readObsDomainInput:north_south')

ios = nf90_inquire_dimension(ftn,ncId, len=lisdaobs(i)%nc)
call LVT_verify(ios,&
'Error in nf90_inquire_dimension in readObsDomainInput:ncId')
ios = nf90_inquire_dimension(ftn,ncId, len=lisdaobs(i)%nc)
call LVT_verify(ios,&
'Error in nf90_inquire_dimension in readObsDomainInput:ncId')

ios = nf90_inquire_dimension(ftn,nrId, len=lisdaobs(i)%nr)
call LVT_verify(ios,&
'Error in nf90_inquire_dimension in readObsDomainInput:nrId')
ios = nf90_inquire_dimension(ftn,nrId, len=lisdaobs(i)%nr)
call LVT_verify(ios,&
'Error in nf90_inquire_dimension in readObsDomainInput:nrId')

ios = nf90_get_att(ftn, NF90_GLOBAL, 'MAP_PROJECTION',map_proj)
call LVT_verify(ios, 'Error in nf90_get_att: MAP_PROJECTION')
!ios = nf90_get_att(ftn, NF90_GLOBAL, 'MAP_PROJECTION',map_proj)
!call LVT_verify(ios, 'Error in nf90_get_att: MAP_PROJECTION')

ios = nf90_get_att(ftn, NF90_GLOBAL, 'SOUTH_WEST_CORNER_LAT',&
stlat)
call LVT_verify(ios, &
'Error in nf90_get_att: SOUTH_WEST_CORNER_LAT')
ios = nf90_get_att(ftn, NF90_GLOBAL, 'SOUTH_WEST_CORNER_LAT',&
stlat)
call LVT_verify(ios, &
'Error in nf90_get_att: SOUTH_WEST_CORNER_LAT')

ios = nf90_get_att(ftn, NF90_GLOBAL, 'SOUTH_WEST_CORNER_LON',&
stlon)
call LVT_verify(ios, &
'Error in nf90_get_att: SOUTH_WEST_CORNER_LON')
ios = nf90_get_att(ftn, NF90_GLOBAL, 'SOUTH_WEST_CORNER_LON',&
stlon)
call LVT_verify(ios, &
'Error in nf90_get_att: SOUTH_WEST_CORNER_LON')

ios = nf90_get_att(ftn, NF90_GLOBAL, 'DX',dx)
call LVT_warning(ios, 'Error in nf90_get_att: DX')
if(ios.ne.0) then
dx = 0
endif
ios = nf90_get_att(ftn, NF90_GLOBAL, 'DX',dx)
call LVT_warning(ios, 'Error in nf90_get_att: DX')
if(ios.ne.0) then
dx = 0
endif

ios = nf90_get_att(ftn, NF90_GLOBAL, 'DY',dy)
call LVT_warning(ios, 'Error in nf90_get_att: DY')
if(ios.ne.0) then
dy = 0
endif
ios = nf90_get_att(ftn, NF90_GLOBAL, 'DY',dy)
call LVT_warning(ios, 'Error in nf90_get_att: DY')
if(ios.ne.0) then
dy = 0
endif

ios = nf90_get_att(ftn,NF90_GLOBAL,'GRIDTYPE',gridtype)
call LVT_warning(ios,'Error in nf90_get_att: GRIDTYPE')
ios = nf90_get_att(ftn,NF90_GLOBAL,'GRIDTYPE',gridtype)
call LVT_warning(ios,'Error in nf90_get_att: GRIDTYPE')

if(map_proj.eq."EQUIDISTANT CYLINDRICAL") then

gridDesci(1) = 0
gridDesci(2) = lisdaobs(i)%nc
gridDesci(3) = lisdaobs(i)%nr
gridDesci(4) = stlat
gridDesci(5) = stlon
gridDesci(7) = stlat + (lisdaobs(i)%nr-1)*dy
gridDesci(8) = stlon + (lisdaobs(i)%nc-1)*dx
gridDesci(9) = dx
if(map_proj.eq."EQUIDISTANT CYLINDRICAL") then

gridDesci(1) = 0
gridDesci(2) = lisdaobs(i)%nc
gridDesci(3) = lisdaobs(i)%nr
gridDesci(4) = stlat
gridDesci(5) = stlon
gridDesci(7) = stlat + (lisdaobs(i)%nr-1)*dy
gridDesci(8) = stlon + (lisdaobs(i)%nc-1)*dx
gridDesci(9) = dx

if(gridDesci(1).eq.0) then
gridDesci(10) = dy
if(gridDesci(1).eq.0) then
gridDesci(10) = dy
gridDesci(6) = 128
gridDesci(11) = 64
gridDesci(20) = 64
endif
if(gridDesci(7).lt.gridDesci(4)) then
write(LVT_logunit,*) '[ERR] lat2 must be greater than lat1'
write(LVT_logunit,*) '[ERR] ',gridDesci(7),&
gridDesci(4)
write(LVT_logunit,*) '[ERR] Stopping run...'
call LVT_endrun
endif
if(gridDesci(8).lt.gridDesci(5)) then
write(LVT_logunit,*) '[ERR] lon2 must be greater than lon1'
write(LVT_logunit,*) '[ERR] ',gridDesci(8),&
gridDesci(5)
write(LVT_logunit,*) '[ERR] Stopping run...'
call LVT_endrun
endif
elseif(map_proj.eq."EASE V2") then

gridDesci(1) = 9
gridDesci(2) = lisdaobs(i)%nc
gridDesci(3) = lisdaobs(i)%nr
gridDesci(4) = stlat
gridDesci(5) = stlon
gridDesci(6) = 128
gridDesci(11) = 64

gridDesci(20) = 64

if(gridtype.eq."M36") then
gridDesci(9) = 4
gridDesci(10) = 0.36
dx = 0.36
dy = 0.36
elseif(gridtype.eq."M09") then
gridDesci(9) = 5
gridDesci(10) = 0.09
dx = 0.09
dy = 0.09
endif
endif
if(gridDesci(7).lt.gridDesci(4)) then
write(LVT_logunit,*) '[ERR] lat2 must be greater than lat1'
write(LVT_logunit,*) '[ERR] ',gridDesci(7),&
gridDesci(4)
write(LVT_logunit,*) '[ERR] Stopping run...'
call LVT_endrun
endif
if(gridDesci(8).lt.gridDesci(5)) then
write(LVT_logunit,*) '[ERR] lon2 must be greater than lon1'
write(LVT_logunit,*) '[ERR] ',gridDesci(8),&
gridDesci(5)
write(LVT_logunit,*) '[ERR] Stopping run...'
call LVT_endrun
endif
elseif(map_proj.eq."EASE V2") then

gridDesci(1) = 9
gridDesci(2) = lisdaobs(i)%nc
gridDesci(3) = lisdaobs(i)%nr
gridDesci(4) = stlat
gridDesci(5) = stlon
gridDesci(6) = 128

gridDesci(20) = 64

if(gridtype.eq."M36") then
gridDesci(9) = 4
gridDesci(10) = 0.36
dx = 0.36
dy = 0.36
elseif(gridtype.eq."M09") then
gridDesci(9) = 5
gridDesci(10) = 0.09
dx = 0.09
dy = 0.09
endif


elseif(map_proj.eq."LAMBERT CONFORMAL") then
ios = nf90_inq_dimid(ftn,"east_west",ncId)
call LVT_verify(ios,&
'Error in nf90_inq_dimid in readObsDomainInput:east_west')

ios = nf90_inq_dimid(ftn,"north_south",nrId)
call LVT_verify(ios,&
'Error in nf90_inq_dimid in readObsDomainInput:north_south')

ios = nf90_inquire_dimension(ftn,ncId, len=lisdaobs(i)%nc)
call LVT_verify(ios,&
'Error in nf90_inquire_dimension in readObsDomainInput:ncId')

ios = nf90_inquire_dimension(ftn,nrId, len=lisdaobs(i)%nr)
call LVT_verify(ios,&
'Error in nf90_inquire_dimension in readObsDomainInput:nrId')

ios = nf90_get_att(ftn, NF90_GLOBAL, 'SOUTH_WEST_CORNER_LAT',&
stlat)
call LVT_verify(ios, &
'Error in nf90_get_att: SOUTH_WEST_CORNER_LAT')

ios = nf90_get_att(ftn, NF90_GLOBAL, 'SOUTH_WEST_CORNER_LON',&
stlon)
call LVT_verify(ios, &
'Error in nf90_get_att: SOUTH_WEST_CORNER_LON')

ios = nf90_get_att(ftn, NF90_GLOBAL, 'TRUELAT1',&
truelat1)
call LVT_verify(ios, &
'Error in nf90_get_att: TRUELAT1')

ios = nf90_get_att(ftn, NF90_GLOBAL, 'TRUELAT2',&
truelat2)
call LVT_verify(ios, &
'Error in nf90_get_att: TRUELAT1')

ios = nf90_get_att(ftn, NF90_GLOBAL, 'STANDARD_LON',&
standard_lon)
call LVT_verify(ios, &
'Error in nf90_get_att: STANDARD_LON')

ios = nf90_get_att(ftn, NF90_GLOBAL, 'DX',&
domain_resolution)
call LVT_verify(ios, &
'Error in nf90_get_att: DX')

ios = nf90_get_att(ftn,NF90_GLOBAL,'GRIDTYPE',gridtype)
call LVT_warning(ios,'Error in nf90_get_att: GRIDTYPE')

gridDesci = 0
gridDesci(1) = 3 ! Lambert conic conformal grid
gridDesci(2) = lisdaobs(i)%nc
gridDesci(3) = lisdaobs(i)%nr
gridDesci(4) = stlat ! latitude of origin -- LL Lat
gridDesci(5) = stlon ! longitude of origin -- LL Lon
gridDesci(6) = 8 ! Set for Lambert in core/LDT_domainMod.F90 (line 2658)
gridDesci(7) = truelat2 ! true lat2
gridDesci(8) = domain_resolution ! grid spacing in km
gridDesci(9) = domain_resolution ! grid spacing in km
gridDesci(10) = truelat1 ! true lat1
gridDesci(11) = standard_lon ! standard long
gridDesci(20) = 0.0


else
write(LVT_logunit,*) '[ERR] Map projection ',trim(map_proj)
write(LVT_logunit,*) '[ERR] not currently supported for LIS DAOBS plugin'
Expand Down

0 comments on commit 9d9c5fe

Please sign in to comment.