!-----------------------------------------------------------------------
!
PROGRAM PREPROC_LATLON_SM
!
!-----------------------------------------------------------------------
!*** Dmitry's coupler needs to read in the geographic lat/lon of
!*** the H and V points on the B-grid that spans the upper parent
!*** domain and is at the resolution of an outer nest. It needs
!*** the sea mask on that grid as well.
!-----------------------------------------------------------------------
!
USE netcdf
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!-----------------------------------------------------------------------
!
INTEGER,PARAMETER :: SINGLE=SELECTED_REAL_KIND(P=6 ,R=37)
INTEGER,PARAMETER :: DOUBLE=SELECTED_REAL_KIND(P=13,R=200)
!
!-----------------------------------------------------------------------
!
!-----------------------
!*** Set the following
!-----------------------
!
INTEGER,PARAMETER :: IM_PARENT=451 & !<-- I extent of the upper parent domain's grid
,JM_PARENT=451 & !<-- J extent of the upper parent domain's grid
,PARENT_CHILD_SPACE_RATIO=3 !<-- Ratio of parent grid increment to its child's
!
!! REAL(kind=DOUBLE),PARAMETER :: TPH0D= 18.4_double, TLM0D=-94.3_double & !<-- Geographic lat/lon (deg) of domain's center
REAL(kind=DOUBLE):: TPH0D, TLM0D !<-- Geographic lat/lon (deg) of domain's center
REAL(kind=DOUBLE), PARAMETER :: SBD =-40.5_double ,WBD =-40.5_double !<-- Rotate lat/lon (deg) of domain's south/west boundaries
!
LOGICAL :: WRITE_SINGLE=.FALSE. !<-- Write lat/lon in single (true) / double (false) precision
!
!-----------------------
!
INTEGER,PARAMETER :: IM=(IM_PARENT-1)*PARENT_CHILD_SPACE_RATIO+1 & !<-- I extent of parent domain at child resolution
,JM=(JM_PARENT-1)*PARENT_CHILD_SPACE_RATIO+1 !<-- J extent of parent domain at child resolution
!
CHARACTER(len=17) :: OUTPUT_FILENAME='coupler_grid_data' !<-- The file that will hold the final data.
!
INTEGER :: I,IRTN,J,N,NUNIT
!
REAL(kind=DOUBLE) :: DEG2RAD,DLM,DLM_HALF,DPH,DPH_HALF &
,I_MID,J_MID,PI,RAD2DEG,TLM0,TPH0
!
REAL(kind=DOUBLE),DIMENSION(IM,JM) :: GLAT_H,GLON_H & !<-- Geographic lat/lon (deg) of H points on child res grid
,GLAT_V,GLON_V & !<-- Geographic lat/lon (deg) of V points on child res grid
,TLAT_H,TLON_H & !<-- Rotated lat/lon (deg) of H points on child res grid
,TLAT_V,TLON_V !<-- Rotated lat/lon (deg) of V points on child res grid
!
REAL(kind=SINGLE),DIMENSION(:,:),ALLOCATABLE :: GLAT_H_S,GLON_H_S &
,GLAT_V_S,GLON_V_S
!
REAL,DIMENSION(IM,JM) :: SEA_MASK !<-- Sea mask (sea=1.) on child res grid
!
LOGICAL :: OPENED
!
!------------------------------------------------------------------------
!read center latitude and longitude
!-----------------------------------------------------------------------
OPEN(100,FILE='storm_center',form='FORMATTED')
READ(100,*)TPH0D
READ(100,*)TLM0D
CLOSE(100)
!-----------------------------------------------------------------------
!***********************************************************************
!-----------------------------------------------------------------------
!
PI=ACOS(-1._double)
DEG2RAD=PI/180._double
RAD2DEG=180._double/PI
!
TPH0=TPH0D*DEG2RAD
TLM0=TLM0D*DEG2RAD
!
!-----------------------------------------------------------------------
!*** Compute the constant grid increments (deg) in rotated lat/lon.
!-----------------------------------------------------------------------
!
DPH=-2._double*SBD/REAL(JM-1)
DLM=-2._double*WBD/REAL(IM-1)
!
!-----------------------------------------------------------------------
!*** Compute the rotated lat/lon of all H and V points.
!-----------------------------------------------------------------------
!
I_MID=(IM+1)*0.5_double
J_MID=(JM+1)*0.5_double
DPH_HALF=0.5_double*DPH
DLM_HALF=0.5_double*DLM
!
DO J=1,JM
DO I=1,IM
TLAT_H(I,J)=(J-J_MID)*DPH
TLAT_V(I,J)=TLAT_H(I,J)+DPH_HALF
TLON_H(I,J)=(I-I_MID)*DLM
TLON_V(I,J)=TLON_H(I,J)+DLM_HALF
ENDDO
ENDDO
!
!-----------------------------------------------------------------------
!*** Compute the geographic lat/lon of all H and V points.
!-----------------------------------------------------------------------
!
CALL GEO(TLAT_H,TLON_H,GLAT_H,GLON_H)
CALL GEO(TLAT_V,TLON_V,GLAT_V,GLON_V)
!
!-----------------------------------------------------------------------
!*** The sea mask field that spans the upper parent domain at the
!*** outer nest's resolution is already available in the NetCDF
!*** file generated for that nest by NPS so read it in.
!-----------------------------------------------------------------------
!
CALL READ_NETCDF
!
!-----------------------------------------------------------------------
!*** Write out the data.
!-----------------------------------------------------------------------
!
select_unit: DO N=51,59
INQUIRE(N,OPENED=OPENED)
IF(.NOT.OPENED)THEN
NUNIT=N
EXIT select_unit
ENDIF
ENDDO select_unit
!
OPEN(unit=NUNIT,file=OUTPUT_FILENAME,status='REPLACE' &
,form='UNFORMATTED',iostat=IRTN)
!
IF(WRITE_SINGLE)THEN
ALLOCATE(GLAT_H_S(IM,JM))
ALLOCATE(GLON_H_S(IM,JM))
ALLOCATE(GLAT_V_S(IM,JM))
ALLOCATE(GLON_V_S(IM,JM))
!
DO J=1,JM
DO I=1,IM
GLAT_H_S(I,J)=GLAT_H(I,J)
GLON_H_S(I,J)=GLON_H(I,J)
GLAT_V_S(I,J)=GLAT_V(I,J)
GLON_V_S(I,J)=GLON_V(I,J)
ENDDO
ENDDO
!
WRITE(NUNIT)GLAT_H_S
WRITE(NUNIT)GLON_H_S
WRITE(NUNIT)GLAT_V_S
WRITE(NUNIT)GLON_V_S
WRITE(NUNIT)SEA_MASK
!
ELSE
WRITE(NUNIT)GLAT_H
WRITE(NUNIT)GLON_H
WRITE(NUNIT)GLAT_V
WRITE(NUNIT)GLON_V
WRITE(NUNIT)SEA_MASK
!
ENDIF
!
CLOSE(NUNIT)
!
!-----------------------------------------------------------------------
!
CONTAINS
!
!-----------------------------------------------------------------------
!
SUBROUTINE GEO(TLAT,TLON,GLAT,GLON)
!
!-----------------------------------------------------------------------
!*** Given the grid specifications and the rotated lat/lon of each
!*** point compute the geographic lat/lon of each point.
!-----------------------------------------------------------------------
IMPLICIT NONE
!-----------------------------------------------------------------------
!
INTEGER,PARAMETER :: DOUBLE=SELECTED_REAL_KIND(P=13,R=200)
!
REAL(kind=DOUBLE),DIMENSION(IM,JM),INTENT(IN) :: TLAT,TLON
REAL(kind=DOUBLE),DIMENSION(IM,JM),INTENT(OUT) :: GLAT,GLON
!
REAL(kind=DOUBLE) :: ANUM,ARG1,DENOM,GLATR,GLONR,RELM,TLATR,TLONR
!
!-----------------------------------------------------------------------
!***********************************************************************
!-----------------------------------------------------------------------
!
DO J=1,JM
DO I=1,IM
TLATR=TLAT(I,J)*DEG2RAD
TLONR=TLON(I,J)*DEG2RAD
ARG1=SIN(TLATR)*COS(TPH0)+COS(TLATR)*SIN(TPH0)*COS(TLONR)
IF(ABS(ARG1)>1._double)THEN
ARG1=SIGN(1._double,ARG1)
ENDIF
GLATR=ASIN(ARG1)
!
ANUM=COS(TLATR)*SIN(TLONR)
DENOM=(COS(TLONR)*COS(TLATR)-SIN(TPH0)*ARG1)/COS(TPH0)
RELM=ATAN2(ANUM,DENOM)
GLONR=RELM+TLM0
IF(GLONR<-PI)GLONR=GLONR+PI*2._double
IF(GLONR> PI)GLONR=GLONR-PI*2._double
!
GLAT(I,J)=GLATR*RAD2DEG
GLON(I,J)=GLONR*RAD2DEG
ENDDO
ENDDO
!
!-----------------------------------------------------------------------
END SUBROUTINE GEO
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
SUBROUTINE READ_NETCDF
!
!-----------------------------------------------------------------------
!*** Read the sea mask data from the NPS-generated file.
!-----------------------------------------------------------------------
IMPLICIT NONE
!-----------------------------------------------------------------------
!
INTEGER :: I_COUNT_DATA,I_END,I_START &
,J_COUNT_DATA,J_END,J_START &
,NCID,NCTYPE,NDIMS,VAR_ID
!
INTEGER,DIMENSION(1:2) :: DIM_IDS
!
CHARACTER(len=2) :: SPACE_RATIO
CHARACTER(len=8) :: FILENAME
CHARACTER(len=15) :: VNAME
!
!-----------------------------------------------------------------------
!***********************************************************************
!-----------------------------------------------------------------------
!
I_START=1
I_END=IM
J_START=1
J_END=JM
!
I_COUNT_DATA=I_END-I_START+1
J_COUNT_DATA=J_END-J_START+1
!
IF(PARENT_CHILD_SPACE_RATIO<=9)THEN
WRITE(SPACE_RATIO,'(I1.1)')PARENT_CHILD_SPACE_RATIO
ELSEIF(PARENT_CHILD_SPACE_RATIO>=10)THEN
WRITE(SPACE_RATIO,'(I2.2)')PARENT_CHILD_SPACE_RATIO
ENDIF
!
FILENAME='SM_'//TRIM(SPACE_RATIO)//'.nc' !<-- Name of the existing sea mask NetCDF file
!
CALL CHECK(NF90_OPEN(FILENAME,NF90_NOWRITE,NCID)) !<-- Open the external NetCDF file with sea mask values.
!
CALL CHECK(NF90_INQUIRE_VARIABLE(NCID,3,VNAME,NCTYPE & !<-- The sea mask is the 3rd variable in the file.
,NDIMS,DIM_IDS))
!
CALL CHECK(NF90_INQ_VARID(NCID,VNAME,VAR_ID))
!
CALL CHECK(NF90_GET_VAR(NCID,VAR_ID & !<-- Extract the
,SEA_MASK(I_START:I_END,J_START:J_END) & ! sea mask values
,start=(/I_START,J_START/) & ! from the external
,count=(/I_COUNT_DATA,J_COUNT_DATA/))) ! NetCDF file.
!
!-----------------------------------------------------------------------
END SUBROUTINE READ_NETCDF
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
SUBROUTINE CHECK(RC)
!
!-----------------------------------------------------------------------
IMPLICIT NONE
!-----------------------------------------------------------------------
!
INTEGER,INTENT(IN) :: RC
!
!-----------------------------------------------------------------------
!***********************************************************************
!-----------------------------------------------------------------------
!
IF(RC/=NF90_NOERR)THEN
WRITE(*,*)TRIM(ADJUSTL(NF90_STRERROR(RC)))
! WRITE(0,11101)RC
11101 FORMAT(' ERROR: RC=',I5)
ENDIF
!
!-----------------------------------------------------------------------
!
END SUBROUTINE CHECK
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!
END PROGRAM PREPROC_LATLON_SM
!
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------