SUBROUTINE INITPOST_GFS_SIGIO(lusig,iunit,iostatusFlux,iostatusD3D,idrt,sighead)
!$$$ SUBPROGRAM DOCUMENTATION BLOCK
! . . .
! SUBPROGRAM: INITPOST INITIALIZE POST FOR RUN
! PRGRMMR: Hui-Ya Chuang DATE: 2007-03-01
!
! ABSTRACT: THIS ROUTINE INITIALIZES CONSTANTS AND
! VARIABLES AT THE START OF AN ETA MODEL OR POST
! PROCESSOR RUN.
!
! REVISION HISTORY
! 2011-02-07 Jun Wang add grib2 option
! 2013-04-19 Jun Wang add changes to read wam tracers
! 2013-05-04 Shrinivas Moorthi: real * 8 for pm1d and pi1d and pt=100hPa and some cosmetic changes
!
! USAGE: CALL INIT
! INPUT ARGUMENT LIST:
! NONE
!
! OUTPUT ARGUMENT LIST:
! NONE
!
! OUTPUT FILES:
! NONE
!
! SUBPROGRAMS CALLED:
! UTILITIES:
! NONE
! LIBRARY:
! COMMON - CTLBLK
! LOOKUP
! SOILDEPTH
!
!
! ATTRIBUTES:
! LANGUAGE: FORTRAN
! MACHINE : CRAY C-90
!$$$
use vrbls3d, only: ZINT, PINT, T, UH, VH, Q, O3, CWM, U, V, QQW, &
OMGA, PMID, PINT, ALPINT, ZMID, QQR, QQS, QQI, Q2, &
CFR, RLWTT, RSWTT, TCUCN, TCUCNS, TRAIN, EL_PBL, &
EXCH_H, VDIFFTT, VDIFFMOIS, DCONVMOIS, SCONVMOIS, &
NRADTT, O3VDIFF, O3PROD, O3TNDY, MWPV, UNKNOWN, &
VDIFFZACCE, ZGDRAG, CNVCTUMMIXING, VDIFFMACCE, &
MGDRAG, CNVCTDETMFLX,NCNVCTCFRAC, CNVCTUMFLX, &
CNVCTVMMIXING, CNVCTDMFLX, CNVCTZGDRAG, CNVCTMGDRAG
use vrbls2d, only: F, PD, FIS, PBLH, USTAR, Z0, THS, QS, TWBS, QWBS, &
AVGCPRATE, CPRATE, AVGPREC, PREC, SR, LSPA, SNO, SI, &
CLDEFI, TH10, Q10, TSHLTR, PSHLTR, QSHLTR, ALBASE, &
AVGALBEDO,AVGTCDC, CZEN, CZMEAN, MXSNAL, RADOT, &
SIGT4, VEGFRC, CFRACL, CFRACM, AVGCFRACH, CFRACH, &
AVGCFRACL, AVGCFRACM, CNVCFR, ISLOPE, CMC, GRNFLX, &
SOILTB, TG, NCFRCV, ACFRCV, ASWINTOA, ACFRST, NCFRST,&
SSROFF, BGROFF, RLWIN, RLWTOA, ALWIN, ALWOUT, ALWTOA,&
RSWIN, RSWINC, RSWOUT, ASWIN, AUVBIN, AUVBINC, &
ASWOUT, ASWTOA, ASWINC, ASWOUTC, ASWTOAC, ASWINTOA, &
AVISBEAMSWIN, AVISDIFFSWIN, AIRBEAMSWIN, AIRDIFFSWIN,&
SFCSHX, SFCLHX, SUBSHX, SNOPCX, SFCUX, SFCVX, SFCUVX,&
SFCUGS, GTAUX, SFCVGS, GTAUY, POTEVP, U10, V10, &
SMSTAV, SMSTOT, IVGTYP, ISLTYP, SFCEVP, SFCEXC, &
ACSNOW, ACSNOM, SST, QZ0, UZ0, VZ0, PTOP, HTOP, PBOT,&
HBOT, PBOT, PTOPL, PBOTL, TTOPL, PTOPM, PBOTM, TTOPM,&
PTOPH, PBOTH, TTOPH, PBLCFR, CLDWORK, RUNOFF, &
MAXTSHLTR, MINTSHLTR, DZICE, SMCWLT, SUNTIME, &
FIELDCAPA, SNOWFALL, HTOPD, HBOTD, HTOPS, HBOTS, &
CUPPT, THZ0, MAXRHSHLTR, MINRHSHLTR, U10H, V10H
use soil, only: SLDPTH, SH2O, SMC, STC
use masks, only: LMV, LMH, HTM, VTM, GDLAT, GDLON, DX, DY, HBM2, SM, SICE
use physcons, only: CON_G, CON_FVIRT, CON_RD, CON_EPS, CON_EPSM1
use masks, only: LMV, LMH, HTM, VTM, GDLAT, GDLON, DX, DY, HBM2, &
SM, SICE
use physcons, only: CON_G, CON_FVIRT, CON_RD, CON_EPS, CON_EPSM1
use params_mod, only: RTD, ERAD, DTR, TFRZ, P1000, CAPA
use lookup_mod, only: THL, PLQ, PTBL, TTBL, RDQ, RDTH, RDP, RDTHE, PL, &
QS0, SQS, STHE, THE0, TTBLQ, RDPQ, RDTHEQ, STHEQ, &
THE0Q
use ctlblk_mod, only: ME, MPI_COMM_COMP, ICNT, IDSP,JEND_M, IHRST, IMIN,&
IDAT, SDAT, IFHR, IFMIN, FILENAME, TPREC, TCLOD, &
TRDLW, TRDSW, TSRFC, TMAXMIN, TD3D, RESTRT, &
IMP_PHYSICS, DT, NUM_PROCS, LP1, PDTOP, SPVAL, PT,&
NPHS, DTQ2, ARDLW, ARDSW, ASRFC, AVRAIN, AVCNVC, &
THEAT, GDSDEGR, SPL, LSM, ALSL, IM, JM, IM_JM, &
LM, JSTA_2L, JEND_2U, NSOIL, JSTA, JEND, ICU_PHYSICS
use gridspec_mod, only: MAPTYPE, GRIDTYPE, LATSTART, LATLAST, LONSTART, &
LONLAST, CENLON, DXVAL, DYVAL, TRUELAT2, &
TRUELAT1, CENLAT
use rqstfld_mod, only: IGDS, AVBL, IQ, IS
use sigio_module, only: SIGIO_HEAD
use sfcio_module, only: sfcio_head, sfcio_data, sfcio_srohdc
! use wrf_io_flags_mod
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
implicit none
!
type(sigio_head):: sighead
!type(sigio_data):: sigdatai
type(sfcio_head):: head
type(sfcio_data):: data
!
! INCLUDE/SET PARAMETERS.
!
INCLUDE "mpif.h"
integer,parameter:: MAXPTS=1000000 ! max im*jm points
!
! real,parameter:: con_g =9.80665e+0! gravity
! real,parameter:: con_rv =4.6150e+2 ! gas constant H2O
! real,parameter:: con_rd =2.8705e+2 ! gas constant air
! real,parameter:: con_fvirt =con_rv/con_rd-1.
! real,parameter:: con_eps =con_rd/con_rv
! real,parameter:: con_epsm1 =con_rd/con_rv-1
!
! real,external::FPVSNEW
! This version of INITPOST shows how to initialize, open, read from, and
! close a NetCDF dataset. In order to change it to read an internal (binary)
! dataset, do a global replacement of _ncd_ with _int_.
integer,intent(in) :: lusig,iostatusFlux,iostatusD3D,iunit,idrt
character(len=20) :: VarName
character(len=20) :: VcoordName
integer :: Status
character startdate*19,SysDepInfo*80,cgar*1
character startdate2(19)*4
!
! NOTE: SOME INTEGER VARIABLES ARE READ INTO DUMMY ( A REAL ). THIS IS OK
! AS LONG AS REALS AND INTEGERS ARE THE SAME SIZE.
!
! ALSO, EXTRACT IS CALLED WITH DUMMY ( A REAL ) EVEN WHEN THE NUMBERS ARE
! INTEGERS - THIS IS OK AS LONG AS INTEGERS AND REALS ARE THE SAME SIZE.
LOGICAL RUNB,SINGLRST,SUBPOST,NEST,HYDRO
LOGICAL IOOMG,IOALL
logical, parameter :: debugprint = .false.
! logical, parameter :: debugprint = .true.
CHARACTER*32 LABEL
CHARACTER*40 CONTRL,FILALL,FILMST,FILTMP,FILTKE,FILUNV &
, FILCLD,FILRAD,FILSFC
CHARACTER*4 RESTHR
CHARACTER FNAME*255,ENVAR*50,sfcfilename*256
INTEGER IDATE(8),JDATE(8)
INTEGER JPDS(200),JGDS(200),KPDS(200),KGDS(200)
LOGICAL*1 LB(IM,JM)
INTEGER IRET
! REAL BUFF(IM_JM)
!
! INCLUDE COMMON BLOCKS.
!
! DECLARE VARIABLES.
!
! REAL fhour
integer nfhour ! forecast hour from nems io file
REAL RINC(5)
REAL u1d(LM), v1d(LM),omga1d(lm)
real*8 pm1d(lm), pi1d(lm+1)
REAL DUM1D (LM+1)
! REAL DUMMY ( IM, JM )
! REAL DUMMY2 ( IM, JM )
REAL, ALLOCATABLE :: dummy_h(:,:),dummy_p(:,:),dummy_px(:,:), &
dummy_py(:,:),dummy_t(:,:,:),dummy_u(:,:,:),dummy_v(:,:,:), &
dummy_d(:,:,:),dummy_trc(:,:,:,:), dummy(:,:), dummy2(:,:)
! dummy18(:,:,:),dummy19(:,:,:)
REAL, ALLOCATABLE :: dummy15(:,:),dummy16(:,:),dummy17(:,:,:)
real, allocatable :: d2d(:,:), u2d(:,:), v2d(:,:), omga2d(:,:), &
p2d(:,:), t2d(:,:), q2d(:,:), qs2d(:,:), &
cw2d(:,:), cfr2d(:,:)
real*8, allocatable :: pm2d(:,:), pi2d(:,:)
REAL FI(IM,JM,2)
! INTEGER IDUMMY(IM,JM)
!jw
integer ii,jj,js,je,iyear,imn,iday,itmp,ioutcount,istatus, &
I,J,L,ll,k,kf,irtn,igdout,n,Index,nframe, &
impf,jmpf,nframed2,iunitd3d
real TSTART,TLMH,TSPH,ES, FACT,soilayert,soilayerb,zhour,dum
real, external :: fpvsnew
real, allocatable:: glat1d(:),glon1d(:),qstl(:)
integer ierr,idum
! integer ntrac,nci,ij,ijl,j1,j2
integer lsta,lend
integer ijmc,ijxc,kna,kxa,kma
real,allocatable :: ri(:),cpi(:)
! integer ibuf(im,jsta_2l:jend_2u)
! real buf(im,jsta_2l:jend_2u),bufsoil(im,nsoil,jsta_2l:jend_2u)
real buf(im,jsta_2l:jend_2u)
real, allocatable :: buf3d(:,:,:), ta(:,:,:,:), tb(:,:,:,:)
real, allocatable :: wrk1(:,:), wrk2(:,:)
! real buf3d(im,lm,jsta:jend)
real timef
real tem,tvll,pmll
integer levs,ntrac,ncld,idvt,jcap,lnt2,ntoz,ntcw,ltrc
!
! DATA BLANK/' '/
!
!***********************************************************************
! START INIT HERE.
!
if (me == 0) WRITE(6,*)'INITPOST: ENTER INITPOST_GFS_SIGIO'
WRITE(6,*)'me=',me,'LMV=',size(LMV,1),size(LMV,2),'LMH=', &
size(LMH,1),size(LMH,2),'jsta_2l=',jsta_2l,'jend_2u=', &
jend_2u,'im=',im
!
!
! STEP 1. READ MODEL OUTPUT FILE
!
!
!***
!
! LMH always = LM for sigma-type vert coord
! LMV always = LM for sigma-type vert coord
!$omp parallel do private(i,j)
do j = jsta_2l, jend_2u
do i = 1, im
LMV ( i, j ) = lm
LMH ( i, j ) = lm
end do
end do
! write(0,*),' LM=',LM,' LP1=',LP1
! HTM VTM all 1 for sigma-type vert coord
!$omp parallel do private(i,j,l)
do l = 1, lm
do j = jsta_2l, jend_2u
do i = 1, im
HTM ( i, j, l ) = 1.0
VTM ( i, j, l ) = 1.0
end do
end do
end do
allocate (dummy(im,jm), dummy2(im,jm))
! The end j row is going to be jend_2u for all variables except for V.
JS = JSTA_2L
JE = JEND_2U
! get start date
if (me == 0)then
idate(1) = sighead%idate(4)
idate(2) = sighead%idate(2)
idate(3) = sighead%idate(3)
idate(4) = sighead%idate(1)
idate(5) = 0
nfhour = nint(sighead%fhour)
allocate(glat1d(jm),glon1d(jm))
! call splat to compute lat for gaussian grid
call splat(idrt,jm,glat1d,glon1d)
!$omp parallel do private(i,j,tem)
do j=1,jm
tem = asin(glat1d(j))*RTD
do i=1,im
dummy(i,j) = tem
dummy2(i,j) = 360./im*(i-1)
end do
end do
deallocate(glat1d,glon1d)
print*,'idate before broadcast = ',(idate(i),i=1,7)
end if
call mpi_bcast(idate(1),7,MPI_INTEGER,0,mpi_comm_comp,iret)
call mpi_bcast(nfhour,1,MPI_INTEGER,0,mpi_comm_comp,iret)
if (me == 0) then
print*,'idate after broadcast = ',(idate(i),i=1,4)
print*,'nfhour = ',nfhour
endif
! sample print point
ii = im/2
jj = jm/2
call mpi_scatterv(dummy(1,1),icnt,idsp,mpi_real &
,gdlat(1,jsta),icnt(me),mpi_real,0,MPI_COMM_COMP,ierr)
call mpi_scatterv(dummy2(1,1),icnt,idsp,mpi_real &
,gdlon(1,jsta),icnt(me),mpi_real,0,MPI_COMM_COMP,ierr)
! write(0,*)'before call EXCH,mype=',me,'max(gdlat)=',maxval(gdlat),&
! 'max(gdlon)=', maxval(gdlon)
CALL EXCH(gdlat(1,JSTA_2L))
print *,'after call EXCH,mype=',me
!$omp parallel do private(i,j)
do j = jsta, jend_m
do i = 1, im-1
DX( i,j) = ERAD*COS(GDLAT(I,J)*DTR)*(GDLON(I+1,J)-GDLON(I,J))*DTR
DY(i,j) = ERAD*(GDLAT(I,J)-GDLAT(I,J+1))*DTR ! like A*DPH
! F(I,J)=1.454441e-4*sin(gdlat(i,j)*DTR) ! 2*omeg*sin(phi)
end do
end do
if (me == 0) print*,'sample LATLON, DY, DY=',ii,jend, &
GDLAT(II,JEND),GDLON(II,JEND),DX(II,JEND),DY(II,JEND)
!$omp parallel do private(i,j)
do j=jsta,jend
do i=1,im
F(I,J) = 1.454441e-4*sin(gdlat(i,j)*DTR) ! 2*omeg*sin(phi)
end do
end do
impf = im
jmpf = jm
if (me == 0) print*,'impf,jmpf= ',impf,jmpf
!Moo print*,'impf,jmpf,nframe= ',impf,jmpf,nframe
! iyear=idate(4)+2000 ! older gfsio only has 2 digit year
iyear = idate(1)
imn = idate(2) ! ask Jun
iday = idate(3) ! ask Jun
ihrst = idate(4)
imin = idate(5)
jdate = 0
idate = 0
!
! read(startdate,15)iyear,imn,iday,ihrst,imin
15 format(i4,1x,i2,1x,i2,1x,i2,1x,i2)
if (me == 0) then
print*,'start yr mo day hr min =',iyear,imn,iday,ihrst,imin
print*,'processing yr mo day hr min=' &
,idat(3),idat(1),idat(2),idat(4),idat(5)
endif
!
idate(1) = iyear
idate(2) = imn
idate(3) = iday
idate(5) = ihrst
idate(6) = imin
SDAT(1) = imn
SDAT(2) = iday
SDAT(3) = iyear
jdate(1) = idat(3)
jdate(2) = idat(1)
jdate(3) = idat(2)
jdate(5) = idat(4)
jdate(6) = idat(5)
!
! CALL W3DIFDAT(JDATE,IDATE,2,RINC)
! ifhr=nint(rinc(2))
!
CALL W3DIFDAT(JDATE,IDATE,0,RINC)
!
ifhr = nint(rinc(2)+rinc(1)*24.)
ifmin = nint(rinc(3))
! if(ifhr /= nint(fhour))print*,'find wrong Grib file';stop
if (me == 0) then
print *,' idate=',idate
print *,' rinc=',rinc
print *,' ifhr=',ifhr
print*,' in INITPOST ifhr ifmin fileName=',ifhr,ifmin,fileName
! GFS has the same accumulation bucket for precipitation and fluxes and it is written to header
! the header has the start hour information so post uses it to recontruct bucket
tprec = 6.
if(ifhr>240)tprec=12.
tclod = tprec
trdlw = tprec
trdsw = tprec
tsrfc = tprec
tmaxmin = tprec
td3d = tprec
print*,'tprec = ',tprec
end if
call mpi_bcast(tprec,1,MPI_REAL,0,mpi_comm_comp,iret)
call mpi_bcast(tclod,1,MPI_REAL,0,mpi_comm_comp,iret)
call mpi_bcast(trdlw,1,MPI_REAL,0,mpi_comm_comp,iret)
call mpi_bcast(trdsw,1,MPI_REAL,0,mpi_comm_comp,iret)
call mpi_bcast(tsrfc,1,MPI_REAL,0,mpi_comm_comp,iret)
call mpi_bcast(tmaxmin,1,MPI_REAL,0,mpi_comm_comp,iret)
call mpi_bcast(td3d,1,MPI_REAL,0,mpi_comm_comp,iret)
! Getting tstart
tstart = 0.
print*,'tstart= ',tstart
! Getiing restart
RESTRT=.TRUE. ! set RESTRT as default
IF(tstart .GT. 1.0E-2)THEN
ifhr = ifhr+NINT(tstart)
rinc = 0
idate = 0
rinc(2) = -1.0*ifhr
call w3movdat(rinc,jdate,idate)
SDAT(1) = idate(2)
SDAT(2) = idate(3)
SDAT(3) = idate(1)
IHRST = idate(5)
print*,'new forecast hours for restrt run= ',ifhr
print*,'new start yr mo day hr min =',sdat(3),sdat(1) &
,sdat(2),ihrst,imin
END IF
imp_physics = 99 !set GFS mp physics to 99 for Zhao scheme
iCU_PHYSICS = 4
print*,'MP_PHYSICS=,cu_physics=',imp_physics,icu_physics
! Initializes constants for Ferrier microphysics
if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)then
CALL MICROINIT(imp_physics)
end if
! waiting to retrieve lat lon infor from raw GFS output
! VarName='DX'
! VarName='DY'
! GFS does not need DT to compute accumulated fields, set it to one
! VarName='DT'
DT = 1
! GFS does not need truelat
! VarName='TRUELAT1'
! VarName='TRUELAT2'
! Specify maptype=4 for Gaussian grid
! maptype=4
! write(6,*) 'maptype is ', maptype
! HBM2 is most likely not in Grib message, set them to ones
HBM2 = 1.0
! try to get kgds from flux grib file and then convert to igds that is used by GRIBIT.f
if(me == 0)then
jpds = -1.0
jgds = -1.0
igds = 0
call getgb(iunit,0,im_jm,0,jpds,jgds,kf &
,k,kpds,kgds,lb,dummy,ierr)
if(ierr == 0)then
call R63W72(KPDS,KGDS,JPDS,IGDS(1:18))
print*,'use IGDS from flux file for GFS= ',(IGDS(I),I=1,18)
else
print*,'no flux file, fill part of kgds with sigma file info'
kgds(1) = idrt
kgds(2) = im
kgds(3) = jm
end if
end if
call mpi_bcast(igds(1),18,MPI_INTEGER,0,mpi_comm_comp,iret)
call mpi_bcast(kgds(1),18,MPI_INTEGER,0,mpi_comm_comp,iret)
print*,'IGDS for GFS= ',(IGDS(I),I=1,18)
! Specigy grid type
! if(iostatusFlux==0)then
if(IGDS(4) /= 0)then
maptype = IGDS(3)
else if((im/2+1) == jm)then
maptype = 0 !latlon grid
else
maptype = 4 ! default gaussian grid
end if
gridtype = 'A'
write(6,*) 'maptype and gridtype is ', maptype,gridtype
if(idrt /= maptype)then
print*,'flux file and sigma file are on different grids - ', &
'post processing isterminated'
call mpi_abort()
stop
end if
levs = sighead%levs
ntrac = sighead%ntrac
ncld = sighead%ncldt
idvt = sighead%idvt
jcap = sighead%jcap
lnt2 = (jcap+1)*(jcap+2)
ntoz = mod(idvt,10) + 1
!jw ntcw = idvt/10 + 1
if( idvt/100 > 0 ) then
ntcw = 3
else
ntcw = idvt/10 + 1
endif
! start reading sigma file
! decompose l to read different levs with different pes
call mptgen(me,num_procs,1,1,lm,lsta,lend,kxa,kma,kna)
write(0,*)' me=',me,' lsta=',lsta,' lend=',lend,' kxa=',kxa, &
' kma=',kma,' kna=',kna
ii = im/2
jj = (jsta+jend)/2
print*,'lsta, lend= ',lsta,lend
if (lsta > lend) lend = lsta
allocate(dummy_h(im,jm),dummy_p(im,jm),dummy_px(im,jm),dummy_py(im,jm) &
,dummy_t(im,jm,lsta:lend),dummy_u(im,jm,lsta:lend) &
,dummy_v(im,jm,lsta:lend),dummy_d(im,jm,lsta:lend) &
,dummy_trc(im,jm,lsta:lend,ntrac))
! ,dummy12(im,jm,lsta:lend), &
! dummy13(im,jm,lsta:lend),dummy14(im,jm,lsta:lend), &
! dummy18(im,jm,lsta:lend),dummy19(im,jm,lsta:lend) )
if (me == 0) then
print*,'calling rtsig with lusig,lsta,lend,im_jm,kgds= ', &
lusig,lsta,lend,im_jm,kgds(1:20)
write(0,*)' levs=',levs,' ntrac=',ntrac,' ncld=',ncld,' jcap=',jcap &
,' lnt2=',lnt2,' ntoz=',ntoz,' ntcw=',ntcw
endif
call rtsig(lusig,sighead,lsta,lend,kgds,im_jm, & ! input
levs,ntrac,jcap,lnt2,me, & ! input
dummy_h,dummy_p,dummy_px,dummy_py, & ! output
dummy_t,dummy_u,dummy_v,dummy_d, & ! output
dummy_trc,iret) ! output
write(0,*)'aft rtsig,iret=',iret
if(iret /= 0)then
print*,'error reading sigma file, stopping'
print*,'error massage is ',iret
call mpi_abort()
end if
if(Debugprint)print*,'done with rtsig, sample t,u,v,q,cwm= ', &
dummy_t(1,1,lsta:lend), dummy_u(1,1,lsta:lend), &
dummy_v(1,1,lsta:lend), dummy_trc(1,1,lsta:lend,1), &
dummy_trc(1,1,lsta:lend,3),' lsta=',lsta,'lend=',lend
! write(0,*)'bf allocate '
! set threads for rest of the code
! call getenv('POST_THREADS',ENVAR)
! read(ENVAR, '(I2)')idum
! idum = max(idum+0,1)
! write(0,*)' post_threads=', idum
! if (idum > 0 .and. idum <= 32) then
! call OMP_SET_NUM_THREADS(idum)
! endif
! scatter to pes
allocate(dummy15(im,jsta_2l:jend_2u), &
dummy16(im,jsta_2l:jend_2u),dummy17(im,jsta_2l:jend_2u,lm))
! write(0,*)'af allocate '
! call mpi_scatterv(dummy_h(1,1),icnt,idsp,mpi_real &
! ,zint(1,jsta,lp1),icnt(me),mpi_real,0,MPI_COMM_COMP,iret)
! call mpi_scatterv(dummy_p(1,1),icnt,idsp,mpi_real &
! ,pint(1,jsta,lp1),icnt(me),mpi_real,0,MPI_COMM_COMP,iret)
! call mpi_scatterv(dummy_px(1,1),icnt,idsp,mpi_real &
! ,dummy15(1,jsta),icnt(me),mpi_real,0,MPI_COMM_COMP,iret)
! call mpi_scatterv(dummy_py(1,1),icnt,idsp,mpi_real &
! ,dummy16(1,jsta),icnt(me),mpi_real,0,MPI_COMM_COMP,iret)
!$omp parallel do private(i,j)
do j=jsta, jend
do i=1, im
zint(i,j,lp1) = dummy_h(i,j)
pint(i,j,lp1) = dummy_p(i,j)
dummy15(i,j) = dummy_px(i,j)
dummy16(i,j) = dummy_py(i,j)
enddo
enddo
deallocate (dummy_h,dummy_p,dummy_px,dummy_py)
write(0,*)'one scattering zs and ps',zint(ii,jj,lp1),pint(ii,jj,lp1)
ijmc = (jm-1)/num_procs+1
ijxc = jend-jsta+1
allocate (ta(im,ijmc,kma,num_procs))
allocate (tb(im,ijmc,kma,num_procs))
allocate (buf3d(im,lm,jsta:jend))
! write(0,*)'be mptranr4'
if(ijxc > ijmc) print*,'ijxc larger than ijmc =',ijxc,ijmc
call mptranr4(MPI_COMM_COMP,num_procs,im,im,im, &
ijmc,jm,ijxc,jm,kma,kxa,lm,lm,dummy_t,buf3d,ta,tb)
! write(0,*)'be set buf3d'
!$omp parallel do private(i,j,l,ll)
do l=1, lm
ll = lm-l+1
do j=jsta, jend
do i=1, im
T(i,j,l) = buf3d(i,ll,j)
end do
end do
end do
if (debugprint) then
do l=1, lm
print*,'sample i,j,l, T = ',ii,jj,l,t(ii,jj,l)
enddo
endif
! write(0,*)'be set uh'
call mptranr4(MPI_COMM_COMP,num_procs,im,im,im, &
ijmc,jm,ijxc,jm,kma,kxa,lm,lm,dummy_u,buf3d,ta,tb)
!$omp parallel do private(i,j,l,ll)
do l=1, lm
ll = lm-l+1
do j=jsta, jend
do i=1, im
uh(i,j,l) = buf3d(i,ll,j)
end do
end do
end do
if (debugprint) then
do l=1, lm
print*,'sample i,j,l,U = ',ii,jj,l,uh(ii,jj,l)
enddo
endif
! write(0,*)'be set vh'
call mptranr4(MPI_COMM_COMP,num_procs,im,im,im, &
ijmc,jm,ijxc,jm,kma,kxa,lm,lm,dummy_v,buf3d,ta,tb)
!$omp parallel do private(i,j,l,ll)
do l=1, lm
ll = lm-l+1
do j=jsta, jend
do i=1, im
vh(i,j,l) = buf3d(i,ll,j)
end do
end do
end do
if (debugprint) then
do l=1, lm
print*,'sample i,j,l,V = ',ii,jj,l,vh(ii,jj,l)
enddo
endif
!
! ltrc = min(lsta,levs) ! For processors greater than levs
! write(0,*)'be set q'
call mptranr4(MPI_COMM_COMP,num_procs,im,im,im,ijmc,jm, &
ijxc,jm,kma,kxa,lm,lm,dummy_trc(1,1,lsta,1),buf3d,ta,tb)
!$omp parallel do private(i,j,l,ll)
do l=1, lm
ll = lm-l+1
do j=jsta, jend
do i=1, im
q(i,j,l) = buf3d(i,ll,j)
end do
end do
end do
if (debugprint) then
do l=1, lm
print*,'sample i,j,l,Q = ',ii,jj,l,q(ii,jj,l)
enddo
endif
! write(0,*)'be set o3'
call mptranr4(MPI_COMM_COMP,num_procs,im,im,im,ijmc,jm, &
ijxc,jm,kma,kxa,lm,lm,dummy_trc(1,1,lsta,ntoz),buf3d,ta,tb)
!$omp parallel do private(i,j,l,ll)
do l=1, lm
ll = lm-l+1
do j=jsta, jend
do i=1, im
o3(i,j,l) = buf3d (i,ll,j)
end do
end do
end do
if (debugprint) then
do l=1, lm
print*,'sample i,j,l,O3 = ',ii,jj,l,o3(ii,jj,l)
enddo
endif
! write(0,*)'be set cld,ntcw=',ntcw
call mptranr4(MPI_COMM_COMP,num_procs,im,im,im,ijmc,jm, &
ijxc,jm,kma,kxa,lm,lm,dummy_trc(1,1,lsta,ntcw),buf3d,ta,tb)
! write(0,*)'aft mptranr4 cwm'
!$omp parallel do private(i,j,l,ll)
do l=1, lm
ll = lm-l+1
do j=jsta, jend
do i=1, im
cwm(i,j,l ) = buf3d (i,ll,j)
end do
end do
end do
if (debugprint) then
do l=1, lm
print*,'sample i,j,l,CWM = ',ii,jj,l,cwm(ii,jj,l)
enddo
endif
! write(0,*)' cwm=',cwm(1,36,100:150)*1000
! write(0,*)'be set div'
call mptranr4(MPI_COMM_COMP,num_procs,im,im,im, &
ijmc,jm,ijxc,jm,kma,kxa,lm,lm,dummy_d,buf3d,ta,tb)
!$omp parallel do private(i,j,l)
do l=1, lm
do j=jsta, jend
do i=1, im
dummy17(i,j,l ) = buf3d(i,l,j)
end do
end do
end do
if (debugprint) then
do l=1, lm
print*,'sample i,j,l,DIV = ',ii,jj,l,dummy17(ii,jj,l)
enddo
endif
deallocate(dummy_t,dummy_u,dummy_v,dummy_d,dummy_trc,ta,tb,buf3d)
! dummy18,dummy19)
!$omp parallel do private(i,j,l)
do l=1,lm
! if(debugprint)print*,'sample T Q U V,CWM',l,VarName,' = ',l,t(ii,jj,l),&
! q(ii,jj,l),u(ii,jj,l),v(ii,jj,l),cwm(ii,jj,l)
do j=jsta,jend
do i=1,im
if(t(i,j,l) < (TFRZ-15.) ) then ! separating cloud water from ice
qqi(i,j,l) = cwm(i,j,l)
else
qqw(i,j,l) = cwm(i,j,l)
end if
end do
end do
end do ! for l loop
! write(0,*)'be set qqi'
! write(0,*)' qqw=',qqw(1,36,100:150)
! write(0,*)' qqi=',qqi(1,36,100:150)
! compute model level pressure and omega
pdtop = spval
! pt = 0.
pt = 10000. ! this is for 100 hPa added by Moorthi
pd = spval ! GFS does not output PD
allocate (d2d(im,lm),u2d(im,lm),v2d(im,lm),pi2d(im,lm+1), &
pm2d(im,lm),omga2d(im,lm))
do j=jsta,jend
!$omp parallel do private(i,l,ll)
do l=1,lm
ll = lm-l+1
do i=1,im
u2d(i,l) = uh(i,j,ll) ! flipping u and v for calling modstuff
v2d(i,l) = vh(i,j,ll)
d2d(i,l) = dummy17(i,j,l)
end do
end do
call modstuff2(im,im,lm, &
sighead%idvc,sighead%idsl,sighead%nvcoord, &
sighead%vcoord,pint(1,j,lp1),dummy15(1,j), &
dummy16(1,j),d2d,u2d,v2d, &
pi2d,pm2d,omga2d,me)
!$omp parallel do private(i,l,ll)
do l=1,lm
ll = lm-l+1
do i=1,im
omga(i,j,l) = omga2d(i,ll)
pmid(i,j,l) = pm2d(i,ll)
pint(i,j,l) = pi2d(i,ll+1)
enddo
enddo
enddo ! end of j loop
! write(0,*)'be set pint'
! write(0,*)' PINT=',pint(ii,jj,:),' me=',me
deallocate(dummy15,dummy16,dummy17)
deallocate (d2d,u2d,v2d,pi2d,pm2d,omga2d)
allocate(wrk1(im,jsta:jend),wrk2(im,jsta:jend))
!$omp parallel do private(i,j)
do j=jsta,jend
do i=1,im
alpint(i,j,lp1) = log(pint(i,j,lp1))
wrk1(i,j) = log(PMID(I,J,LM))
wrk2(i,j) = T(I,J,LM)*(Q(I,J,LM)*con_fvirt+1.0)
fis(i,j) = zint(i,j,lp1)*con_G
FI(I,J,1) = FIS(I,J) &
+ wrk2(i,j)*con_rd*(ALPINT(I,J,Lp1)-wrk1(i,j))
ZMID(I,J,LM) = FI(I,J,1)/con_G
enddo
enddo
! write(0,*)'be set zmid'
! SECOND, INTEGRATE HEIGHT HYDROSTATICLY, GFS integrate height on mid-layer
ii = im/2
jj = (jsta+jend)/2
DO L=LM,2,-1 ! omit computing model top height because it's infinity
ll = l - 1
!$omp parallel do private(i,j,tvll,pmll,fact)
do j = jsta, jend
! write(0,*)' j=',j,' l=',l,' T=',T(1,j,l),' Q=',q(1,j,l), &
! ' pmid=',pmid(1,j,ll),pmid(1,j,l),' l=',l
do i = 1, im
! if (me == 40) &
! write(0,*)' i=',i,' j=',j,' l=',l,' T=',T(i,j,l),' Q=',q(i,j,l),&
! ' pmid=',pmid(i,j,ll),pmid(i,j,l),' l=',l
alpint(i,j,l) = log(pint(i,j,l))
tvll = T(I,J,LL)*(Q(I,J,LL)*con_fvirt+1.0)
pmll = log(PMID(I,J,LL))
FI(I,J,2) = FI(I,J,1) + (0.5*con_rd)*(wrk2(i,j)+tvll) &
* (wrk1(i,j)-pmll)
ZMID(I,J,LL) = FI(I,J,2)/con_G
!
FACT = (ALPINT(I,J,L)-wrk1(i,j)) / (pmll-wrk1(i,j))
ZINT(I,J,L) = ZMID(I,J,L) + (ZMID(I,J,LL)-ZMID(I,J,L)) * FACT
! if(i.eq.ii.and.j.eq.jj) &
! print*,'L,sample T,Q,ALPMID(L+1),ALPMID(L),ZMID= ' &
! ,l,T(I,J,L),Q(I,J,L),LOG(PMID(I,J,L+1)), &
! LOG(PMID(I,J,L)),ZMID(I,J,L)
FI(I,J,1) = FI(I,J,2)
wrk1(i,J) = pmll
wrk2(i,j) = tvll
ENDDO
ENDDO
if (me == 0) print*,'L ZINT= ',l,zint(ii,jj,l), &
'alpint=',ALPINT(ii,jj,l),'pmid=',LOG(PMID(Ii,Jj,L)),'pmid(l-1)=', &
LOG(PMID(Ii,Jj,L-1)),'zmd=',ZMID(Ii,Jj,L),'zmid(l-1)=',ZMID(Ii,Jj,L-1)
ENDDO
deallocate(wrk1,wrk2)
! write(0,*)'af sigma file'
! start retrieving data using getgb, first land/sea mask
Index = 50
VarName = avbl(index)
jpds = -1
jgds = -1
jpds(5) = iq(index)
jpds(6) = is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,sm)
where(sm /= spval) sm = 1.0 - sm ! convert to sea mask
! write(0,*)'get land-see mask'
if(debugprint)print*,'sample ',VarName,' = ',sm(im/2,(jsta+jend)/2)
! sea ice mask using getgb
Index = 51
VarName = avbl(index)
jpds = -1.0
jgds = -1.0
jpds(5) = iq(index)
jpds(6) = is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,sice)
where(sm/=spval .and. sm==0.0) sice = 0.0 !specify sea ice=0 at land
if(debugprint)print*,'sample ',VarName,' = ',sice(im/2,(jsta+jend)/2)
! Zhao scheme does not produce suspended rain and snow
qqr = 0.
qqs = 0.
! PBL height using getgb
Index = 221
VarName = avbl(index)
jpds = -1.0
jgds = -1.0
jpds(5) = iq(index)
jpds(6) = is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,pblh)
if(debugprint)print*,'sample ',VarName,' = ',pblh(im/2,(jsta+jend)/2)
! frictional velocity using getgb
Index=45
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,ustar)
if(debugprint)print*,'sample ',VarName,' = ',ustar(im/2,(jsta+jend)/2)
! roughness length using getgb
Index=44
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,z0)
if(debugprint)print*,'sample ',VarName,' = ',z0(im/2,(jsta+jend)/2)
! surface potential T using getgb
Index=26
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,ths)
! where(ths/=spval)ths=ths*(p1000/pint(:,:,lp1))**CAPA ! convert to THS
!$omp parallel do private(i,j)
do j=jsta,jend
do i=1,im
if (ths(i,j) /= spval) then
! write(0,*)' i=',i,' j=',j,' ths=',ths(i,j),' pint=',pint(i,j,lp1)
ths(i,j) = ths(i,j) * (p1000/pint(i,j,lp1))**capa
endif
enddo
enddo
if(debugprint)print*,'sample ',VarName,' = ',ths(im/2,(jsta+jend)/2)
QS = SPVAL ! GFS does not have surface specific humidity
twbs = SPVAL ! GFS does not have inst sensible heat flux
qwbs = SPVAL ! GFS does not have inst latent heat flux
! GFS does not have time step and physics time step, make up ones since they
! are not really used anyway
NPHS = 2.
DT = 80.
DTQ2 = DT * NPHS !MEB need to get physics DT
TSPH = 3600./DT !MEB need to get DT
! convective precip in m per physics time step using getgb
Index=272
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
! jpds(16)=3 ! CFSRR uses 1 for fhr>1532
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,avgcprate)
where(avgcprate /= spval)avgcprate=avgcprate*dtq2/1000. ! convert to m
if(debugprint)print*,'sample ',VarName,' = ',avgcprate(im/2,(jsta+jend)/2)
cprate=avgcprate
! construct tprec from flux grib massage
! comment this out because you can't get precip bucket from flux file
! if(me==0 .and. iostatusFlux==0)then
! if(kpds(16)==3)then ! Grib1 can't specify accumulated field fhr>1532
! if(KPDS(13)==1)then
! TPREC=float(KPDS(15)-KPDS(14))
! else if(KPDS(13)==10)then
! TPREC=float(KPDS(15)-KPDS(14))*3.0
! else if(KPDS(13)==11)then
! TPREC=float(KPDS(15)-KPDS(14))*6.0
! else if(KPDS(13)==12)then
! TPREC=float(KPDS(15)-KPDS(14))*12.0
! else if(KPDS(13)==2)then
! TPREC=float(KPDS(15)-KPDS(14))*24.0
! else
! TPREC=float(KPDS(15)-KPDS(14))
! end if
! else
! CALL GETENV('FHZER',ENVAR)
! read(ENVAR, '(I2)')idum
! tprec = idum * 1.0
! print*,'TPREC from FHZER= ',tprec
! end if
! end if
! call mpi_bcast(tprec,1,MPI_REAL,0,mpi_comm_comp,iret)
! precip rate in m per physics time step using getgb
Index=271
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,avgprec)
where(avgprec /= spval) avgprec = avgprec*dtq2/1000. ! convert to m
if(debugprint)print*,'sample ',VarName,' = ',avgprec(im/2,(jsta+jend)/2)
! inst precip rate in m per physics time step using sfcio
if(me==0)then
call getenv('SFCINPUT',sfcfilename)
print*,'opening sfcfile to read',sfcfilename
call sfcio_srohdc(35,sfcfilename,head,data,iret)
if(iret /=0 )then
print*,'fail to read ',sfcfilename
dummy = spval
dummy2 = spval
else
dummy = data%tprcp
print '(f8.2)',dummy(1,1)
! dummy2 = data%srflag
end if
end if
call mpi_scatterv(dummy(1,1),icnt,idsp,mpi_real &
,prec(1,jsta),icnt(me),mpi_real,0,MPI_COMM_COMP,iret)
print*,'sampe inst precip= ',prec(im/2,jsta)
where(prec /= spval) prec = prec*dtq2/1000. ! convert to m
! call mpi_scatterv(dummy2(1,1),icnt,idsp,mpi_real &
! ,sr(1,jsta),icnt(me),mpi_real,0,MPI_COMM_COMP,iret)
! print*,'sampe GFS sr= ',sr(im/2,jsta)
deallocate(dummy2)
! prec=avgprec !set avg cprate to inst one to derive other fields
! GFS does not have accumulated total, gridscale, and convective precip, will use inst precip to derive in SURFCE.f
lspa = spval ! GFS does not have similated precip
! inst snow water eqivalent using getgb
Index=119
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,sno)
if(debugprint)print*,'sample ',VarName,' = ',sno(im/2,(jsta+jend)/2)
! snow depth in mm using getgb
Index=224
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,si)
where(si /= spval) si = si*1000. ! convert to mm
if(debugprint)print*,'sample ',VarName,' = ',si(im/2,(jsta+jend)/2)
CLDEFI = SPVAL ! GFS does not have convective cloud efficiency
TH10 = SPVAL ! GFS does not have 10 m theta
Q10 = SPVAL ! GFS does not have 10 m humidity
! 2m T using nemsio
Index=106
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=2
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,tshltr)
if(debugprint)print*,'sample ',VarName,' = ',tshltr(im/2,(jsta+jend)/2)
! GFS does not have 2m pres, estimate it, also convert t to theta
!$omp parallel do private(i,j)
Do j=jsta,jend
Do i=1,im
if(tshltr(i,j) /= spval)then
! if (me == 127) write(0,*)' i=',i,' j=',j,' tshltr=',tshltr(i,j) &
! ,' pint=',pint(I,J,lm+1)
PSHLTR(I,J) = pint(I,J,lm+1)*EXP(-0.068283/tshltr(i,j))
tshltr(i,j) = tshltr(i,j)*(p1000/PSHLTR(I,J))**CAPA ! convert to theta
else
PSHLTR(I,J) = spval
end if
! if (j.eq.jm/2 .and. mod(i,50).eq.0)
! + print*,'sample 2m T and P after scatter= '
! + ,i,j,tshltr(i,j),pshltr(i,j)
end do
end do
! 2m specific humidity using gfsio
! VarName='spfh'
! VcoordName='2m above gnc'
! l=1
! if(me == 0)then
! call gfsio_readrecvw34(gfile,trim(VarName),trim(VcoordName) &
! + ,l,dummy,iret=iret)
! if (iret /= 0) then
! print*,VarName," not found in file-Assigned missing values"
! dummy=spval
! end if
! end if
! call mpi_scatterv(dummy,icnt,idsp,mpi_real &
! + ,qshltr(1,jsta),icnt(me),mpi_real,0,MPI_COMM_COMP,iret)
! if (iret /= 0)print*,'Error scattering array';stop
! 2m specific humidity using nemsio
Index=112
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,qshltr)
if(debugprint)print*,'sample ',VarName,' = ',qshltr(im/2,(jsta+jend)/2)
Q2 = SPVAL ! GFS does not have TKE because it uses GFS scheme
! GFS does not have surface exchange coeff
ALBASE = SPVAL ! GFS does not have snow free albedo
! mid day avg albedo in fraction using nemsio
Index=266
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,avgalbedo)
where(avgalbedo /= spval)avgalbedo=avgalbedo/100. ! convert to fraction
if(debugprint)print*,'sample ',VarName,' = ',avgalbedo(im/2,(jsta+jend)/2)
! time averaged column cloud fractionusing nemsio
Index=144
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,avgtcdc)
where(avgtcdc /= spval)avgtcdc=avgtcdc/100. ! convert to fraction
if(debugprint)print*,'sample ',VarName,' = ',avgtcdc(im/2,(jsta+jend)/2)
! if(me==0 .and. iostatusFlux==0)then
! if(KPDS(13)==1)then
! TCLOD=float(KPDS(15)-KPDS(14))
! else if(KPDS(13)==10)then
! TCLOD=float(KPDS(15)-KPDS(14))*3.0
! else if(KPDS(13)==11)then
! TCLOD=float(KPDS(15)-KPDS(14))*6.0
! else if(KPDS(13)==12)then
! TCLOD=float(KPDS(15)-KPDS(14))*12.0
! else if(KPDS(13)==2)then
! TCLOD=float(KPDS(15)-KPDS(14))*24.0
! else
! TCLOD=float(KPDS(15)-KPDS(14))
! end if
! end if
! call mpi_bcast(tclod,1,MPI_REAL,0,mpi_comm_comp,iret)
! print*,'TCLOD from flux grib massage= ',TCLOD
Czen = spval ! GFS probably does not use zenith angle (What???)
CZMEAN = SPVAL
! maximum snow albedo in fraction using nemsio
Index=227
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,mxsnal)
where(mxsnal /= spval)mxsnal=mxsnal/100. ! convert to fraction
if(debugprint)print*,'sample ',VarName,' = ',mxsnal(im/2,(jsta+jend)/2)
radot = spval ! GFS does not have inst surface outgoing longwave
! GFS probably does not use sigt4, set it to sig*t^4
!$omp parallel do private(i,j,tlmh)
Do j=jsta,jend
Do i=1,im
TLMH = T(I,J,LM)
Sigt4(I,j) = 5.67E-8*TLMH*TLMH*TLMH*TLMH
End do
End do
! TG is not used, skip it for now
! allocate(qstl(lm))
allocate(p2d(im,lm),t2d(im,lm),q2d(im,lm),cw2d(im,lm), &
qs2d(im,lm),cfr2d(im,lm))
do j=jsta,jend
!$omp parallel do private(i,k,es)
do k=1,lm
do i=1,im
p2d(i,k) = pmid(i,j,k)*0.01
t2d(i,k) = t(i,j,k)
q2d(i,k) = q(i,j,k)
cw2d(i,k) = cwm(i,j,k)
es = min(fpvsnew(t(i,j,k)),pmid(i,j,k))
qs2d(i,k) = con_eps*es/(pmid(i,j,k)+con_epsm1*es)!saturation q for GFS
enddo
enddo
call progcld1 &
!...................................
! --- inputs:
( p2d,t2d,q2d,qs2d,cw2d,im,lm,0, &
! --- outputs:
cfr2d &
)
!$omp parallel do private(i,k)
do k=1,lm
do i=1,im
cfr(i,j,k) = cfr2d(i,k)
enddo
end do
end do
deallocate(p2d,t2d,q2d,qs2d,cw2d,cfr2d)
! ask moorthi if there is snow rate in GFS
! varname='SR'
! call retrieve_index(index,VarName,varname_all,nrecs,iret)
! if (iret /= 0) then
! print*,VarName," not found in file-Assigned missing values"
! SR=SPVAL
! else
! this_offset=file_offset(index+1)+(jsta_2l-1)*4*im
! this_length=im*(jend_2u-jsta_2l+1)
! call mpi_file_read_at(iunit,this_offset
! + ,sr,this_length,mpi_real4
! + , mpi_status_ignore, ierr)
! if (ierr /= 0) then
! print*,"Error reading ", VarName,"Assigned missing values"
! SR=SPVAL
! end if
! end if
! GFS does not have inst cloud fraction for high, middle, and low cloud
cfrach = spval
cfracl = spval
cfracm = spval
! ave high cloud fraction using nemsio
Index=302
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,avgcfrach)
where(avgcfrach /= spval)avgcfrach=avgcfrach/100. ! convert to fraction
if(debugprint)print*,'sample ',VarName,' = ',avgcfrach(im/2,(jsta+jend)/2)
! ave low cloud fraction using nemsio
VarName='tcdc'
Index=300
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,avgcfracl)
where(avgcfracl /= spval)avgcfracl=avgcfracl/100. ! convert to fraction
if(debugprint)print*,'sample ',VarName,' = ',avgcfracl(im/2,(jsta+jend)/2)
! ave middle cloud fraction using nemsio
VarName='tcdc'
VcoordName='mid cld lay'
Index=301
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,avgcfracm)
where(avgcfracm /= spval)avgcfracm=avgcfracm/100. ! convert to fraction
if(debugprint)print*,'sample ',VarName,' = ',avgcfracm(im/2,(jsta+jend)/2)
! inst convective cloud fraction using nemsio
VarName='tcdc'
VcoordName='convect-cld laye'
Index=196
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=71
jpds(6)=244
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,cnvcfr)
where(cnvcfr /= spval)cnvcfr=cnvcfr/100. ! convert to fraction
if(debugprint)print*,'sample ',VarName,' = ',cnvcfr(im/2,(jsta+jend)/2)
! slope type using nemsio
Index=223
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,buf)
where(buf /= spval)islope=nint(buf)
if(debugprint)print*,'sample ',VarName,' = ',islope(im/2,(jsta+jend)/2)
! plant canopy sfc wtr in m using nemsio
VarName='cnwat'
VcoordName='sfc'
Index=118
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,cmc)
where(cmc /= spval)cmc=cmc/1000. ! convert from kg*m^2 to m
if(debugprint)print*,'sample ',VarName,' = ',cmc(im/2,(jsta+jend)/2)
grnflx = spval ! GFS does not have inst ground heat flux
! GFS does not have snow cover yet
! VarName='gflux'
! VcoordName='sfc'
! l=1
! if(me == 0)then
! call gfsio_readrecvw34(gfile,trim(VarName),trim(VcoordName) &
! + ,l,dummy,iret=iret)
! if (iret /= 0) then
! print*,VarName," not found in file-Assigned missing values"
! dummy=spval
! end if
! end if
! call mpi_scatterv(dummy,icnt,idsp,mpi_real &
! + , pctsno(1,jsta),icnt(me),mpi_real,0,MPI_COMM_COMP,iret)
! if (iret /= 0)print*,'Error scattering array';stop
soiltb = spval
tg = spval
! vegetation fraction in fraction. using nemsio
VarName='veg'
VcoordName='sfc'
Index=170
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,vegfrc)
where(vegfrc /= spval)
vegfrc = vegfrc/100. ! convert to fraction
elsewhere (vegfrc == spval)
vegfrc = 0. ! set to zero to be reasonable input for crtm
end where
if(debugprint)print*,'sample ',VarName,' = ',vegfrc(im/2,(jsta+jend)/2)
! GFS doesn not yet output soil layer thickness, assign SLDPTH to be the same as nam
SLDPTH(1) = 0.10
SLDPTH(2) = 0.3
SLDPTH(3) = 0.6
SLDPTH(4) = 1.0
! liquid volumetric soil mpisture in fraction using nemsio
VarName='soill'
VcoordName='0-10 cm down'
Index=225
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
do l=1,nsoil
if(l == 1)then
jpds(7)=nint(sldpth(1)*100.)
else
soilayert=0
do n=1,l-1
soilayert=soilayert+sldpth(n)*100.
end do
soilayerb=soilayert+sldpth(l)*100.
jpds(7)=nint(soilayert*256.+soilayerb)
end if
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,sh2o(1,jsta_2l,l))
if(debugprint)print*,'sample l',VarName,' = ',l,sh2o(im/2,(jsta+jend)/2,1)
End do ! do loop for l
! volumetric soil moisture using nemsio
VarName='soilw'
VcoordName='0-10 cm down'
Index=117
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
do l=1,nsoil
if(l == 1)then
jpds(7)=nint(sldpth(1)*100.)
else
soilayert=0
do n=1,l-1
soilayert=soilayert+sldpth(n)*100.
end do
soilayerb=soilayert+sldpth(l)*100.
jpds(7)=nint(soilayert*256.+soilayerb)
end if
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,smc(1,jsta_2l,l))
if(debugprint)print*,'sample l',VarName,' = ',l,smc(im/2,(jsta+jend)/2,1)
End do ! do loop for l
! soil temperature using nemsio
VarName='tmp'
VcoordName='0-10 cm down'
Index=116
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=11 ! GFS used 11 for soil T instead of 85
jpds(6)=is(index)
do l=1,nsoil
if(l == 1)then
jpds(7)=nint(sldpth(1)*100.)
else
soilayert=0
do n=1,l-1
soilayert=soilayert+sldpth(n)*100.
end do
soilayerb=soilayert+sldpth(l)*100.
jpds(7)=nint(soilayert*256.+soilayerb)
end if
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,stc(1,jsta_2l,l))
if(debugprint)print*,'sample l','stc',' = ',1,stc(im/2,(jsta+jend)/2,1)
End do ! do loop for l
! GFS does not output time averaged convective and strat cloud fraction, set acfrcv to spval, ncfrcv to 1
acfrcv = spval
ncfrcv = 1.0
! GFS does not output time averaged cloud fraction, set acfrst to spval, ncfrst to 1
acfrst = spval
ncfrst = 1.0
! GFS does not have storm runoff
ssroff = spval
! GFS does not have UNDERGROUND RUNOFF
bgroff = spval
! GFS incoming sfc longwave has been averaged over 6 hr bucket, set ARDLW to 1
ardlw = 1.0
! trdlw=6.0
! GFS does not have inst incoming sfc longwave
rlwin = spval
! GFS does not have inst model top outgoing longwave
rlwtoa = spval
! time averaged incoming sfc longwave using nemsio
VarName='dlwrf'
VcoordName='sfc'
Index=127
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,alwin)
! if(me==0 .and. iostatusFlux==0)then
! if(KPDS(13)==1)then
! TRDLW=float(KPDS(15)-KPDS(14))
! else if(KPDS(13)==10)then
! TRDLW=float(KPDS(15)-KPDS(14))*3.0
! else if(KPDS(13)==11)then
! TRDLW=float(KPDS(15)-KPDS(14))*6.0
! else if(KPDS(13)==12)then
! TRDLW=float(KPDS(15)-KPDS(14))*12.0
! else if(KPDS(13)==2)then
! TRDLW=float(KPDS(15)-KPDS(14))*24.0
! else
! TRDLW=float(KPDS(15)-KPDS(14))
! end if
! end if
! call mpi_bcast(TRDLW,1,MPI_REAL,0,mpi_comm_comp,iret)
! print*,'TRDLW from flux grib massage= ',TRDLW
! time averaged outgoing sfc longwave using gfsio
VarName='ulwrf'
VcoordName='sfc'
Index=129
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,alwout)
where(alwout /= spval) alwout=-alwout ! CLDRAD puts a minus sign before gribbing
if(debugprint)print*,'sample l',VarName,' = ',1,alwout(im/2,(jsta+jend)/2)
! time averaged outgoing model top longwave using gfsio
VarName='ulwrf'
VcoordName='nom. top'
Index=131
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,alwtoa)
if(debugprint)print*,'sample l',VarName,' = ',1,alwtoa(im/2,(jsta+jend)/2)
! GFS does not have inst incoming sfc shortwave
rswin = spval
! GFS does not have inst incoming clear sky sfc shortwave
rswinc = spval
! GFS does not have inst outgoing sfc shortwave
rswout = spval
! GFS incoming sfc longwave has been averaged, set ARDLW to 1
ardsw = 1.0
! trdsw=6.0
! time averaged incoming sfc shortwave using gfsio
Index=126
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,aswin)
if(debugprint)print*,'sample l',VarName,' = ',1,aswin(im/2,(jsta+jend)/2)
! if(me==0 .and. iostatusFlux==0)then
! if(KPDS(13)==1)then
! TRDSW=float(KPDS(15)-KPDS(14))
! else if(KPDS(13)==10)then
! TRDSW=float(KPDS(15)-KPDS(14))*3.0
! else if(KPDS(13)==11)then
! TRDSW=float(KPDS(15)-KPDS(14))*6.0
! else if(KPDS(13)==12)then
! TRDSW=float(KPDS(15)-KPDS(14))*12.0
! else if(KPDS(13)==2)then
! TRDSW=float(KPDS(15)-KPDS(14))*24.0
! else
! TRDSW=float(KPDS(15)-KPDS(14))
! end if
! end if
! call mpi_bcast(trdsw,1,MPI_REAL,0,mpi_comm_comp,iret)
! print*,'TRDSW from flux grib massage= ',trdsw
! time averaged incoming sfc uv-b using getgb
VarName='duvb'
VcoordName='sfc'
Index=298
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,auvbin)
if(debugprint)print*,'sample l',VarName,' = ',1,auvbin(im/2,(jsta+jend)/2)
! time averaged incoming sfc clear sky uv-b using getgb
VarName='cduvb'
VcoordName='sfc'
Index=297
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,auvbinc)
if(debugprint)print*,'sample l',VarName,' = ',1,auvbinc(im/2,(jsta+jend)/2)
! time averaged outgoing sfc shortwave using gfsio
VarName='uswrf'
VcoordName='sfc'
Index=128
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,aswout)
where(aswout /= spval) aswout=-aswout ! CLDRAD puts a minus sign before gribbing
if(debugprint)print*,'sample l',VarName,' = ',1,aswout(im/2,(jsta+jend)/2)
! time averaged model top outgoing shortwave
VarName='uswrf'
VcoordName='nom. top'
Index=130
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,aswtoa)
if(debugprint)print*,'sample l',VarName,' = ',1,aswtoa(im/2,(jsta+jend)/2)
! time averaged incoming clear sky sfc shortwave using getgb
Index=383
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
jpds(13)=3
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,aswinc)
! time averaged outgoing clear sky sfc shortwave using getgb
Index=386
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
jpds(13)=3
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,aswoutc)
! time averaged outgoing clear sky toa shortwave using getgb
Index=387
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
jpds(13)=3
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,aswtoac)
! time averaged model top incoming shortwave
Index=388
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,aswintoa)
! time averaged surface visible beam downward solar flux
Index=401
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,avisbeamswin)
! time averaged surface visible diffuse downward solar flux
Index=402
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,avisdiffswin)
! time averaged surface near IR beam downward solar flux
Index=403
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,airbeamswin)
! time averaged surface near IR diffuse downward solar flux
Index=404
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,airdiffswin)
! time averaged surface sensible heat flux, multiplied by -1 because wrf model flux
! has reversed sign convention using gfsio
VarName='shtfl'
VcoordName='sfc'
Index=43
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,sfcshx)
where (sfcshx /= spval)sfcshx=-sfcshx
if(debugprint)print*,'sample l',VarName,' = ',1,sfcshx(im/2,(jsta+jend)/2)
! if(me==0 .and. iostatusFlux==0)then
! if(KPDS(13)==1)then
! TSRFC=float(KPDS(15)-KPDS(14))
! else if(KPDS(13)==10)then
! TSRFC=float(KPDS(15)-KPDS(14))*3.0
! else if(KPDS(13)==11)then
! TSRFC=float(KPDS(15)-KPDS(14))*6.0
! else if(KPDS(13)==12)then
! TSRFC=float(KPDS(15)-KPDS(14))*12.0
! else if(KPDS(13)==2)then
! TSRFC=float(KPDS(15)-KPDS(14))*24.0
! else
! TSRFC=float(KPDS(15)-KPDS(14))
! end if
! end if
! call mpi_bcast(tsrfc,1,MPI_REAL,0,mpi_comm_comp,iret)
! print*,'TSRFC from flux grib massage= ',tsrfc
! GFS surface flux has been averaged, set ASRFC to 1
asrfc=1.0
! tsrfc=6.0
! time averaged surface latent heat flux, multiplied by -1 because wrf model flux
! has reversed sign vonvention using gfsio
Index=42
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,sfclhx)
where (sfclhx /= spval)sfclhx=-sfclhx
if(debugprint)print*,'sample l',VarName,' = ',1,sfclhx(im/2,(jsta+jend)/2)
! time averaged ground heat flux using nemsio
VarName='gflux'
VcoordName='sfc'
Index=135
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,subshx)
if(debugprint)print*,'sample l',VarName,' = ',1,subshx(im/2,(jsta+jend)/2)
! GFS does not have snow phase change heat flux
snopcx=spval
! time averaged zonal momentum flux using gfsio
VarName='uflx'
VcoordName='sfc'
Index=269
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,sfcux)
if(debugprint)print*,'sample l',VarName,' = ',1,sfcux(im/2,(jsta+jend)/2)
! time averaged meridional momentum flux using nemsio
VarName='vflx'
VcoordName='sfc'
Index=270
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,sfcvx)
if(debugprint)print*,'sample l',VarName,' = ',1,sfcvx(im/2,(jsta+jend)/2)
! GFS does not use total momentum flux
sfcuvx=spval
! time averaged zonal gravity wave stress using nemsio
VarName='u-gwd'
VcoordName='sfc'
Index=315
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,sfcugs)
gtaux=sfcugs
if(debugprint)print*,'sample l',VarName,' = ',1,gtaux(im/2,(jsta+jend)/2)
! time averaged meridional gravity wave stress using getgb
VarName='v-gwd'
VcoordName='sfc'
Index=316
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,sfcvgs)
gtauy=sfcvgs
if(debugprint)print*,'sample l',VarName,' = ',1,gtauy(im/2,(jsta+jend)/2)
! time averaged accumulated potential evaporation
VarName='pevpr'
VcoordName='sfc'
Index=242
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,potevp)
if(debugprint)print*,'sample l',VarName,' = ',1,potevp(im/2,(jsta+jend)/2)
! GFS does not have temperature tendency due to long wave radiation
rlwtt=spval
! GFS does not have temperature tendency due to solar radiation
rswtt=spval
! GFS does not have temperature tendency due to latent heating from convection
tcucn=spval
tcucns=spval
! set avrain to 1
avrain=1.0
avcnvc=1.0
theat=6.0 ! just in case GFS decides to output T tendency
! GFS does not have temperature tendency due to latent heating from grid scale
train=spval
! 10 m u using nemsio
VarName='ugrd'
VcoordName='10 m above gnd'
Index=64
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=10
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,u10)
if(debugprint)print*,'sample l',VarName,' = ',1,u10(im/2,(jsta+jend)/2)
u10h=u10
! 10 m v using gfsio
VarName='vgrd'
VcoordName='10 m above gnd'
Index=65
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=10
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,v10)
if(debugprint)print*,'sample l',VarName,' = ',1,v10(im/2,(jsta+jend)/2)
v10h=v10
! GFS does not have soil moisture availability
smstav=spval
! GFS does not have total soil moisture
smstot=spval
! vegetation type, it's in GFS surface file, hopefully will merge into gfsio soon
VarName='vgtyp'
VcoordName='sfc'
Index=218
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,buf)
where (buf /= spval)
ivgtyp=nint(buf)
elsewhere
ivgtyp=0 !need to feed reasonable value to crtm
end where
if(debugprint)print*,'sample l',VarName,' = ',1,ivgtyp(im/2,(jsta+jend)/2)
! soil type, it's in GFS surface file, hopefully will merge into gfsio soon
VarName='sotyp'
VcoordName='sfc'
Index=219
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,buf)
where (buf /= spval)
isltyp=nint(buf)
elsewhere
isltyp=0 !need to feed reasonable value to crtm
end where
if(debugprint)print*,'sample l',VarName,' = ',1,isltyp(im/2,(jsta+jend)/2)
! GFS does not have accumulated surface evaporation
sfcevp=spval
! GFS does not have surface exchange coeefficient
sfcexc=spval
! GFS does not have averaged accumulated snow
acsnow=spval
! GFS does not have snow melt
acsnom=spval
! GFS does not have sst????
sst=spval
! GFS does not have mixing length
EL_PBL=spval
! GFS does not output exchange coefficient
exch_h=spval
! GFS does not have THZ0, use THS to substitute
thz0=ths
if(debugprint)print*,'sample l',VarName,' = ',1,thz0(im/2,(jsta+jend)/2)
! GFS does not output humidity at roughness length
qz0=spval
! GFS does not output u at roughness length
uz0=spval
! GFS does not output humidity at roughness length
vz0=spval
! retrieve inst convective cloud top, GFS has cloud top pressure instead of index,
! will need to modify CLDRAD.f to use pressure directly instead of index
VarName='pres'
VcoordName='convect-cld top'
Index=189
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,ptop)
if(debugprint)print*,'sample l',VarName,' = ',1,ptop(im/2,(jsta+jend)/2)
htop = spval
do j=jsta,jend
do i=1,im
if(ptop(i,j) <= 0.0) ptop(i,j) = spval
if(ptop(i,j) < spval) then
do l=1,lm
if(ptop(i,j) <= pmid(i,j,l))then
htop(i,j) = l
exit
end if
end do
end if
end do
end do
if (me == 0) then
l = lm/2
print*,'sample ptop,pmid pmid-1,pint= ', &
ptop(ii,jj),pmid(ii,jj,l),pmid(ii,jj,l-1),pint(ii,jj,l),htop(ii,jj)
endif
! retrieve inst convective cloud bottom, GFS has cloud top pressure instead of index,
! will need to modify CLDRAD.f to use pressure directly instead of index
VarName='pres'
VcoordName='convect-cld bot'
Index=188
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,pbot)
if(debugprint)print*,'sample l',VarName,VcoordName,' = ',1,pbot(im/2,(jsta+jend)/2)
hbot = spval
do j=jsta,jend
do i=1,im
if(pbot(i,j) <= 0.0) pbot(i,j) = spval
! if(.not.lb(i,j))print*,'false bitmask for pbot at '
! + ,i,j,pbot(i,j)
if(pbot(i,j) .lt. spval)then
do l=lm,1,-1
if(pbot(i,j) >= pmid(i,j,l))then
hbot(i,j) = l
exit
end if
end do
end if
end do
end do
print*,'sample pbot,pmid= ', pbot(ii,jj),pmid(ii,jj,l),hbot(ii,jj)
! retrieve time averaged low cloud top pressure using nemsio
VarName='pres'
VcoordName='low cld top'
Index=304
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,ptopl)
if(debugprint)print*,'sample l',VarName,' = ',1,ptopl(im/2,(jsta+jend)/2)
! retrieve time averaged low cloud bottom pressure using nemsio
VarName='pres'
VcoordName='low cld bot'
Index=303
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,pbotl)
if(debugprint)print*,'sample l',VarName,' = ',1,pbotl(im/2,(jsta+jend)/2)
! retrieve time averaged low cloud top temperature using nemsio
VarName='tmp'
VcoordName='low cld top'
Index=305
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,Ttopl)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,Ttopl(im/2,(jsta+jend)/2)
! retrieve time averaged middle cloud top pressure using nemsio
Index=307
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,ptopm)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,ptopm(im/2,(jsta+jend)/2)
! retrieve time averaged middle cloud bottom pressure using nemsio
VarName='pres'
VcoordName='mid cld bot'
Index=306
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,pbotm)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,pbotm(im/2,(jsta+jend)/2)
! retrieve time averaged middle cloud top temperature using nemsio
VarName='tmp'
VcoordName='mid cld top'
Index=308
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,Ttopm)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,Ttopm(im/2,(jsta+jend)/2)
! retrieve time averaged high cloud top pressure using nemsio *********
VarName='pres'
VcoordName='high cld top'
Index=310
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,ptoph)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,ptoph(im/2,(jsta+jend)/2)
! retrieve time averaged high cloud bottom pressure using nemsio
VarName='pres'
VcoordName='high cld bot'
Index=309
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,pboth)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,pboth(im/2,(jsta+jend)/2)
! retrieve time averaged high cloud top temperature using nemsio
VarName='tmp'
VcoordName='high cld top'
Index=311
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,Ttoph)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,Ttoph(im/2,(jsta+jend)/2)
! retrieve boundary layer cloud cover using nemsio
VarName='tcdc'
VcoordName='bndary-layer cld'
Index=342
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,pblcfr)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,pblcfr(im/2,(jsta+jend)/2)
where (pblcfr /= spval)pblcfr=pblcfr/100. ! convert to fraction
! retrieve cloud work function using nemsio
Index=313
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,cldwork)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,cldwork(im/2,(jsta+jend)/2)
! retrieve water runoff using nemsio
VarName='watr'
VcoordName='sfc'
Index=343
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=0
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,runoff)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,runoff(im/2,(jsta+jend)/2)
! retrieve shelter max temperature using nemsio
VarName='tmax'
VcoordName='2 m above gnd'
Index=345
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=2
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,maxtshltr)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,maxtshltr(im/2,(jsta+jend)/2)
! retrieve shelter max temperature using nemsio
Index=346
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
jpds(7)=2
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,mintshltr)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,mintshltr(im/2,(jsta+jend)/2)
MAXRHSHLTR=SPVAL
MINRHSHLTR=SPVAL
! bucket for max and min temperature and RH
! if(me==0 .and. iostatusFlux==0)then
! if(KPDS(13)==1)then
! TMAXMIN=float(KPDS(15)-KPDS(14))
! else if(KPDS(13)==10)then
! TMAXMIN=float(KPDS(15)-KPDS(14))*3.0
! else if(KPDS(13)==11)then
! TMAXMIN=float(KPDS(15)-KPDS(14))*6.0
! else if(KPDS(13)==12)then
! TMAXMIN=float(KPDS(15)-KPDS(14))*12.0
! else if(KPDS(13)==2)then
! TMAXMIN=float(KPDS(15)-KPDS(14))*24.0
! else
! TMAXMIN=float(KPDS(15)-KPDS(14))
! end if
! end if
! call mpi_bcast(TMAXMIN,1,MPI_REAL,0,mpi_comm_comp,iret)
! print*,'TMAXMIN from flux grib massage= ',TMAXMIN
! retrieve ice thickness using nemsio
VarName='icetk'
VcoordName='sfc'
Index=349
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,dzice)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,dzice(im/2,(jsta+jend)/2)
! retrieve wilting point using nemsio
VarName='wilt'
VcoordName='sfc'
Index=236
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,smcwlt)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,smcwlt(im/2,(jsta+jend)/2)
! retrieve sunshine duration using nemsio
VarName='sunsd'
VcoordName='sfc'
Index=396
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,suntime)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,suntime(im/2,(jsta+jend)/2)
! retrieve field capacity using nemsio
VarName='fldcp'
VcoordName='sfc'
Index=397
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,fieldcapa)
if(debugprint)print*,'sample l',VcoordName,VarName,' = ', &
1,fieldcapa(im/2,(jsta+jend)/2)
! retrieve snowfall rate using getgb
Index=405
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,snowfall)
! retrieve frozen precipitation fraction using getgb
Index=172
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
call getgbandscatter(me,iunit,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,sr)
print*,'sampe GFS sr= ',sr(im/2,jsta)
! GFS does not have deep convective cloud top and bottom fields
HTOPD = SPVAL
HBOTD = SPVAL
HTOPS = SPVAL
HBOTS = SPVAL
CUPPT = SPVAL
!!!! DONE GETTING
! Will derive isobaric OMEGA from continuity equation later.
! OMGA=SPVAL
! retrieve d3d fields if it's listed
if (me == 0) print*,'iostatus for d3d file= ',iostatusD3D
if(iostatusD3D == 0) then ! start reading d3d file
! retrieve longwave tendency using getgb
Index=41
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
do l=1,lm
jpds(7)=l
ll = lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,rlwtt(1,jsta_2l,ll))
end do
! bucket for max and min temperature and RH
! if(me==0 .and. iostatusFlux==0)then
! if(KPDS(13)==1)then
! TD3D=float(KPDS(15)-KPDS(14))
! else if(KPDS(13)==10)then
! TD3D=float(KPDS(15)-KPDS(14))*3.0
! else if(KPDS(13)==11)then
! TD3D=float(KPDS(15)-KPDS(14))*6.0
! else if(KPDS(13)==12)then
! TD3D=float(KPDS(15)-KPDS(14))*12.0
! else if(KPDS(13)==2)then
! TD3D=float(KPDS(15)-KPDS(14))*24.0
! else
! TD3D=float(KPDS(15)-KPDS(14))
! end if
! end if
! call mpi_bcast(TD3D,1,MPI_REAL,0,mpi_comm_comp,iret)
! print*,'TD3D from D3D grib massage= ',TD3D
! retrieve shortwave tendency using getgb
Index=40
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
do l=1,lm
jpds(7)=l
ll = lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,rswtt(1,jsta_2l,ll))
end do
! retrieve vertical diffusion tendency using getgb
Index=356
VarName='VDIFF TNDY'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll = lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,vdifftt(1,jsta_2l,ll))
end do
! retrieve deep convective tendency using getgb
Index=79
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,tcucn(1,jsta_2l,ll))
end do
! retrieve shallow convective tendency using getgb
Index=358
VarName='S CNVCT TNDY'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,tcucns(1,jsta_2l,ll))
end do
! retrieve grid scale latent heat tendency using getgb
Index=78
VarName=avbl(index)
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=is(index)
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,train(1,jsta_2l,ll))
end do
! retrieve vertical diffusion moistening using getgb
Index=360
VarName='Vertical diffusion moistening'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,vdiffmois(1,jsta_2l,ll))
end do
! retrieve deep convection moistening using getgb
Index=361
VarName='deep convection moistening'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,dconvmois(1,jsta_2l,ll))
end do
! retrieve shallow convection moistening using getgb
Index=362
VarName='shallow convection moistening'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,sconvmois(1,jsta_2l,ll))
end do
! retrieve non-radiation tendency using getgb
Index=363
VarName='non-radiation tendency'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,nradtt(1,jsta_2l,ll))
end do
! retrieve Vertical diffusion of ozone using getgb
Index=364
VarName='Vertical diffusion of ozone'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,o3vdiff(1,jsta_2l,ll))
end do
! retrieve ozone production using getgb
Index=365
VarName='Ozone production'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,o3prod(1,jsta_2l,ll))
end do
! retrieve ozone tendency using getgb
Index=366
VarName='Ozone tendency'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,o3tndy(1,jsta_2l,ll))
end do
! retrieve mass weighted PV using getgb
Index=367
VarName='Mass weighted PV'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,mwpv(1,jsta_2l,ll))
end do
! retrieve OZONE TNDY using getgb
Index=368
VarName='?'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,unknown(1,jsta_2l,ll))
end do
! retrieve vertical diffusion zonal acceleration
Index=369
VarName='VDIFF Z ACCE'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,vdiffzacce(1,jsta_2l,ll))
end do
! retrieve gravity drag zonal acceleration
Index=370
VarName='G DRAG Z ACCE'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,zgdrag(1,jsta_2l,ll))
end do
! retrieve convective U momemtum mixing
Index=371
VarName='CNVCT U M MIX'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,cnvctummixing(1,jsta_2l,ll))
end do
! retrieve vertical diffusion meridional acceleration
Index=372
VarName='VDIFF M ACCE'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,vdiffmacce(1,jsta_2l,ll))
end do
! retrieve gravity drag meridional acceleration
Index=373
VarName='G DRAG M ACCE'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,mgdrag(1,jsta_2l,ll))
end do
! retrieve convective V momemtum mixing
Index=374
VarName='CNVCT V M MIX'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,cnvctvmmixing(1,jsta_2l,ll))
end do
! retrieve nonconvective cloud fraction
Index=375
VarName='N CNVCT CLD FRA'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,ncnvctcfrac(1,jsta_2l,ll))
end do
! retrieve convective upward mass flux
Index=391
VarName='CNVCT U M FLX'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,cnvctumflx(1,jsta_2l,ll))
end do
! retrieve convective downward mass flux
Index=392
VarName='CNVCT D M FLX'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,cnvctdmflx(1,jsta_2l,ll))
end do
! retrieve nonconvective detraintment flux
Index=393
VarName='CNVCT DET M FLX'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,cnvctdetmflx(1,jsta_2l,ll))
end do
! retrieve cnvct gravity drag zonal acceleration
Index=394
VarName='CNVCT G DRAG Z ACCE'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,cnvctzgdrag(1,jsta_2l,ll))
end do
! retrieve cnvct gravity drag meridional acceleration
Index=395
VarName='CNVCT G DRAG M ACCE'
jpds=-1.0
jgds=-1.0
jpds(5)=iq(index)
jpds(6)=109
do l=1,lm
jpds(7)=l
ll=lm-l+1 !flip 3d fields to count from top down
call getgbandscatter(me,iunitd3d,im,jm,im_jm,jsta,jsta_2l &
,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName &
,jpds,jgds,kpds,cnvctmgdrag(1,jsta_2l,ll))
end do
call baclose(iunitd3d,status)
if (me == 0) print*,'done reading D3D fields'
end if ! end of d3d file read
if (me == 0) print *,'after d3d files reading,mype=',me
! pos east
call collect_loc(gdlat,dummy)
if(me == 0) then
latstart=nint(dummy(1,1)*gdsdegr)
latlast=nint(dummy(im,jm)*gdsdegr)
print*,'laststart,latlast B bcast= ',latstart,latlast,'gdsdegr=',gdsdegr,&
'dummy(1,1)=',dummy(1,1),dummy(im,jm),'gdlat=',gdlat(1,1)
end if
call mpi_bcast(latstart,1,MPI_INTEGER,0,mpi_comm_comp,irtn)
call mpi_bcast(latlast,1,MPI_INTEGER,0,mpi_comm_comp,irtn)
write(6,*) 'laststart,latlast,me A calling bcast=',latstart,latlast,me
call collect_loc(gdlon,dummy)
if(me.eq.0)then
lonstart=nint(dummy(1,1)*gdsdegr)
lonlast=nint(dummy(im,jm)*gdsdegr)
end if
call mpi_bcast(lonstart,1,MPI_INTEGER,0,mpi_comm_comp,irtn)
call mpi_bcast(lonlast,1,MPI_INTEGER,0,mpi_comm_comp,irtn)
write(6,*)'lonstart,lonlast A calling bcast=',lonstart,lonlast
!
! ncdump -h
!!
!!
!!
write(6,*) 'filename in INITPOST=', filename,' is'
! status=nf_open(filename,NF_NOWRITE,ncid)
! write(6,*) 'returned ncid= ', ncid
! status=nf_get_att_real(ncid,varid,'DX',tmp)
! dxval=int(tmp)
! status=nf_get_att_real(ncid,varid,'DY',tmp)
! dyval=int(tmp)
! status=nf_get_att_real(ncid,varid,'CEN_LAT',tmp)
! cenlat=int(1000.*tmp)
! status=nf_get_att_real(ncid,varid,'CEN_LON',tmp)
! cenlon=int(1000.*tmp)
! status=nf_get_att_real(ncid,varid,'TRUELAT1',tmp)
! truelat1=int(1000.*tmp)
! status=nf_get_att_real(ncid,varid,'TRUELAT2',tmp)
! truelat2=int(1000.*tmp)
! status=nf_get_att_real(ncid,varid,'MAP_PROJ',tmp)
! maptype=int(tmp)
! status=nf_close(ncid)
! dxval=30000.
! dyval=30000.
!
! write(6,*) 'dxval= ', dxval
! write(6,*) 'dyval= ', dyval
! write(6,*) 'cenlat= ', cenlat
! write(6,*) 'cenlon= ', cenlon
! write(6,*) 'truelat1= ', truelat1
! write(6,*) 'truelat2= ', truelat2
! write(6,*) 'maptype is ', maptype
!
! close up shop
! call ext_int_ioclose ( DataHandle, Status )
! generate look up table for lifted parcel calculations
THL = 210.
PLQ = 70000.
! write(0,*)' bef TABLE PT=',PT,' THL=',THL
CALL TABLE(PTBL,TTBL,PT, &
RDQ,RDTH,RDP,RDTHE,PL,THL,QS0,SQS,STHE,THE0)
CALL TABLEQ(TTBLQ,RDPQ,RDTHEQ,PLQ,THL,STHEQ,THE0Q)
! write(0,*)'end ini_gfs_sigio'
!
!
IF(ME.EQ.0)THEN
WRITE(6,*)' SPL (POSTED PRESSURE LEVELS) BELOW: '
WRITE(6,51) (SPL(L),L=1,LSM)
50 FORMAT(14(F4.1,1X))
51 FORMAT(8(F8.1,1X))
ENDIF
!
! COMPUTE DERIVED TIME STEPPING CONSTANTS.
!
!MEB need to get DT
! DT = 120. !MEB need to get DT
! NPHS = 4 !MEB need to get physics DT
! TPREC=float(ifhr)
!MEB need to get DT
!how am i going to get this information?
! NPREC = INT(TPREC *TSPH+D50)
! NHEAT = INT(THEAT *TSPH+D50)
! NCLOD = INT(TCLOD *TSPH+D50)
! NRDSW = INT(TRDSW *TSPH+D50)
! NRDLW = INT(TRDLW *TSPH+D50)
! NSRFC = INT(TSRFC *TSPH+D50)
!how am i going to get this information?
!
! IF(ME.EQ.0)THEN
! WRITE(6,*)' '
! WRITE(6,*)'DERIVED TIME STEPPING CONSTANTS'
! WRITE(6,*)' NPREC,NHEAT,NSRFC : ',NPREC,NHEAT,NSRFC
! WRITE(6,*)' NCLOD,NRDSW,NRDLW : ',NCLOD,NRDSW,NRDLW
! ENDIF
!
! COMPUTE DERIVED MAP OUTPUT CONSTANTS.
DO L = 1,LSM
ALSL(L) = LOG(SPL(L))
END DO
!
!HC WRITE IGDS OUT FOR WEIGHTMAKER TO READ IN AS KGDSIN
if(me.eq.0)then
print*,'writing out igds'
igdout=110
! open(igdout,file='griddef.out',form='unformatted'
! + ,status='unknown')
if(maptype .eq. 1)THEN ! Lambert conformal
WRITE(igdout)3
WRITE(6,*)'igd(1)=',3
WRITE(igdout)im
WRITE(igdout)jm
WRITE(igdout)LATSTART
WRITE(igdout)LONSTART
WRITE(igdout)8
WRITE(igdout)CENLON
WRITE(igdout)DXVAL
WRITE(igdout)DYVAL
WRITE(igdout)0
WRITE(igdout)64
WRITE(igdout)TRUELAT2
WRITE(igdout)TRUELAT1
WRITE(igdout)255
ELSE IF(MAPTYPE .EQ. 2)THEN !Polar stereographic
WRITE(igdout)5
WRITE(igdout)im
WRITE(igdout)jm
WRITE(igdout)LATSTART
WRITE(igdout)LONSTART
WRITE(igdout)8
WRITE(igdout)CENLON
WRITE(igdout)DXVAL
WRITE(igdout)DYVAL
WRITE(igdout)0
WRITE(igdout)64
WRITE(igdout)TRUELAT2 !Assume projection at +-90
WRITE(igdout)TRUELAT1
WRITE(igdout)255
ELSE IF(MAPTYPE .EQ. 3)THEN !Mercator
WRITE(igdout)1
WRITE(igdout)im
WRITE(igdout)jm
WRITE(igdout)LATSTART
WRITE(igdout)LONSTART
WRITE(igdout)8
WRITE(igdout)latlast
WRITE(igdout)lonlast
WRITE(igdout)TRUELAT1
WRITE(igdout)0
WRITE(igdout)64
WRITE(igdout)DXVAL
WRITE(igdout)DYVAL
WRITE(igdout)255
ELSE IF(MAPTYPE.EQ.0 .OR. MAPTYPE.EQ.203)THEN !A STAGGERED E-GRID
WRITE(igdout)203
WRITE(igdout)im
WRITE(igdout)jm
WRITE(igdout)LATSTART
WRITE(igdout)LONSTART
WRITE(igdout)136
WRITE(igdout)CENLAT
WRITE(igdout)CENLON
WRITE(igdout)DXVAL
WRITE(igdout)DYVAL
WRITE(igdout)64
WRITE(igdout)0
WRITE(igdout)0
WRITE(igdout)0
END IF
end if
!
! close all files
call baclose(iunit,status)
deallocate(dummy)
RETURN
END