Skip to content

Commit

Permalink
Post wafs (NOAA-EMC#28)
Browse files Browse the repository at this point in the history 
* Output 3D WAFS products to standard atmospheric pressuer levels, not on
model pressure levels.
1. 3D WAFS products include U/V/T/RH/VVEL/ABSV, icing and GTG.
2. GTG processing becomes part of WAFS product processing.
3. New code sorc/ncep_post.fd/MDL2STD_P.f is in charge of this.
4. sorc/ncep_post.fd/FDLVL.f is modified to process multiple 3D fields
   on the same vertical levels.
5. Modified for control file in sorc/ncep_post.fd/SET_LVLSXML.f
6. New control files for WAFS analysis, FF 03-36 and FF >=42
7. Icing and GTG is moved from MDL2P to MDL2STD_P.
8. sorc/ncep_post.fd/CALVOR.f is modified to take of SPVAL input values.
9. parm/post_avblflds.xml moves icing and GTG from ISOBARIC_SFC to
   STD_ISOBARIC_SFC, add U/V/T/RH/VVEL/ABSV on STD_ISOBARIC_SFC
10. scripts/exgfs_nceppost.sh.ecf removed GTG and added WAFS
11. Added sorc/ncep_post_gtg.fd and modified sorc/build_ncep_post.sh
    to keep GTG code private

* Make a clean copy of UPP WAFS branch for general UPP

* Add underground interpolation to WAFS related fields and
FDLVL_MASS subroutine is modified
  • Loading branch information
@YaliMao-NOAA @WenMeng-NOAA
YaliMao-NOAA authored and WenMeng-NOAA committed Dec 10, 2019
1 parent dc48bdb commit 3b7f408
Show file tree
Hide file tree
Showing 2 changed files with 138 additions and 14 deletions.
116 changes: 106 additions & 10 deletions sorc/ncep_post.fd/FDLVL.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -690,6 +690,9 @@ SUBROUTINE FDLVL_UV(ITYPE,NFD,HTFD,UFD,VFD)
ELSEIF (L == LM) THEN
UFD(I,J,IFD)=UH(I,J,L)
VFD(I,J,IFD)=VH(I,J,L)
ELSE ! Underground
UFD(I,J,IFD)=UH(I,J,LM)
VFD(I,J,IFD)=VH(I,J,LM)
ENDIF
!
enddo ! end of i loop
Expand Down Expand Up @@ -792,7 +795,7 @@ SUBROUTINE FDLVL_UV(ITYPE,NFD,HTFD,UFD,VFD)
RETURN
END

SUBROUTINE FDLVL_MASS(ITYPE,NFD,HTFD,NIN,QIN,QFD)
SUBROUTINE FDLVL_MASS(ITYPE,NFD,PTFD,HTFD,NIN,QIN,QTYPE,QFD)
!$$$ SUBPROGRAM DOCUMENTATION BLOCK
! . . .
! SUBPROGRAM: FDLVL_MASS COMPUTES FD LEVEL FOR MASS VARIABLES
Expand Down Expand Up @@ -837,15 +840,22 @@ SUBROUTINE FDLVL_MASS(ITYPE,NFD,HTFD,NIN,QIN,QFD)
! WITH THE SAME LEVELS AT ONE TIME
! DUST=>AERFD CAN BE PROCESSED WHEN NIN=NBIN_DU
!
! USAGE: CALL FDLVL_MASS(ITYPE,NFD,HTFD,NIN,QIN,QFD)
! USAGE: CALL FDLVL_MASS(ITYPE,NFD,PTFD,HTFD,NIN,QIN,QTYPE,QFD)
! INPUT ARGUMENT LIST:
! ITYPE - FLAG THAT DETERMINES WHETHER MSL (1) OR AGL (2)
! LEVELS ARE USED.
! NFD - NUMBER OF FD LEVELS
! HTFD - FD LEVELS
! PTFD - FD PRESSURE LEVELS
! HTFD - FD HEIGHT LEVELS
! NIN - NUMBER OF INPUT FIELDS
! QIN - ARRAY OF MASS POINT VALUE ON MODEL LEVELS
!
! QIN - ARRAY OF MASS POINT VALUE ON MODEL LEVELS
! QTYPE - CHARACTER ARRAY OF VARIABLE TYPE TO DIFFERENTIATE UNDERGROUND INTERPOLATION
! C-5 Cloud Species
! K-TURBULENT KINETIC ENERGY
! Q-Specific Humidity
! T-Temperature,
! W-Vertical Velocity or Omega
!
! OUTPUT ARGUMENT LIST:
! QFD - ARRAY OF MASS POINT VALUE ON FD LEVELS.
!
Expand Down Expand Up @@ -884,25 +894,31 @@ SUBROUTINE FDLVL_MASS(ITYPE,NFD,HTFD,NIN,QIN,QFD)
! endif
!

use vrbls3d, only: ZMID
use vrbls3d, only: T,Q,ZMID,PMID,PINT,ZINT
use vrbls2d, only: FIS
use masks, only: LMH
use params_mod, only: GI, G
use params_mod, only: GI, G, GAMMA,PQ0, A2, A3, A4, RHMIN,RGAMOG
use ctlblk_mod, only: JSTA, JEND, SPVAL, JSTA_2L, JEND_2U, LM, JSTA_M, &
JEND_M, IM, JM,global
JEND_M, IM, JM,global,MODELNAME
use gridspec_mod, only: GRIDTYPE
use physcons_post,only: CON_FVIRT, CON_ROG, CON_EPS, CON_EPSM1
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
implicit none
!
! SET NUMBER OF FD LEVELS.
!
! DECLARE VARIABLES
!
real,parameter:: zshul=75.,tvshul=290.66
real,external :: fpvsnew

integer,intent(in) :: ITYPE(NFD)
integer,intent(in) :: NFD ! coming from calling subroutine
real, intent(in) :: PTFD(NFD)
real,intent(in) :: HTFD(NFD)
integer,intent(in) :: NIN
real,intent(in) :: QIN(IM,JSTA:JEND,LM,NIN)
character, intent(in) :: QTYPE(NIN)
real,intent(out) :: QFD(IM,JSTA:JEND,NFD,NIN)

!
Expand All @@ -912,6 +928,9 @@ SUBROUTINE FDLVL_MASS(ITYPE,NFD,HTFD,NIN,QIN,QFD)
integer I,J,L,LLMH,IFD,N
integer ISTART,ISTOP,JSTART,JSTOP
real htt,htsfc,dz,rdz,delq,htabh

real :: tvu,tvd,gammas,part,ES,QSAT,RHL,PL,ZL,TL,QL
real :: TVRL,TVRBLO,TBLO,QBLO
!
!****************************************************************
! START FDLVL_MASS HERE
Expand Down Expand Up @@ -942,6 +961,7 @@ SUBROUTINE FDLVL_MASS(ITYPE,NFD,HTFD,NIN,QIN,QFD)
END IF

DO IFD = 1, NFD

!
! MSL FD LEVELS
!
Expand Down Expand Up @@ -975,7 +995,7 @@ SUBROUTINE FDLVL_MASS(ITYPE,NFD,HTFD,NIN,QIN,QFD)
exit
END IF

enddo ! end of l loop
ENDDO ! end of L loop
!
! COMPUTE Q AT FD LEVELS.
!
Expand All @@ -997,7 +1017,83 @@ SUBROUTINE FDLVL_MASS(ITYPE,NFD,HTFD,NIN,QIN,QFD)
DO N = 1, NIN
QFD(I,J,IFD,N) = QIN(I,J,L,N)
ENDDO
ENDIF
ELSE ! Underground
DO N = 1, NIN
! Deduce T and Q differently by different models
IF(MODELNAME == 'GFS')THEN ! GFS deduce T using Shuell
if(QTYPE(N) == "T" .or. QTYPE(N) == "Q") then
tvu = T(I,J,LM) * (1.+con_fvirt*Q(I,J,LM))
if(ZMID(I,J,LM) > zshul) then
tvd = tvu + gamma*ZMID(I,J,LM)
if(tvd > tvshul) then
if(tvu > tvshul) then
tvd = tvshul - 5.e-3*(tvu-tvshul)*(tvu-tvshul)
else
tvd = tvshul
endif
endif
gammas = (tvu-tvd)/ZMID(I,J,LM)
else
gammas = 0.
endif
part = con_rog*(LOG(PTFD(IFD))-LOG(PMID(I,J,LM)))
part = ZMID(I,J,LM) - tvu*part/(1.+0.5*gammas*part)
part = T(I,J,LM) - gamma*(part-ZMID(I,J,LM))

if(QTYPE(N) == "T") QFD(I,J,IFD,N) = part

if(QTYPE(N) == "Q") then

! Compute RH at lowest model layer because Iredell and Chuang decided to compute
! underground GFS Q to maintain RH
ES = min(FPVSNEW(T(I,J,LM)), PMID(I,J,LM))
QSAT = CON_EPS*ES/(PMID(I,J,LM)+CON_EPSM1*ES)
RHL = Q(I,J,LM)/QSAT
! compute saturation water vapor at isobaric level
ES = min(FPVSNEW(part), PTFD(IFD))
QSAT = CON_EPS*ES/(PTFD(IFD)+CON_EPSM1*ES)
! Q at isobaric level is computed by maintaining constant RH
QFD(I,J,IFD,N) = RHL*QSAT
endif
endif

ELSE
if(QTYPE(N) == "T" .or. QTYPE(N) == "Q") then
PL = PINT(I,J,LM-1)
ZL = ZINT(I,J,LM-1)
TL = 0.5*(T(I,J,LM-2)+T(I,J,LM-1))
QL = 0.5*(Q(I,J,LM-2)+Q(I,J,LM-1))

QSAT = PQ0/PL*EXP(A2*(TL-A3)/(TL-A4))
RHL = QL/QSAT
!
IF(RHL > 1.)THEN
RHL = 1.
QL = RHL*QSAT
ENDIF
!
IF(RHL < RHmin)THEN
RHL = RHmin
QL = RHL*QSAT
ENDIF
!
TVRL = TL*(1.+0.608*QL)
TVRBLO = TVRL*(PTFD(IFD)/PL)**RGAMOG
TBLO = TVRBLO/(1.+0.608*QL)

QSAT = PQ0/PTFD(IFD)*EXP(A2*(TBLO-A3)/(TBLO-A4))
if(QTYPE(N) == "T") QFD(I,J,IFD,N) = TBLO
QBLO = RHL*QSAT
if(QTYPE(N) == "Q") QFD(I,J,IFD,N) = MAX(1.E-12,QBLO)
endif
END IF ! endif loop for deducing T and Q differently for GFS

if(QTYPE(N) == "W") QFD(I,J,IFD,N)=QIN(I,J,LM,N) ! W OMGA
if(QTYPE(N) == "K") QFD(I,J,IFD,N)= max(0.0,0.5*(QIN(I,J,LM,N)+QIN(I,J,LM-1,N))) ! TKE
if(QTYPE(N) == "C") QFD(I,J,IFD,N)=0.0 ! Hydrometeor fields
END DO

ENDIF ! Underground

!
! COMPUTE FD LEVEL Q AT NEXT K.
Expand Down
36 changes: 32 additions & 4 deletions sorc/ncep_post.fd/MDL2STD_P.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,7 @@ SUBROUTINE MDL2STD_P()
REAL, allocatable :: HTFDCTL(:)
integer, allocatable :: ITYPEFDLVLCTL(:)
real, allocatable :: QIN(:,:,:,:), QFD(:,:,:,:)
character, allocatable :: QTYPE(:)
real, allocatable :: VAR3D1(:,:,:), VAR3D2(:,:,:)

integer, parameter :: NFDMAX=50 ! Max number of fields with the same HTFDCTL
Expand Down Expand Up @@ -118,6 +119,8 @@ SUBROUTINE MDL2STD_P()
DO i = 1, NFDCTL
HTFDCTL(i)=P2H(HTFDCTL(i)/100.)
ENDDO
if(allocated(VAR3D1)) deallocate(VAR3D1)
if(allocated(VAR3D2)) deallocate(VAR3D2)
allocate(VAR3D1(IM,JSTA_2L:JEND_2U,NFDCTL))
allocate(VAR3D2(IM,JSTA_2L:JEND_2U,NFDCTL))
VAR3D1=SPVAL
Expand Down Expand Up @@ -203,74 +206,90 @@ SUBROUTINE MDL2STD_P()
! HGT(TO BE FIXED VALUES)
! RH ABSV (TO BE CACULATED)

if(allocated(QIN)) deallocate(QIN)
if(allocated(QTYPE)) deallocate(QTYPE)
ALLOCATE(QIN(IM,JSTA:JEND,LM,NFDMAX))
ALLOCATE(QTYPE(NFDMAX))

! INITIALIZE INPUTS
nFDS = 0
IF(IGET(450) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 450
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=icing_gfip(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="O"
end if
IF(IGET(480) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 480
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=icing_gfis(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="O"
end if
IF(IGET(464) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 464
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=gtg(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="O"
end if
IF(IGET(465) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 465
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=catedr(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="O"
end if
IF(IGET(466) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 466
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=mwt(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="O"
end if
IF(IGET(519) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 519
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=T(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="T"
end if
IF(IGET(522) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 522
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=Q(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="Q"
end if
IF(IGET(524) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 524
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=OMGA(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="W"
end if
IF(IGET(526) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 526
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=QQW(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="C"
end if
IF(IGET(527) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 527
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=QQR(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="C"
end if
IF(IGET(528) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 528
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=QQS(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="C"
end if
IF(IGET(529) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 529
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=QQG(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="C"
end if
IF(IGET(530) > 0) THEN
nFDS = nFDS + 1
IDS(nFDS) = 530
QIN(1:IM,JSTA:JEND,1:LM,nFDS)=QQI(1:IM,JSTA:JEND,1:LM)
QTYPE(nFDS)="C"
end if

! FOR WAFS, ALL LEVLES OF DIFFERENT VARIABLES ARE THE SAME, USE ANY
Expand All @@ -289,10 +308,11 @@ SUBROUTINE MDL2STD_P()
HTFDCTL(i)=P2H(HTFDCTL(i)/100.)
ENDDO

if(allocated(QFD)) deallocate(QFD)
ALLOCATE(QFD(IM,JSTA:JEND,NFDCTL,nFDS))
QFD=SPVAL

call FDLVL_MASS(ITYPEFDLVLCTL,NFDCTL,HTFDCTL,nFDS,QIN,QFD)
call FDLVL_MASS(ITYPEFDLVLCTL,NFDCTL,pset%param(N)%level,HTFDCTL,nFDS,QIN,QTYPE,QFD)

! Adjust values before output
N1 = -1
Expand Down Expand Up @@ -415,6 +435,7 @@ SUBROUTINE MDL2STD_P()
ENDDO

DEALLOCATE(QIN,QFD)
DEALLOCATE(QTYPE)

! STEP 3 - MASS FIELDS CALCULATION
! HGT(TO BE FIXED VALUES)
Expand Down Expand Up @@ -472,16 +493,23 @@ SUBROUTINE MDL2STD_P()
HTFDCTL(i)=P2H(HTFDCTL(i)/100.)
ENDDO

if(allocated(QIN)) deallocate(QIN)
if(allocated(QTYPE)) deallocate(QTYPE)
ALLOCATE(QIN(IM,JSTA:JEND,LM,2))
ALLOCATE(QTYPE(2))
QIN(1:IM,JSTA:JEND,1:LM,1)=T(1:IM,JSTA:JEND,1:LM)
QIN(1:IM,JSTA:JEND,1:LM,2)=Q(1:IM,JSTA:JEND,1:LM)
QTYPE(1)="T"
QTYPE(2)="Q"

if(allocated(QFD)) deallocate(QFD)
ALLOCATE(QFD(IM,JSTA:JEND,NFDCTL,2))
QFD=SPVAL

call FDLVL_MASS(ITYPEFDLVLCTL,NFDCTL,HTFDCTL,2,QIN,QFD)
HTFDCTL=pset%param(N)%level ! Save back to pressure
print *, "wafs levels",pset%param(N)%level
call FDLVL_MASS(ITYPEFDLVLCTL,NFDCTL,pset%param(N)%level,HTFDCTL,2,QIN,QTYPE,QFD)

HTFDCTL=pset%param(N)%level ! Save back to pressure

DO IFD = 1,NFDCTL
IF (LVLS(IFD,IGET(iID)) > 0) THEN
Expand Down Expand Up @@ -533,7 +561,7 @@ SUBROUTINE MDL2STD_P()
ENDIF
ENDDO
deallocate(QIN,QFD)

deallocate(QTYPE)
ENDIF

ENDIF
Expand Down

0 comments on commit 3b7f408

Please sign in to comment.