subroutine COF2GRD(LUN1,IMAX,JMAX,KMAX,TGRID,ZGRID, &
UGRID,VGRID,QGRID,CWMGRID,PRESGRID,PINTGRID)
!
!$$$ SUBPROGRAM DOCUMENTATION BLOCK
!C . . . .
! SUBPROGRAM: COF2GRD CONVERT ONE RECORD OF SIGMA COEFF FILE
!C TO LAT/LON GRID
! PRGMMR: ROGERS ORG: W/NP22 DATE: 99-01-28
!
! ABSTRACT: CONVERT SIGMA COEFFICIENT RECORD TO GRID SPACE USING
! SPLIB ROUTINES. THESE ROUTINES WILL RETURN A GLOBAL
! LAT/LON GRID WHOSE RESOLUTION IS DETERMINED BY THE
! NUMBER OF GRID POINTS. THEN, THE RELEVENT SUBSET FOR
! WHICH WE HAVE HIGH-RES OROGRAPHY IS EXTRACTED (DIMENSION
! OF BOTH THE EXTRACTED GRID AND GLOBAL GRID SET IN
! parmanl FILE)
!
! PROGRAM HISTORY LOG:
! 99-01-28 ROGERS
!
! USAGE: CALL COF2GRD(LUN1,JROMB,JCAP,KMAX)
! JCAP)
!
! INPUT ARGUMENT LIST:
! LUN1 - FORTRAN UNIT FOR SIGMA FILE
! JROMB - SPECTRAL DOMAIN SHAPE (0 FOR TRIANGULAR,
! 1 FOR RHOMBOIDAL)
! JCAP - SPECTRAL TRUNCATION
!
! OUTPUT FILES:
! KMAX - NUMBER OF SIGMA LEVELS IN GLOBAL MODEL
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN-90
! MACHINE: CRAY C-90
!
!$$$
USE SIGIO_MODULE
IMPLICIT NONE
integer,parameter::real_32=selected_real_kind(6,30)
INTEGER, INTENT(IN) :: LUN1
INTEGER, INTENT(IN) :: IMAX,JMAX, KMAX
REAL(REAL_32), DIMENSION(IMAX,JMAX,KMAX), INTENT(OUT) :: TGRID, UGRID, VGRID
REAL(REAL_32), DIMENSION(IMAX,JMAX,KMAX), INTENT(OUT) :: QGRID, CWMGRID, PRESGRID
REAL(REAL_32), DIMENSION(IMAX,JMAX,KMAX+1), INTENT(OUT) :: PINTGRID, ZGRID
TYPE(sigio_head), SAVE:: HEAD
TYPE(sigio_data), SAVE:: DATA
INTEGER :: I,J,IRET,L,IER,K
INTEGER :: II,JJ,LL
real, parameter:: G=9.81
real, parameter:: R=287.04
real, parameter:: P00=1.0E5
real, parameter:: CP=1004.6
real, parameter:: CPOG=CP/G
real, parameter:: CAPA=R/CP
CHARACTER FILENAME*80
REAL(REAL_32) :: EXNT
REAL(REAL_32),DIMENSION(IMAX,JMAX) :: EXNL
REAL(REAL_32),DIMENSION(IMAX,JMAX,KMAX) :: PD
REAL(REAL_32),DIMENSION(IMAX,JMAX,3) :: PGRID
write(filename,633) LUN1
633 format('fort.',I2.2)
CLOSE(LUN1)
CALL sigio_srohdc(LUN1,trim(filename),head,data,iret)
if (iret /= 0) then
write(0,*) 'iret from sigio_srohdc: ', iret, LUN1, trim(filename)
call mpi_abort(iret)
end if
if (KMAX /= head%levs) then
write(0,*) ' KMAX /= head%levs : KMAX=',KMAX,' head%levs-',head%levs
call mpi_abort(iret)
end if
! hs (surface topo)
CALL SPTRAN(0,head%JCAP,0,IMAX,JMAX,1,0,0,-IMAX,IMAX, &
0,0,0,0,1,data%hs,PGRID(1,JMAX,2),PGRID(1,1,2),1)
! ps and midlayer P values
CALL SPTRAN(0,head%JCAP,0,IMAX,JMAX,1,0,0,-IMAX,IMAX, &
0,0,0,0,1,data%ps,PGRID(1,JMAX,1),PGRID(1,1,1),1)
DO J=1,JMAX
DO I=1,IMAX
PGRID(I,J,1)=1000.*EXP(PGRID(I,J,1))
ENDDO
ENDDO
CALL sigio_modpr(IMAX*JMAX,IMAX*JMAX,head%levs, &
head%nvcoord,head%idvc,head%idsl, &
head%vcoord,iret,ps=PGRID(1,1,1),pm=PRESGRID, &
pd=PD)
! temperatures
DO J=1,JMAX
DO I=1,IMAX
PINTGRID(I,J,1)=PGRID(I,J,1)
ENDDO
ENDDO
DO L=1,head%levs
DO JJ=1,JMAX
DO II=1,IMAX
PINTGRID(II,JJ,L+1)=PINTGRID(II,JJ,L)-PD(II,JJ,L)
IF (L .eq. head%levs .and. PINTGRID(II,JJ,L+1) .lt. 20.) then
PINTGRID(II,JJ,L+1)=20.
ENDIF
ENDDO
ENDDO
ENDDO
DO LL=1,KMAX
CALL SPTRAN(0,head%JCAP,0,IMAX,JMAX,1,0,0, -IMAX,IMAX, &
0,0,0,0,1,data%t(:,LL),TGRID(1,JMAX,LL),TGRID(1,1,LL),1)
END DO
!
! READ DIVERGENCE AND VORTICITY COEFFICIENTS
!
DO L=1,KMAX
CALL SPTRANV(0,head%jcap,0,IMAX,JMAX,1,0,0,-IMAX,IMAX, &
0,0,0,0,1,data%d(:,L),data%z(:,L),UGRID(1,JMAX,L),UGRID(1,1,L), &
VGRID(1,JMAX,L),VGRID(1,1,L),1)
ENDDO
!
! READ SPECIFIC HUMIDITY COEFFICIENTS
!
DO L=1,KMAX
CALL SPTRAN(0,head%JCAP,0,IMAX,JMAX,1,0,0,-IMAX,IMAX, &
0,0,0,0,1,data%q(:,L,1),QGRID(1,JMAX,L),QGRID(1,1,L),1)
ENDDO
DO K = 1, KMAX
DO J = 1, JMAX
DO I = 1, IMAX
QGRID(I,J,K) = AMAX1(QGRID(I,J,K),1.0E-8)
ENDDO
ENDDO
ENDDO
CWMGRID=-9999.
DO L=1,KMAX
CALL SPTRAN(0,head%JCAP,0,IMAX,JMAX,1,0,0,-IMAX,IMAX, &
0,0,0,0,1,data%q(:,L,3),CWMGRID(1,JMAX,L),CWMGRID(1,1,L),1)
ENDDO
!!! compute ZGRID following the approach in SIG2HYB.f of the mkbnd code
CALL FLIP(PRESGRID,IMAX,JMAX,KMAX)
CALL FLIP(TGRID,IMAX,JMAX,KMAX)
CALL FLIP(UGRID,IMAX,JMAX,KMAX)
CALL FLIP(VGRID,IMAX,JMAX,KMAX)
CALL FLIP(QGRID,IMAX,JMAX,KMAX)
CALL FLIP(PINTGRID,IMAX,JMAX,KMAX+1)
!! has zsfc
DO 250 LL=1,KMAX
L=KMAX+1-LL
DO 250 J=1,JMAX
DO 250 I=1,IMAX
IF(L.EQ.KMAX) then
EXNL(I,J)=(PGRID(I,J,1)/P00)**CAPA
ZGRID(I,J,L+1)=PGRID(I,J,2)
ENDIF
EXNT=(PINTGRID(I,J,L)/P00)**CAPA
ZGRID(I,J,L)=ZGRID(I,J,L+1)-CPOG*TGRID(I,J,L)*(EXNT-EXNL(I,J)) &
*(P00/PRESGRID(I,J,L))**CAPA
EXNL(I,J)=EXNT
250 CONTINUE
!!!! end ZGRID compute
! Devirtualize GFS temperature
DO K = 1, KMAX
DO J = 1, JMAX
DO I = 1, IMAX
TGRID(I,J,K)=TGRID(I,J,K)/(1.0+0.602*QGRID(I,J,K))
ENDDO
ENDDO
ENDDO
END subroutine COF2GRD
SUBROUTINE FLIP (AE,IMAX,JMAX,KMAX)
DIMENSION AE(IMAX,JMAX,KMAX),TEMP(IMAX,JMAX,KMAX)
!************************************************************************
DO 100 KN=1,KMAX
K=KMAX-KN+1
DO 50 J=1,JMAX
DO 50 I=1,IMAX
TEMP(I,J,K)=AE(I,J,KN)
50 CONTINUE
100 CONTINUE
DO KN=1,KMAX
DO J=1,JMAX
DO I=1,IMAX
AE(I,J,KN)=TEMP(I,J,KN)
ENDDO
ENDDO
ENDDO
RETURN
END SUBROUTINE FLIP