diff --git a/sorc/ncep_post.fd/ALLOCATE_ALL.f b/sorc/ncep_post.fd/ALLOCATE_ALL.f index 4643f84a5..afe9c7eed 100644 --- a/sorc/ncep_post.fd/ALLOCATE_ALL.f +++ b/sorc/ncep_post.fd/ALLOCATE_ALL.f @@ -18,6 +18,7 @@ !! - 21-04-06 Wen Meng - Initializing all allocated arrays !! - 21-04-16 Wen Meng - Initializing aextc55 and extc55 as 0. These !! two arrays are involved in GSL visibility computation. +!! - 22-03-22 Wen Meng - Initializing pwat. !! !! OUTPUT FILES: !! - STDOUT - RUN TIME STANDARD OUT. @@ -970,6 +971,7 @@ SUBROUTINE ALLOCATE_ALL() allocate(tedir(im,jsta_2l:jend_2u)) allocate(twa(im,jsta_2l:jend_2u)) allocate(fdnsst(im,jsta_2l:jend_2u)) + allocate(pwat(im,jsta_2l:jend_2u)) !Initialization !$omp parallel do private(i,j) do j=jsta_2l,jend_2u @@ -1020,6 +1022,7 @@ SUBROUTINE ALLOCATE_ALL() tedir(i,j)=spval twa(i,j)=spval fdnsst(i,j)=spval + pwat(i,j)=spval enddo enddo ! diff --git a/sorc/ncep_post.fd/CLDRAD.f b/sorc/ncep_post.fd/CLDRAD.f index df5799171..461fcac48 100644 --- a/sorc/ncep_post.fd/CLDRAD.f +++ b/sorc/ncep_post.fd/CLDRAD.f @@ -107,18 +107,19 @@ SUBROUTINE CLDRAD HBOT, HBOTD, HBOTS, HTOP, HTOPD, HTOPS, FIS, PBLH, & PBOT, PBOTL, PBOTM, PBOTH, CNVCFR, PTOP, PTOPL, & PTOPM, PTOPH, TTOPL, TTOPM, TTOPH, PBLCFR, CLDWORK, & - ASWIN, AUVBIN, AUVBINC, ASWIN, ASWOUT,ALWOUT, ASWTOA,& + ASWIN, AUVBIN, AUVBINC, ASWOUT,ALWOUT, ASWTOA, & RLWTOA, CZMEAN, CZEN, RSWIN, ALWIN, ALWTOA, RLWIN, & SIGT4, RSWOUT, RADOT, RSWINC, ASWINC, ASWOUTC, & ASWTOAC, ALWOUTC, ASWTOAC, AVISBEAMSWIN, & - AVISDIFFSWIN, ASWINTOA, ASWINC, ASWTOAC, AIRBEAMSWIN,& + AVISDIFFSWIN, ASWINTOA, ASWTOAC, AIRBEAMSWIN, & AIRDIFFSWIN, DUSMASS, DUSMASS25, DUCMASS, DUCMASS25, & ALWINC, ALWTOAC, SWDDNI, SWDDIF, SWDNBC, SWDDNIC, & SWDDIFC, SWUPBC, LWDNBC, LWUPBC, SWUPT, & TAOD5502D, AERSSA2D, AERASY2D, MEAN_FRP, LWP, IWP, & AVGCPRATE, & DUSTCB,SSCB,BCCB,OCCB,SULFCB,DUSTPM,SSPM,aod550, & - du_aod550,ss_aod550,su_aod550,oc_aod550,bc_aod550 + du_aod550,ss_aod550,su_aod550,oc_aod550,bc_aod550, & + PWAT use masks, only: LMH, HTM use params_mod, only: TFRZ, D00, H99999, QCLDMIN, SMALL, D608, H1, ROG, & GI, RD, QCONV, ABSCOEFI, ABSCOEF, STBOL, PQ0, A2, & @@ -254,6 +255,7 @@ SUBROUTINE CLDRAD data INDX_EXT / 610, 611, 612, 613, 614 / data INDX_SCA / 651, 652, 653, 654, 655 / logical, parameter :: debugprint = .false. + logical :: Model_Pwat ! ! !************************************************************************* @@ -417,12 +419,29 @@ SUBROUTINE CLDRAD IF (IGET(080) > 0) THEN ! dong GRID1 = spval + Model_Pwat = .false. + DO J=JSTA,JEND + DO I=1,IM + IF(ABS(PWAT(I,J)-SPVAL)>SMALL) THEN + Model_Pwat = .true. + exit + ENDIF + END DO + END DO + IF (Model_Pwat) THEN + DO J=JSTA,JEND + DO I=1,IM + GRID1(I,J) = PWAT(I,J) + END DO + END DO + ELSE CALL CALPW(GRID1(1,jsta),1) DO J=JSTA,JEND DO I=1,IM IF(FIS(I,J) >= SPVAL) GRID1(I,J)=spval END DO END DO + ENDIF CALL BOUND(GRID1,D00,H99999) if(grib == "grib2" )then cfld = cfld + 1 diff --git a/sorc/ncep_post.fd/CMakeLists.txt b/sorc/ncep_post.fd/CMakeLists.txt index 3215bad5f..470c02f6a 100644 --- a/sorc/ncep_post.fd/CMakeLists.txt +++ b/sorc/ncep_post.fd/CMakeLists.txt @@ -151,8 +151,6 @@ list(APPEND EXE_SRC GFSPOSTSIG.F INITPOST.F INITPOST_GFS_NEMS_MPIIO.f - INITPOST_GFS_NETCDF.f - INITPOST_GFS_NETCDF_PARA.f INITPOST_NEMS.f INITPOST_NETCDF.f WRFPOST.f diff --git a/sorc/ncep_post.fd/DEALLOCATE.f b/sorc/ncep_post.fd/DEALLOCATE.f index c6cf34ad9..f16cf8b6e 100644 --- a/sorc/ncep_post.fd/DEALLOCATE.f +++ b/sorc/ncep_post.fd/DEALLOCATE.f @@ -386,6 +386,7 @@ SUBROUTINE DE_ALLOCATE deallocate(tedir) deallocate(twa) deallocate(fdnsst) + deallocate(pwat) ! GSD deallocate(rainc_bucket) deallocate(rainc_bucket1) diff --git a/sorc/ncep_post.fd/INITPOST_GFS_NETCDF.f b/sorc/ncep_post.fd/INITPOST_GFS_NETCDF.f deleted file mode 100644 index b61732212..000000000 --- a/sorc/ncep_post.fd/INITPOST_GFS_NETCDF.f +++ /dev/null @@ -1,2761 +0,0 @@ -!> @file -! . . . -!> SUBPROGRAM: INITPOST_NETCDF INITIALIZE POST FOR RUN -!! PRGRMMR: Hui-Ya Chuang DATE: 2016-03-04 -!! -!! ABSTRACT: THIS ROUTINE INITIALIZES CONSTANTS AND -!! VARIABLES AT THE START OF GFS MODEL OR POST -!! PROCESSOR RUN. -!! -!! REVISION HISTORY -!! 2017-08-11 H Chuang start from INITPOST_GFS_NEMS_MPIIO.f -!! 2021-03-11 Bo Cui change local arrays to dimension (im,jsta:jend) -!! -!! USAGE: CALL INITPOST_NETCDF -!! 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 -!! - SUBROUTINE INITPOST_GFS_NETCDF(ncid3d) - - - use netcdf - use vrbls4d, only: dust, SALT, SUSO, SOOT, WASO, PP25, PP10 - use vrbls3d, only: t, q, uh, vh, pmid, pint, alpint, dpres, zint, zmid, o3, & - qqr, qqs, cwm, qqi, qqw, omga, rhomid, q2, cfr, rlwtt, rswtt, tcucn, & - tcucns, train, el_pbl, exch_h, vdifftt, vdiffmois, dconvmois, nradtt, & - o3vdiff, o3prod, o3tndy, mwpv, unknown, vdiffzacce, zgdrag,cnvctummixing, & - vdiffmacce, mgdrag, cnvctvmmixing, ncnvctcfrac, cnvctumflx, cnvctdmflx, & - cnvctzgdrag, sconvmois, cnvctmgdrag, cnvctdetmflx, duwt, duem, dusd, dudp, & - wh, qqg, ref_10cm - use vrbls2d, only: f, pd, fis, pblh, ustar, z0, ths, qs, twbs, qwbs, avgcprate, & - cprate, avgprec, prec, lspa, sno, si, cldefi, th10, q10, tshltr, pshltr, & - tshltr, albase, avgalbedo, avgtcdc, czen, czmean, mxsnal, radot, sigt4, & - cfrach, cfracl, cfracm, avgcfrach, qshltr, avgcfracl, avgcfracm, cnvcfr, & - islope, cmc, grnflx, vegfrc, acfrcv, ncfrcv, acfrst, ncfrst, ssroff, & - bgroff, rlwin, rlwtoa, cldwork, alwin, alwout, alwtoa, rswin, rswinc, & - rswout, aswin, auvbin, auvbinc, aswout, aswtoa, sfcshx, sfclhx, subshx, & - snopcx, sfcux, sfcvx, sfcuxi, sfcvxi, sfcuvx, gtaux, gtauy, potevp, u10, v10, smstav, & - smstot, ivgtyp, isltyp, sfcevp, sfcexc, acsnow, acsnom, sst, thz0, qz0, & - uz0, vz0, ptop, htop, pbot, hbot, ptopl, pbotl, ttopl, ptopm, pbotm, ttopm, & - ptoph, pboth, pblcfr, ttoph, runoff, maxtshltr, mintshltr, maxrhshltr, & - minrhshltr, dzice, smcwlt, suntime, fieldcapa, htopd, hbotd, htops, hbots, & - cuppt, dusmass, ducmass, dusmass25, ducmass25, aswintoa,rel_vort_maxhy1, & - maxqshltr, minqshltr, acond, sr, u10h, v10h,refd_max, w_up_max, w_dn_max, & - up_heli_max,up_heli_min,up_heli_max03,up_heli_min03,rel_vort_max01,u10max, v10max, & - avgedir,avgecan,avgetrans,avgesnow,avgprec_cont,avgcprate_cont,rel_vort_max, & - avisbeamswin,avisdiffswin,airbeamswin,airdiffswin,refdm10c_max,wspd10max, & - alwoutc,alwtoac,aswoutc,aswtoac,alwinc,aswinc,avgpotevp,snoavg, & - ti,aod550,du_aod550,ss_aod550,su_aod550,oc_aod550,bc_aod550 - use soil, only: sldpth, sh2o, smc, stc - use masks, only: lmv, lmh, htm, vtm, gdlat, gdlon, dx, dy, hbm2, sm, sice - use physcons_post, only: grav => con_g, fv => con_fvirt, rgas => con_rd, & - eps => con_eps, epsm1 => con_epsm1 - use params_mod, only: erad, dtr, tfrz, h1, d608, rd, p1000, capa,pi - use lookup_mod, only: thl, plq, ptbl, ttbl, rdq, rdth, rdp, rdthe, pl, qs0, sqs, sthe, & - ttblq, rdpq, rdtheq, stheq, the0q, the0 - use ctlblk_mod, only: me, mpi_comm_comp, icnt, idsp, jsta, jend, ihrst, idat, sdat, ifhr, & - ifmin, filename, tprec, tclod, trdlw, trdsw, tsrfc, tmaxmin, td3d, restrt, sdat, & - jend_m, imin, imp_physics, dt, spval, pdtop, pt, qmin, nbin_du, nphs, dtq2, ardlw,& - ardsw, asrfc, avrain, avcnvc, theat, gdsdegr, spl, lsm, alsl, im, jm, im_jm, lm, & - jsta_2l, jend_2u, nsoil, lp1, icu_physics, ivegsrc, novegtype, nbin_ss, nbin_bc, & - nbin_oc, nbin_su, gocart_on, pt_tbl, hyb_sigp, filenameFlux, fileNameAER,rdaod - use gridspec_mod, only: maptype, gridtype, latstart, latlast, lonstart, lonlast, cenlon, & - dxval, dyval, truelat2, truelat1, psmapf, cenlat,lonstartv, lonlastv, cenlonv, & - latstartv, latlastv, cenlatv,latstart_r,latlast_r,lonstart_r,lonlast_r - use upp_physics, only: fpvsnew -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none -! -! type(nemsio_gfile) :: nfile,ffile,rfile - integer,parameter :: nvar2d=48 -! character(nemsio_charkind) :: name2d(nvar2d) - integer :: nvar3d, numDims -! character(nemsio_charkind), allocatable :: name3din(:), name3dout(:) -! character(nemsio_charkind) :: varname,levtype -! -! 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 -! -! 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_. - - real, parameter :: gravi = 1.0/grav - character(len=20) :: VarName, VcoordName - integer :: Status, fldsize, fldst, recn, recn_vvel - 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,IOOMG,IOALL -! logical, parameter :: debugprint = .true., zerout = .false. - logical, parameter :: debugprint = .false., zerout = .false. - logical :: convert_rad_to_deg=.false. - CHARACTER*32 varcharval -! CHARACTER*40 CONTRL,FILALL,FILMST,FILTMP,FILTKE,FILUNV,FILCLD,FILRAD,FILSFC - CHARACTER*4 RESTHR - CHARACTER FNAME*255,ENVAR*50 - INTEGER IDATE(8),JDATE(8),JPDS(200),JGDS(200),KPDS(200),KGDS(200) -! LOGICAL*1 LB(IM,JM) -! -! INCLUDE COMMON BLOCKS. -! -! DECLARE VARIABLES. -! -! REAL fhour -! integer nfhour ! forecast hour from nems io file - integer fhzero !bucket - real dtp !physics time step - REAL RINC(5) - -! REAL FI(IM,JM,2) - REAL DUMMY(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, & - nframed2,iunitd3d,ierr,idum,iret,nrec,idrt - integer ncid3d,ncid2d,varid,nhcas - real TSTART,TLMH,TSPH,ES,FACT,soilayert,soilayerb,zhour,dum, & - tvll,pmll,tv, tx1, tx2 - - character*20,allocatable :: recname(:) - integer, allocatable :: reclev(:), kmsk(:,:) - real, allocatable :: glat1d(:), glon1d(:), qstl(:) - real, allocatable :: wrk1(:,:), wrk2(:,:) - real, allocatable :: p2d(:,:), t2d(:,:), q2d(:,:), & - qs2d(:,:), cw2d(:,:), cfr2d(:,:) - real(kind=4),allocatable :: vcoord4(:,:,:) - real, dimension(lm+1) :: ak5, bk5 - real*8, allocatable :: pm2d(:,:), pi2d(:,:) - real, allocatable :: tmp(:) - real :: buf(im,jsta_2l:jend_2u) - real :: buf3d(im,jsta_2l:jend_2u,lm) - -! real buf(im,jsta_2l:jend_2u),bufsoil(im,nsoil,jsta_2l:jend_2u) & -! ,buf3d(im,jsta_2l:jend_2u,lm),buf3d2(im,lp1,jsta_2l:jend_2u) - - real LAT - integer isa, jsa, latghf, jtem, idvc, idsl, nvcoord, ip1, nn, npass - - integer, parameter :: npass2=5, npass3=30 - real, parameter :: third=1.0/3.0 - INTEGER, DIMENSION(2) :: ij4min, ij4max - REAL :: omgmin, omgmax - real, allocatable :: d2d(:,:), u2d(:,:), v2d(:,:), omga2d(:,:) - REAL, ALLOCATABLE :: ps2d(:,:),psx2d(:,:),psy2d(:,:) - real, allocatable :: div3d(:,:,:) - real(kind=4),allocatable :: vcrd(:,:) - real :: dum_const - -!*********************************************************************** -! START INIT HERE. -! - WRITE(6,*)'INITPOST: ENTER INITPOST_GFS_NETCDF' - WRITE(6,*)'me=',me, & - 'jsta_2l=',jsta_2l,'jend_2u=', & - jend_2u,'im=',im -! - isa = im / 2 - jsa = (jsta+jend) / 2 - -!$omp parallel do private(i,j) - do j = jsta_2l, jend_2u - do i=1,im - buf(i,j) = spval - enddo - enddo - - Status=nf90_get_att(ncid3d,nf90_global,'ak',ak5) - if(Status/=0)then - print*,'ak not found; assigning missing value' - ak5=spval - else - if(me==0)print*,'ak5= ',ak5 - end if - Status=nf90_get_att(ncid3d,nf90_global,'idrt',idrt) - if(Status/=0)then - print*,'idrt not in netcdf file,reading grid' - Status=nf90_get_att(ncid3d,nf90_global,'grid',varcharval) - if(Status/=0)then - print*,'idrt and grid not in netcdf file, set default to latlon' - idrt=0 - MAPTYPE=0 - else - if(trim(varcharval)=='rotated_latlon')then - MAPTYPE=207 - idrt=207 - Status=nf90_get_att(ncid3d,nf90_global,'cen_lon',dum_const) - if(Status/=0)then - print*,'cen_lon not found; assigning missing value' - cenlon=spval - else - if(dum_const<0.)then - cenlon=nint((dum_const+360.)*gdsdegr) - else - cenlon=dum_const*gdsdegr - end if - end if - Status=nf90_get_att(ncid3d,nf90_global,'cen_lat',dum_const) - if(Status/=0)then - print*,'cen_lat not found; assigning missing value' - cenlat=spval - else - cenlat=dum_const*gdsdegr - end if - - Status=nf90_get_att(ncid3d,nf90_global,'lon1',dum_const) - if(Status/=0)then - print*,'lonstart_r not found; assigning missing value' - lonstart_r=spval - else - if(dum_const<0.)then - lonstart_r=nint((dum_const+360.)*gdsdegr) - else - lonstart_r=dum_const*gdsdegr - end if - end if - Status=nf90_get_att(ncid3d,nf90_global,'lat1',dum_const) - if(Status/=0)then - print*,'latstart_r not found; assigning missing value' - latstart_r=spval - else - latstart_r=dum_const*gdsdegr - end if - - Status=nf90_get_att(ncid3d,nf90_global,'lon2',dum_const) - if(Status/=0)then - print*,'lonlast_r not found; assigning missing value' - lonlast_r=spval - else - if(dum_const<0.)then - lonlast_r=nint((dum_const+360.)*gdsdegr) - else - lonlast_r=dum_const*gdsdegr - end if - end if - Status=nf90_get_att(ncid3d,nf90_global,'lat2',dum_const) - if(Status/=0)then - print*,'latlast_r not found; assigning missing value' - latlast_r=spval - else - latlast_r=dum_const*gdsdegr - end if - - Status=nf90_get_att(ncid3d,nf90_global,'dlon',dum_const) - if(Status/=0)then - print*,'dlmd not found; assigning missing value' - dxval=spval - else - dxval=dum_const*gdsdegr - end if - Status=nf90_get_att(ncid3d,nf90_global,'dlat',dum_const) - if(Status/=0)then - print*,'dphd not found; assigning missing value' - dyval=spval - else - dyval=dum_const*gdsdegr - end if - - print*,'lonstart,latstart,cenlon,cenlat,dyval,dxval', & - lonstart,latstart,cenlon,cenlat,dyval,dxval - -! Jili Dong add support for regular lat lon (2019/03/22) start - else if(trim(varcharval)=='latlon')then - MAPTYPE=0 - idrt=0 - - Status=nf90_get_att(ncid3d,nf90_global,'lon1',dum_const) - if(Status/=0)then - print*,'lonstart not found; assigning missing value' - lonstart=spval - else - if(dum_const<0.)then - lonstart=nint((dum_const+360.)*gdsdegr) - else - lonstart=dum_const*gdsdegr - end if - end if - Status=nf90_get_att(ncid3d,nf90_global,'lat1',dum_const) - if(Status/=0)then - print*,'latstart not found; assigning missing value' - latstart=spval - else - latstart=dum_const*gdsdegr - end if - - Status=nf90_get_att(ncid3d,nf90_global,'lon2',dum_const) - if(Status/=0)then - print*,'lonlast not found; assigning missing value' - lonlast=spval - else - if(dum_const<0.)then - lonlast=nint((dum_const+360.)*gdsdegr) - else - lonlast=dum_const*gdsdegr - end if - end if - Status=nf90_get_att(ncid3d,nf90_global,'lat2',dum_const) - if(Status/=0)then - print*,'latlast not found; assigning missing value' - latlast=spval - else - latlast=dum_const*gdsdegr - end if - - Status=nf90_get_att(ncid3d,nf90_global,'dlon',dum_const) - if(Status/=0)then - print*,'dlmd not found; assigning missing value' - dxval=spval - else - dxval=dum_const*gdsdegr - end if - Status=nf90_get_att(ncid3d,nf90_global,'dlat',dum_const) - if(Status/=0)then - print*,'dphd not found; assigning missing value' - dyval=spval - else - dyval=dum_const*gdsdegr - end if - - print*,'lonstart,latstart,dyval,dxval', & - lonstart,lonlast,latstart,latlast,dyval,dxval - -! Jili Dong add support for regular lat lon (2019/03/22) end - - else if(trim(varcharval)=='gaussian')then - MAPTYPE=4 - idrt=4 - else ! setting default maptype - MAPTYPE=0 - idrt=0 - end if - end if !end reading grid - end if !end reading idrt - if(me==0)print*,'idrt MAPTYPE= ',idrt,MAPTYPE -! STEP 1. READ MODEL OUTPUT FILE -! -! -!*** -! -! LMH and 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 - -! 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 - - Status=nf90_get_att(ncid3d,nf90_global,'nhcas',nhcas) - if(Status/=0)then - print*,'nhcas not in netcdf file, set default to nonhydro' - nhcas=0 - end if - if(me==0)print*,'nhcas= ',nhcas - if (nhcas == 0 ) then !non-hydrostatic case - nrec=15 - allocate (recname(nrec)) - recname=[character(len=20) :: 'ugrd','vgrd','spfh','tmp','o3mr', & - 'presnh','dzdt', 'clwmr','dpres', & - 'delz','icmr','rwmr', & - 'snmr','grle','cld_amt'] - else - nrec=8 - allocate (recname(nrec)) - recname=[character(len=20) :: 'ugrd','vgrd','tmp','spfh','o3mr', & - 'hypres', 'clwmr','dpres'] - endif - -! write(0,*)'nrec=',nrec - !allocate(recname(nrec),reclevtyp(nrec),reclev(nrec)) - allocate(glat1d(jm),glon1d(im)) - allocate(vcoord4(lm+1,3,2)) - -! hardwire idate for now -! idate=(/2017,08,07,00,0,0,0,0/) -! get cycle start time - Status=nf90_inq_varid(ncid3d,'time',varid) - if(Status/=0)then - print*,'time not in netcdf file, stopping' - STOP 1 - else - Status=nf90_get_att(ncid3d,varid,'units',varcharval) - if(Status/=0)then - print*,'time unit not available' - else - print*,'time unit read from netcdf file= ',varcharval -! assume use hours as unit -! idate_loc=index(varcharval,'since')+6 - read(varcharval,101)idate(1),idate(2),idate(3),idate(4),idate(5) - end if -! Status=nf90_inquire_dimension(ncid3d,varid,len=ntimes) -! allocate(fhours(ntimes)) -! status = nf90_inq_varid(ncid3d,varid,fhours) -! Status=nf90_get_var(ncid3d,varid,nfhour,start=(/1/), & -! count=(/1/)) -! if(Status/=0)then -! print*,'forecast hour not in netcdf file, stopping' -! STOP 1 -! end if - end if - 101 format(T13,i4,1x,i2,1x,i2,1x,i2,1x,i2) - print*,'idate= ',idate(1:5) -! get longitude - Status=nf90_inq_varid(ncid3d,'grid_xt',varid) - Status=nf90_inquire_variable(ncid3d,varid,ndims = numDims) - if(debugprint)print*,'number of dim for gdlon ',numDims - if(numDims==1)then - Status=nf90_get_var(ncid3d,varid,glon1d) - do j=jsta,jend - do i=1,im - gdlon(i,j) = real(glon1d(i),kind=4) - end do - end do - lonstart = nint(glon1d(1)*gdsdegr) - lonlast = nint(glon1d(im)*gdsdegr) - dxval = nint(abs(glon1d(1)-glon1d(2))*gdsdegr) - else if(numDims==2)then - Status=nf90_get_var(ncid3d,varid,dummy) - if(maxval(abs(dummy))<2.0*pi)convert_rad_to_deg=.true. - if(convert_rad_to_deg)then - do j=jsta,jend - do i=1,im - gdlon(i,j) = real(dummy(i,j),kind=4)*180./pi - end do - end do - else - do j=jsta,jend - do i=1,im - gdlon(i,j) = real(dummy(i,j),kind=4) - end do - end do - end if - if(convert_rad_to_deg)then - lonstart = nint(dummy(1,1)*gdsdegr)*180./pi - lonlast = nint(dummy(im,jm)*gdsdegr)*180./pi - dxval = nint(abs(dummy(1,1)-dummy(2,1))*gdsdegr)*180./pi - else - lonstart = nint(dummy(1,1)*gdsdegr) - lonlast = nint(dummy(im,jm)*gdsdegr) - dxval = nint(abs(dummy(1,1)-dummy(2,1))*gdsdegr) - end if - -! Jili Dong add support for regular lat lon (2019/03/22) start - if (MAPTYPE == 0) then - if(lonstart<0.)then - lonstart=lonstart+360.*gdsdegr - end if - if(lonlast<0.)then - lonlast=lonlast+360.*gdsdegr - end if - end if -! Jili Dong add support for regular lat lon (2019/03/22) end - - end if - print*,'lonstart,lonlast,dxval ',lonstart,lonlast,dxval -! get latitude - Status=nf90_inq_varid(ncid3d,'grid_yt',varid) - Status=nf90_inquire_variable(ncid3d,varid,ndims = numDims) - if(debugprint)print*,'number of dim for gdlat ',numDims - if(numDims==1)then - Status=nf90_get_var(ncid3d,varid,glat1d) - do j=jsta,jend - do i=1,im - gdlat(i,j) = real(glat1d(j),kind=4) - end do - end do - latstart = nint(glat1d(1)*gdsdegr) - latlast = nint(glat1d(jm)*gdsdegr) - dyval = nint(abs(glat1d(1)-glat1d(2))*gdsdegr) - else if(numDims==2)then - Status=nf90_get_var(ncid3d,varid,dummy) - if(maxval(abs(dummy))1.0e5)print*,'bad dpres ',i,j,dpres(i,j,l) -!make sure delz is positive -! if(dpres(i,j,l)/=spval .and. t(i,j,l)/=spval .and. & -! q(i,j,l)/=spval .and. buf3d(i,j,l)/=spval)then -! pmid(i,j,l)=rgas*dpres(i,j,l)* & -! t(i,j,l)*(q(i,j,l)*fv+1.0)/grav/abs(buf3d(i,j,l)) -! else -! pmid(i,j,l)=spval -! end if -! dong add missing value - if (wh(i,j,l) < spval) then - omga(i,j,l)=(-1.)*wh(i,j,l)*dpres(i,j,l)/abs(buf3d(i,j,l)) - else - omga(i,j,l) = spval - end if -! if(t(i,j,l)>1000.)print*,'bad T ',t(i,j,l) - enddo - enddo - enddo - call read_netcdf_3d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,recname(11) & - ,lm,qqi(1,jsta_2l,1)) - call read_netcdf_3d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,recname(12) & - ,lm,qqr(1,jsta_2l,1)) - call read_netcdf_3d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,recname(13) & - ,lm,qqs(1,jsta_2l,1)) - call read_netcdf_3d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,recname(14) & - ,lm,qqg(1,jsta_2l,1)) - call read_netcdf_3d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,recname(15) & - ,lm,cfr(1,jsta_2l,1)) -! calculate CWM from FV3 output - do l=1,lm - do j=jsta,jend - do i=1,im - cwm(i,j,l)=qqg(i,j,l)+qqs(i,j,l)+qqr(i,j,l)+qqi(i,j,l)+qqw(i,j,l) - enddo - enddo - if(debugprint)print*,'sample l,t,q,u,v,w,= ',isa,jsa,l & - ,t(isa,jsa,l),q(isa,jsa,l),uh(isa,jsa,l),vh(isa,jsa,l) & - ,wh(isa,jsa,l) - if(debugprint)print*,'sample l cwm for FV3',l, & - cwm(isa,jsa,l) - end do -! max hourly updraft velocity -! VarName='upvvelmax' -! call read_netcdf_2d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,w_up_max) -! if(debugprint)print*,'sample ',VarName,' = ',w_up_max(isa,jsa) - -! max hourly downdraft velocity -! VarName='dnvvelmax' -! call read_netcdf_2d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,w_dn_max) -! if(debugprint)print*,'sample ',VarName,' = ',w_dn_max(isa,jsa) -! max hourly updraft helicity -! VarName='uhmax25' -! call read_netcdf_2d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,up_heli_max) -! if(debugprint)print*,'sample ',VarName,' = ',up_heli_max(isa,jsa) -! min hourly updraft helicity -! VarName='uhmin25' -! call read_netcdf_2d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,up_heli_min) -! if(debugprint)print*,'sample ',VarName,' = ',up_heli_min(isa,jsa) -! max hourly 0-3km updraft helicity -! VarName='uhmax03' -! call read_netcdf_2d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,up_heli_max03) -! if(debugprint)print*,'sample ',VarName,' = ',up_heli_max03(isa,jsa) -! min hourly 0-3km updraft helicity -! VarName='uhmin03' -! call read_netcdf_2d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,up_heli_min03) -! if(debugprint)print*,'sample ',VarName,' = ',up_heli_min03(isa,jsa) - -! max 0-1km relative vorticity max -! VarName='maxvort01' -! call read_netcdf_2d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,rel_vort_max01) -! if(debugprint)print*,'sample ',VarName,' = ',rel_vort_max01(isa,jsa) -! max 0-2km relative vorticity max -! VarName='maxvort02' -! call read_netcdf_2d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,rel_vort_max) -! if(debugprint)print*,'sample ',VarName,' =',rel_vort_max(isa,jsa) -! max hybrid lev 1 relative vorticity max -! VarName='maxvorthy1' -! call read_netcdf_2d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,rel_vort_maxhy1) -! if(debugprint)print*,'sample ',VarName,' =',rel_vort_maxhy1(isa,jsa) -! surface pressure - VarName='pressfc' - call read_netcdf_2d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName & - ,pint(1,jsta_2l,lp1)) - do j=jsta,jend - do i=1,im -! if(pint(i,j,lp1)>1.0E6 .or. pint(1,jsta_2l,lp1)<50000.) & -! print*,'bad psfc ',i,j,pint(i,j,lp1) - end do - end do - if(debugprint)print*,'sample ',VarName,' =',pint(isa,jsa,lp1) - - pt = ak5(1) - - do j=jsta,jend - do i=1,im - pint(i,j,1)= pt - end do - end do - - do l=2,lp1 - do j=jsta,jend - do i=1,im - pint(i,j,l) = pint(i,j,l-1) + dpres(i,j,l-1) - enddo - enddo -! if (me == 0) print*,'sample model pint,pmid' ,ii,jj,l & -! ,pint(ii,jj,l),pmid(ii,jj,l) - end do - -!compute pmid from averaged two layer pint - do l=lm,1,-1 - do j=jsta,jend - do i=1,im - pmid(i,j,l) = 0.5*(pint(i,j,l)+pint(i,j,l+1)) - enddo - enddo - enddo - -! do l=lm,1,-1 -! do j=jsta,jend -! do i=1,im -! if(pint(i,j,l+1)/=spval .and. dpres(i,j,l)/=spval)then -! pint(i,j,l)=pint(i,j,l+1)-dpres(i,j,l) -! else -! pint(i,j,l)=spval -! end if -! end do -! end do -! print*,'sample pint= ',isa,jsa,l,pint(isa,jsa,l) -! end do - -! surface height from FV3 -! dong set missing value for zint -! zint=spval - VarName='hgtsfc' - call read_netcdf_2d_scatter(me,ncid3d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName & - ,zint(1,jsta_2l,lp1)) - if(debugprint)print*,'sample ',VarName,' =',zint(isa,jsa,lp1) - do j=jsta,jend - do i=1,im - if (zint(i,j,lp1) /= spval) then - fis(i,j) = zint(i,j,lp1) * grav - else - fis(i,j) = spval - endif - enddo - enddo - - do l=lm,1,-1 - do j=jsta,jend - do i=1,im - if(zint(i,j,l+1)/=spval .and. buf3d(i,j,l)/=spval)then -!make sure delz is positive - zint(i,j,l)=zint(i,j,l+1)+abs(buf3d(i,j,l)) -! if(zint(i,j,l)>1.0E6)print*,'bad H ',i,j,l,zint(i,j,l) - else - zint(i,j,l)=spval - end if - end do - end do - print*,'sample zint= ',isa,jsa,l,zint(isa,jsa,l) - end do - - do l=lp1,1,-1 - do j=jsta,jend - do i=1,im - alpint(i,j,l)=log(pint(i,j,l)) - end do - end do - end do - - do l=lm,1,-1 - do j=jsta,jend - do i=1,im - if(zint(i,j,l+1)/=spval .and. zint(i,j,l)/=spval & - .and. pmid(i,j,l)/=spval)then - zmid(i,j,l)=zint(i,j,l+1)+(zint(i,j,l)-zint(i,j,l+1))* & - (log(pmid(i,j,l))-alpint(i,j,l+1))/ & - (alpint(i,j,l)-alpint(i,j,l+1)) - if(zmid(i,j,l)>1.0E6)print*,'bad Hmid ',i,j,l,zmid(i,j,l) - else - zmid(i,j,l)=spval - endif - end do - end do - end do - - - pt = ak5(1) - -! else -! do l=2,lm -!!$omp parallel do private(i,j) -! do j=jsta,jend -! do i=1,im -! pint(i,j,l) = pint(i,j,l-1) + dpres(i,j,l-1) -! enddo -! enddo -! if (me == 0) print*,'sample pint,pmid' ,ii,jj,l,pint(ii,jj,l),pmid(ii,jj,l) -! end do -! endif -! - - deallocate (vcoord4) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - -! - -! done with 3d file, close it for now - Status=nf90_close(ncid3d) - deallocate(recname) - -! open flux file - Status = nf90_open(trim(fileNameFlux),NF90_NOWRITE, ncid2d) - - if ( Status /= 0 ) then - print*,'error opening ',fileNameFlux, ' Status = ', Status - print*,'skip reading of flux file' - endif - -! IVEGSRC=1 for IGBP, 0 for USGS, 2 for UMD - VarName='IVEGSRC' - Status=nf90_get_att(ncid2d,nf90_global,'IVEGSRC',IVEGSRC) - if (Status /= 0) then - print*,VarName,' not found-Assigned 1 for IGBP as default' - IVEGSRC=1 - end if - if (me == 0) print*,'IVEGSRC= ',IVEGSRC - -! set novegtype based on vegetation classification - if(ivegsrc==2)then - novegtype=13 - else if(ivegsrc==1)then - novegtype=20 - else if(ivegsrc==0)then - novegtype=24 - end if - if (me == 0) print*,'novegtype= ',novegtype - - Status=nf90_get_att(ncid2d,nf90_global,'imp_physics',imp_physics) - if (Status /= 0) then - print*,VarName,' not found-Assigned 11 GFDL as default' - imp_physics=11 - end if - if (me == 0) print*,'MP_PHYSICS= ',imp_physics -! - Status=nf90_get_att(ncid2d,nf90_global,'fhzero',fhzero) - if (Status /= 0) then - print*,VarName,' not found-Assigned 3 hours as default' - fhzero=3 - end if - if (me == 0) print*,'fhzero= ',fhzero -! - Status=nf90_get_att(ncid2d,nf90_global,'dtp',dtp) - if (Status /= 0) then - print*,VarName,' not found-Assigned 90s as default' - dtp=90 - end if - if (me == 0) print*,'dtp= ',dtp -! Initializes constants for Ferrier microphysics - if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95) then - CALL MICROINIT(imp_physics) - end if - -! Chuang: zhour is when GFS empties bucket last so using this -! to compute buket will result in changing bucket with forecast time. -! set default bucket for now - -! call nemsio_getheadvar(ffile,'zhour',zhour,iret=iret) -! if(iret == 0) then -! tprec = 1.0*ifhr-zhour -! tclod = tprec -! trdlw = tprec -! trdsw = tprec -! tsrfc = tprec -! tmaxmin = tprec -! td3d = tprec -! print*,'tprec from flux file header= ',tprec -! else -! print*,'Error reading accumulation bucket from flux file', & -! 'header - will try to read from env variable FHZER' -! CALL GETENV('FHZER',ENVAR) -! read(ENVAR, '(I2)')idum -! tprec = idum*1.0 -! tclod = tprec -! trdlw = tprec -! trdsw = tprec -! tsrfc = tprec -! tmaxmin = tprec -! td3d = tprec -! print*,'TPREC from FHZER= ',tprec -! end if - - - tprec = float(fhzero) - if(ifhr>240)tprec=12. - tclod = tprec - trdlw = tprec - trdsw = tprec - tsrfc = tprec - tmaxmin = tprec - td3d = tprec - print*,'tprec = ',tprec - - -! start reading 2d netcdf file -! surface pressure -! VarName='pressfc' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName & -! ,pint(1,jsta_2l,lp1)) -! if(debugprint)print*,'sample ',VarName,' =',pint(isa,jsa,lp1) -! do l=lm,1,-1 -! do j=jsta,jend -! do i=1,im -! pint(i,j,l)=pint(i,j,l+1)-dpres(i,j,l) -! if(pint(i,j,l)>1.0E6)print*,'bad P ',i,j,l,pint(i,j,l) & -! ,pint(i,j,l+1),dpres(i,j,l) -! end do -! end do -! print*,'sample pint= ',isa,jsa,l,pint(isa,jsa,l) -! end do -! surface height from FV3 already multiplied by G -! VarName='orog' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,fis) -! if(debugprint)print*,'sample ',VarName,' =',fis(isa,jsa) -! do j=jsta,jend -! do i=1,im -! if (fis(i,j) /= spval) then -! zint(i,j,lp1) = fis(i,j) -! fis(i,j) = fis(i,j) * grav -! else -! zint(i,j,lp1) = spval -! fis(i,j) = spval -! endif -! enddo -! enddo - -! do l=lm,1,-1 -! do j=jsta,jend -! do i=1,im -! if(zint(i,j,l+1)/=spval .and. buf3d(i,j,l)/=spval)then -! zint(i,j,l)=zint(i,j,l+1)+buf3d(i,j,l) -! if(zint(i,j,l)>1.0E6)print*,'bad H ',i,j,l,zint(i,j,l) -! else -! zint(i,j,l)=spval -! end if -! end do -! end do -! print*,'sample zint= ',isa,jsa,l,zint(isa,jsa,l) -! end do - -! Per communication with Fanglin, P from model in not monotonic -! so compute P using ak and bk for now Sep. 2017 -! do l=lm,1,-1 -!!!$omp parallel do private(i,j) -! do j=jsta,jend -! do i=1,im -! pint(i,j,l) = ak5(lm+2-l) + bk5(lm+2-l)*pint(i,j,lp1) -! pmid(i,j,l) = 0.5*(pint(i,j,l)+pint(i,j,l+1)) ! for now - -! enddo -! enddo -! print*,'sample pint,pmid' & -! ,l,pint(isa,jsa,l),pmid(isa,jsa,l) -! enddo - -! allocate(wrk1(im,jsta:jend),wrk2(im,jsta:jend)) -! do j=jsta,jend -! do i=1,im -! pd(i,j) = spval ! GFS does not output PD -! pint(i,j,1) = PT -! 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)*fv+1.0) -! FI(I,J,1) = FIS(I,J) & -! + wrk2(i,j)*rgas*(ALPINT(I,J,Lp1)-wrk1(i,j)) -! ZMID(I,J,LM) = FI(I,J,1) * gravi -! end do -! end do - -! SECOND, INTEGRATE HEIGHT HYDROSTATICLY, GFS integrate height on -! mid-layer - -! DO L=LM,2,-1 ! omit computing model top height -! ll = l - 1 -! do j = jsta, jend -! do i = 1, im -! alpint(i,j,l) = log(pint(i,j,l)) -! tvll = T(I,J,LL)*(Q(I,J,LL)*fv+1.0) -! pmll = log(PMID(I,J,LL)) - -! FI(I,J,2) = FI(I,J,1) + (0.5*rgas)*(wrk2(i,j)+tvll) & -! * (wrk1(i,j)-pmll) -! ZMID(I,J,LL) = FI(I,J,2) * gravi -! -! 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 -! FI(I,J,1) = FI(I,J,2) -! wrk1(i,J) = pmll -! wrk2(i,j) = tvll -! ENDDO -! ENDDO - -! print*,'L ZINT= ',l,zint(isa,jsa,l),ZMID(isa,jsa,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) - -! do l=lp1,2,-1 -! do j=jsta,jend -! do i=1,im -! alpint(i,j,l)=log(pint(i,j,l)) -! end do -! end do -! end do - -! do l=lm,2,-1 -! do j=jsta,jend -! do i=1,im -! zmid(i,j,l)=zint(i,j,l+1)+(zint(i,j,l)-zint(i,j,l+1))* & -! (log(pmid(i,j,l))-alpint(i,j,l+1))/ & -! (alpint(i,j,l)-alpint(i,j,l+1)) -! if(zmid(i,j,l)>1.0E6)print*,'bad Hmid ',i,j,l,zmid(i,j,l) -! end do -! end do -! end do - -! VarName='refl_10cm' -! do l=1,lm -! call read_netcdf_3d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName & -! ,lm,REF_10CM(1,jsta_2l,1)) -! if(debugprint)print*,'sample ',VarName,'isa,jsa,l =' & -! ,REF_10CM(isa,jsa,l),isa,jsa,l -! enddo -!Set REF_10CM as missning since gfs doesn't ouput it - do l=1,lm - do j=jsta,jend - do i=1,im - REF_10CM(i,j,l)=spval - enddo - enddo - enddo - - VarName='land' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sm) - if(debugprint)print*,'sample ',VarName,' =',sm(im/2,(jsta+jend)/2) - -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= spval) sm(i,j) = 1.0 - sm(i,j) - enddo - enddo - -! sea ice mask - - VarName = 'icec' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sice) - if(debugprint)print*,'sample ',VarName,' = ',sice(isa,jsa) - -! where(sice /=spval .and. sice >=1.0)sm=0.0 !sea ice has sea -! mask=0 -! GFS flux files have land points with non-zero sea ice, per Iredell, -! these -! points have sea ice changed to zero, i.e., trust land mask more than -! sea ice -! where(sm/=spval .and. sm==0.0)sice=0.0 !specify sea ice=0 at land - -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= spval .and. sm(i,j) == 0.0) sice(i,j) = 0.0 - enddo - enddo - - -! PBL height using nemsio - VarName = 'hpbl' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,pblh) - if(debugprint)print*,'sample ',VarName,' = ',pblh(isa,jsa) - -! frictional velocity using nemsio - VarName='fricv' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,ustar) -! if(debugprint)print*,'sample ',VarName,' = ',ustar(isa,jsa) - -! roughness length using getgb - VarName='sfcr' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,z0) -! if(debugprint)print*,'sample ',VarName,' = ',z0(isa,jsa) - -! sfc exchange coeff - VarName='sfexc' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,SFCEXC) - -! aerodynamic conductance - VarName='acond' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,acond) - if(debugprint)print*,'sample ',VarName,' = ',acond(isa,jsa) -! mid day avg albedo - VarName='albdo_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgalbedo) - if(debugprint)print*,'sample ',VarName,' = ',avgalbedo(isa,jsa) - do j=jsta,jend - do i=1,im - if (avgalbedo(i,j) /= spval) avgalbedo(i,j) = avgalbedo(i,j) * 0.01 - enddo - enddo - -! surface potential T using getgb - VarName='tmpsfc' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,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 - QS(i,j) = SPVAL ! GFS does not have surface specific humidity - twbs(i,j) = SPVAL ! GFS does not have inst sensible heat flux - qwbs(i,j) = SPVAL ! GFS does not have inst latent heat flux -!assign sst - if (sm(i,j) /= 0.0) then - sst(i,j) = ths(i,j) * (pint(i,j,lp1)/p1000)**capa - else - sst(i,j) = spval - endif - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',ths(isa,jsa) - - -! GFS does not have time step and physics time step, make up ones since they -! are not really used anyway -! NPHS=1. -! DT=90. -! DTQ2 = DT * NPHS !MEB need to get physics DT - DTQ2 = DTP !MEB need to get physics DT - NPHS=1 - DT = DTQ2/NPHS !MEB need to get DT - TSPH = 3600./DT - -! convective precip in m per physics time step using getgb -! read 3 hour bucket - VarName='cpratb_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgcprate) -! where(avgcprate /= spval)avgcprate=avgcprate*dtq2/1000. ! convert to m -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgcprate(i,j) /= spval) avgcprate(i,j) = avgcprate(i,j) * (dtq2*0.001) - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',avgcprate(isa,jsa) - -! print*,'maxval CPRATE: ', maxval(CPRATE) - -! read continuous bucket - VarName='cprat_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgcprate_cont) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgcprate_cont(i,j) /= spval) avgcprate_cont(i,j) = & - avgcprate_cont(i,j) * (dtq2*0.001) - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',avgcprate_cont(isa,jsa) - -! print*,'maxval CPRATE: ', maxval(CPRATE) - -! precip rate in m per physics time step using getgb - VarName='prateb_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgprec) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if(avgprec(i,j) /= spval)avgprec(i,j)=avgprec(i,j)*(dtq2*0.001) - enddo - enddo - - if(debugprint)print*,'sample ',VarName,' = ',avgprec(isa,jsa) - -! prec = avgprec !set avg cprate to inst one to derive other fields - - VarName='prate_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgprec_cont) -! where(avgprec /= spval)avgprec=avgprec*dtq2/1000. ! convert to m -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgprec_cont(i,j) /=spval)avgprec_cont(i,j)=avgprec_cont(i,j) & - * (dtq2*0.001) - enddo - enddo - - if(debugprint)print*,'sample ',VarName,' = ',avgprec_cont(isa,jsa) -! precip rate in m per physics time step - VarName='tprcp' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,prec) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (prec(i,j) /= spval) prec(i,j)=prec(i,j)* (dtq2*0.001) & - * 1000. / dtp - enddo - enddo - -! convective precip rate in m per physics time step - VarName='cnvprcp' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,cprate) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (cprate(i,j) /= spval) then - cprate(i,j) = max(0.,cprate(i,j)) * (dtq2*0.001) & - * 1000. / dtp - else - cprate(i,j) = 0. - endif - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',cprate(isa,jsa) - -! GFS does not have accumulated total, gridscale, and convective precip, will use inst precip to derive in SURFCE.f - -! max hourly 1-km agl reflectivity -! VarName='refdmax' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,refd_max) -! if(debugprint)print*,'sample ',VarName,' = ',refd_max(isa,jsa) -! max hourly -10C reflectivity -! VarName='refdmax263k' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,refdm10c_max) -! if(debugprint)print*,'sample ',VarName,' = ',refdm10c_max(isa,jsa) - -! max hourly u comp of 10m agl wind -! VarName='u10max' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,u10max) -! if(debugprint)print*,'sample ',VarName,' = ',u10max(isa,jsa) -! max hourly v comp of 10m agl wind -! VarName='v10max' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,v10max) -! if(debugprint)print*,'sample ',VarName,' = ',v10max(isa,jsa) -! max hourly 10m agl wind speed -! VarName='spd10max' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,wspd10max) -! if(debugprint)print*,'sample ',VarName,' = ',wspd10max(isa,jsa) - - -! 2m T using nemsio - VarName='tmp2m' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,tshltr) - if(debugprint)print*,'sample ',VarName,' = ',tshltr(isa,jsa) - -! inst snow water eqivalent using nemsio - VarName='weasd' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sno) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j)==0.) sno(i,j) = spval - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',sno(isa,jsa) - -! ave snow cover - VarName='snowc_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,snoavg) -! snow cover is multipled by 100 in SURFCE before writing it out - do j=jsta,jend - do i=1,im - if (sm(i,j)==1.0 .and. sice(i,j)==0.) snoavg(i,j)=spval - if(snoavg(i,j)/=spval)snoavg(i,j)=snoavg(i,j)/100. - end do - end do - -! snow depth in mm using nemsio - VarName='snod' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,si) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j)==1.0 .and. sice(i,j)==0.) si(i,j)=spval - if (si(i,j) /= spval) si(i,j) = si(i,j) * 1000.0 - CLDEFI(i,j) = SPVAL ! GFS does not have convective cloud efficiency - lspa(i,j) = spval ! GFS does not have similated precip - TH10(i,j) = SPVAL ! GFS does not have 10 m theta - TH10(i,j) = SPVAL ! GFS does not have 10 m theta - Q10(i,j) = SPVAL ! GFS does not have 10 m humidity - ALBASE(i,j) = SPVAL ! GFS does not have snow free albedo - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',si(isa,jsa) - -! 2m T using nemsio - VarName='tmp2m' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,tshltr) - if(debugprint)print*,'sample ',VarName,' = ',tshltr(isa,jsa) - -! GFS does not have 2m pres, estimate it, also convert t to theta - Do j=jsta,jend - Do i=1,im - 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 -! if (j == jm/2 .and. mod(i,50) == 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 nemsio - VarName='spfh2m' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,qshltr) - if(debugprint)print*,'sample ',VarName,' = ',qshltr(isa,jsa) - -! mid day avg albedo in fraction using nemsio -! VarName='albdosfc' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgalbedo) -!! where(avgalbedo /= spval)avgalbedo=avgalbedo/100. ! convert to fraction -!!$omp parallel do private(i,j) -! do j=jsta,jend -! do i=1,im -! if (avgalbedo(i,j) /= spval) avgalbedo(i,j) = avgalbedo(i,j) * 0.01 -! enddo -! enddo -! if(debugprint)print*,'sample ',VarName,' = ',avgalbedo(isa,jsa) - -! time averaged column cloud fractionusing nemsio - VarName='tcdc_aveclm' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgtcdc) -! where(avgtcdc /= spval)avgtcdc=avgtcdc/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgtcdc(i,j) /= spval) avgtcdc(i,j) = avgtcdc(i,j) * 0.01 - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',avgtcdc(isa,jsa) - -! GFS probably does not use zenith angle -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - Czen(i,j) = spval - CZMEAN(i,j) = SPVAL - enddo - enddo - -! maximum snow albedo in fraction using nemsio - VarName='snoalb' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,mxsnal) -! where(mxsnal /= spval)mxsnal=mxsnal/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (mxsnal(i,j) /= spval) mxsnal(i,j) = mxsnal(i,j) * 0.01 - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',mxsnal(isa,jsa) - -! 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) * T(I,J,LM) - Sigt4(i,j) = 5.67E-8 * TLMH * TLMH - End do - End do - -! TG is not used, skip it for now - -! GFS does not have inst cloud fraction for high, middle, and low cloud -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - cfrach(i,j) = spval - cfracl(i,j) = spval - cfracm(i,j) = spval - enddo - enddo - -! ave high cloud fraction using nemsio - VarName='tcdc_avehcl' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgcfrach) -! where(avgcfrach /= spval)avgcfrach=avgcfrach/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgcfrach(i,j) /= spval) avgcfrach(i,j) = avgcfrach(i,j) * 0.01 - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',avgcfrach(isa,jsa) - -! ave low cloud fraction using nemsio - VarName='tcdc_avelcl' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgcfracl) -! where(avgcfracl /= spval)avgcfracl=avgcfracl/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgcfracl(i,j) /= spval) avgcfracl(i,j) = avgcfracl(i,j) * 0.01 - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',avgcfracl(isa,jsa) - -! ave middle cloud fraction using nemsio - VarName='tcdc_avemcl' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgcfracm) -! where(avgcfracm /= spval)avgcfracm=avgcfracm/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgcfracm(i,j) /= spval) avgcfracm(i,j) = avgcfracm(i,j) * 0.01 - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',avgcfracm(isa,jsa) - -! inst convective cloud fraction using nemsio - VarName='tcdccnvcl' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,cnvcfr) -! where(cnvcfr /= spval)cnvcfr=cnvcfr/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (cnvcfr(i,j) /= spval) cnvcfr (i,j)= cnvcfr(i,j) * 0.01 - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',cnvcfr(isa,jsa) - -! slope type using nemsio - VarName='sltyp' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,buf) -!$omp parallel do private(i,j) - do j = jsta_2l, jend_2u - do i=1,im - if (buf(i,j) < spval) then - islope(i,j) = nint(buf(i,j)) - else - islope(i,j) = 0 - endif - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',islope(isa,jsa) - -! plant canopy sfc wtr in m - VarName='cnwat' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,cmc) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (cmc(i,j) /= spval) cmc(i,j) = cmc(i,j) * 0.001 - if (sm(i,j) /= 0.0) cmc(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',cmc(isa,jsa) - -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - grnflx(i,j) = spval ! GFS does not have inst ground heat flux - enddo - enddo - -! frozen precip fraction using nemsio - VarName='cpofp' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sr) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if(sr(i,j) /= spval) then -!set range within (0,1) - sr(i,j)=min(1.,max(0.,sr(i,j))) - endif - enddo - enddo - -! sea ice skin temperature - VarName='tisfc' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,ti) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sice(i,j) == spval .or. sice(i,j) == 0.) ti(i,j)=spval - enddo - enddo - -! vegetation fraction in fraction. using nemsio - VarName='veg' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,vegfrc) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (vegfrc(i,j) /= spval) then - vegfrc(i,j) = vegfrc(i,j) * 0.01 - else - vegfrc(i,j) = 0.0 - endif - enddo - enddo -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) vegfrc(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',vegfrc(isa,jsa) - -! 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='soill1' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sh2o(1,jsta_2l,1)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) sh2o(i,j,1) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,sh2o(isa,jsa,1) - - VarName='soill2' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sh2o(1,jsta_2l,2)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) sh2o(i,j,2) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,sh2o(isa,jsa,2) - - VarName='soill3' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sh2o(1,jsta_2l,3)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) sh2o(i,j,3) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,sh2o(isa,jsa,3) - - VarName='soill4' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sh2o(1,jsta_2l,4)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) sh2o(i,j,4) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,sh2o(isa,jsa,4) - -! volumetric soil moisture using nemsio - VarName='soilw1' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,smc(1,jsta_2l,1)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) smc(i,j,1) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,1) - - VarName='soilw2' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,smc(1,jsta_2l,2)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) smc(i,j,2) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,2) - - VarName='soilw3' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,smc(1,jsta_2l,3)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) smc(i,j,3) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,3) - - VarName='soilw4' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,smc(1,jsta_2l,4)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) smc(i,j,4) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,4) - -! soil temperature using nemsio - VarName='soilt1' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,stc(1,jsta_2l,1)) -! mask open water areas, combine with sea ice tmp -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,1) = spval - !if (sm(i,j) /= 0.0) stc(i,j,1) = spval - enddo - enddo - if(debugprint)print*,'sample l','stc',' = ',1,stc(isa,jsa,1) - - VarName='soilt2' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,stc(1,jsta_2l,2)) -! mask open water areas, combine with sea ice tmp -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,2) = spval - !if (sm(i,j) /= 0.0) stc(i,j,2) = spval - enddo - enddo - if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,2) - - VarName='soilt3' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,stc(1,jsta_2l,3)) -! mask open water areas, combine with sea ice tmp -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,3) = spval - !if (sm(i,j) /= 0.0) stc(i,j,3) = spval - enddo - enddo - if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,3) - - VarName='soilt4' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,stc(1,jsta_2l,4)) -! mask open water areas, combine with sea ice tmp -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,4) = spval - !if (sm(i,j) /= 0.0) stc(i,j,4) = spval - enddo - enddo - if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,4) - -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - acfrcv(i,j) = spval ! GFS does not output time averaged convective and strat cloud fraction, set acfrcv to spval, ncfrcv to 1 - ncfrcv(i,j) = 1.0 - acfrst(i,j) = spval ! GFS does not output time averaged cloud fraction, set acfrst to spval, ncfrst to 1 - ncfrst(i,j) = 1.0 - bgroff(i,j) = spval ! GFS does not have UNDERGROUND RUNOFF - rlwtoa(i,j) = spval ! GFS does not have inst model top outgoing longwave - enddo - enddo -! trdlw(i,j) = 6.0 - ardlw = 1.0 ! GFS incoming sfc longwave has been averaged over 6 hr bucket, set ARDLW to 1 - -! time averaged incoming sfc longwave - VarName='dlwrf_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,alwin) - -! inst incoming sfc longwave - VarName='dlwrf' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,rlwin) - -! time averaged outgoing sfc longwave - VarName='ulwrf_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,alwout) -! inst outgoing sfc longwave - VarName='ulwrf' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,radot) - -! where(alwout /= spval) alwout=-alwout ! CLDRAD puts a minus sign before gribbing -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (alwout(i,j) /= spval) alwout(i,j) = -alwout(i,j) - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,alwout(isa,jsa) - -! time averaged outgoing model top longwave using gfsio - VarName='ulwrf_avetoa' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,alwtoa) -! if(debugprint)print*,'sample l',VarName,' = ',1,alwtoa(isa,jsa) - -! GFS incoming sfc longwave has been averaged, set ARDLW to 1 - ardsw=1.0 -! trdsw=6.0 - -! time averaged incoming sfc shortwave - VarName='dswrf_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,aswin) -! if(debugprint)print*,'sample l',VarName,' = ',1,aswin(isa,jsa) - -! inst incoming sfc shortwave - VarName='dswrf' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,rswin) - -! inst incoming clear sky sfc shortwave - VarName='csdlf' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,rswinc) - -! time averaged incoming sfc uv-b using getgb - VarName='duvb_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,auvbin) -! if(debugprint)print*,'sample l',VarName,' = ',1,auvbin(isa,jsa) - -! time averaged incoming sfc clear sky uv-b using getgb - VarName='cduvb_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,auvbinc) -! if(debugprint)print*,'sample l',VarName,' = ',1,auvbinc(isa,jsa) - -! time averaged outgoing sfc shortwave using gfsio - VarName='uswrf_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,aswout) -! where(aswout /= spval) aswout=-aswout ! CLDRAD puts a minus sign before gribbing -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (aswout(i,j) /= spval) aswout(i,j) = -aswout(i,j) - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,aswout(isa,jsa) - -! inst outgoing sfc shortwave using gfsio - VarName='uswrf' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,rswout) - -! time averaged model top incoming shortwave - VarName='dswrf_avetoa' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,aswintoa) -! if(debugprint)print*,'sample l',VarName,' = ',1,aswintoa(isa,jsa) - -! time averaged model top outgoing shortwave - VarName='uswrf_avetoa' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,aswtoa) -! if(debugprint)print*,'sample l',VarName,' = ',1,aswtoa(isa,jsa) - -! time averaged surface sensible heat flux, multiplied by -1 because wrf model flux -! has reversed sign convention using gfsio - VarName='shtfl_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sfcshx) -! where (sfcshx /= spval)sfcshx=-sfcshx -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sfcshx(i,j) /= spval) sfcshx(i,j) = -sfcshx(i,j) - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,sfcshx(isa,jsa) - -! inst surface sensible heat flux - VarName='shtfl' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,twbs) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (twbs(i,j) /= spval) twbs(i,j) = -twbs(i,j) - enddo - enddo - -! 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 - VarName='lhtfl_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sfclhx) -! where (sfclhx /= spval)sfclhx=-sfclhx -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sfclhx(i,j) /= spval) sfclhx(i,j) = -sfclhx(i,j) - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,sfclhx(isa,jsa) - -! inst surface latent heat flux - VarName='lhtfl' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,qwbs) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (qwbs(i,j) /= spval) qwbs(i,j) = -qwbs(i,j) - enddo - enddo - - if(me==0)print*,'rdaod= ',rdaod -! inst aod550 optical depth - if(rdaod) then - VarName='aod550' - call read_netcdf_2d_scatter(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,aod550) - - VarName='du_aod550' - call read_netcdf_2d_scatter(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,du_aod550) - - VarName='ss_aod550' - call read_netcdf_2d_scatter(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,ss_aod550) - - VarName='su_aod550' - call read_netcdf_2d_scatter(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,su_aod550) - - VarName='oc_aod550' - call read_netcdf_2d_scatter(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,oc_aod550) - - VarName='bc_aod550' - call read_netcdf_2d_scatter(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,bc_aod550) - end if - -! time averaged ground heat flux using nemsio - VarName='gflux_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,subshx) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) subshx(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,subshx(isa,jsa) - -! inst ground heat flux using nemsio - VarName='gflux' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,grnflx) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) grnflx(i,j) = spval - enddo - enddo - -! time averaged zonal momentum flux using gfsio - VarName='uflx_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sfcux) -! if(debugprint)print*,'sample l',VarName,' = ',1,sfcux(isa,jsa) - -! time averaged meridional momentum flux using nemsio - VarName='vflx_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sfcvx) -! if(debugprint)print*,'sample l',VarName,' = ',1,sfcvx(isa,jsa) - -! dong read in inst surface flux -! inst zonal momentum flux using gfsio -! VarName='uflx' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sfcuxi) -! if(debugprint)print*,'sample l',VarName,' = ',1,sfcuxi(isa,jsa) - -! inst meridional momentum flux using nemsio -! VarName='vflx' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sfcvxi) -! if(debugprint)print*,'sample l',VarName,' = ',1,sfcvxi(isa,jsa) - - -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - sfcuvx(i,j) = spval ! GFS does not use total momentum flux - enddo - enddo - -! time averaged zonal gravity wave stress using nemsio - VarName='u-gwd_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,gtaux) -! if(debugprint)print*,'sample l',VarName,' = ',1,gtaux(isa,jsa) - -! time averaged meridional gravity wave stress using getgb - VarName='v-gwd_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,gtauy) -! if(debugprint)print*,'sample l',VarName,' = ',1,gtauy(isa,jsa) - -! time averaged accumulated potential evaporation - VarName='pevpr_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgpotevp) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) avgpotevp(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,potevp(isa,jsa) - -! inst potential evaporation - VarName='pevpr' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,potevp) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) potevp(i,j) = spval - enddo - enddo - - do l=1,lm -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im -! GFS does not have temperature tendency due to long wave radiation - rlwtt(i,j,l) = spval -! GFS does not have temperature tendency due to short wave radiation - rswtt(i,j,l) = spval -! GFS does not have temperature tendency due to latent heating from convection - tcucn(i,j,l) = spval - tcucns(i,j,l) = spval -! GFS does not have temperature tendency due to latent heating from grid scale - train(i,j,l) = spval - enddo - enddo - enddo - -! set avrain to 1 - avrain=1.0 - avcnvc=1.0 - theat=6.0 ! just in case GFS decides to output T tendency - -! 10 m u using nemsio - VarName='ugrd10m' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,u10) - - do j=jsta,jend - do i=1,im - u10h(i,j)=u10(i,j) - end do - end do -! if(debugprint)print*,'sample l',VarName,' = ',1,u10(isa,jsa) - -! 10 m v using gfsio - VarName='vgrd10m' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,v10) - - do j=jsta,jend - do i=1,im - v10h(i,j)=v10(i,j) - end do - end do -! if(debugprint)print*,'sample l',VarName,' = ',1,v10(isa,jsa) - -! vegetation type, it's in GFS surface file, hopefully will merge into gfsio soon - VarName='vtype' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,buf) -! where (buf /= spval) -! ivgtyp=nint(buf) -! elsewhere -! ivgtyp=0 !need to feed reasonable value to crtm -! end where -!$omp parallel do private(i,j) - do j = jsta_2l, jend_2u - do i=1,im - if (buf(i,j) < spval) then - ivgtyp(i,j) = nint(buf(i,j)) - else - ivgtyp(i,j) = 0 - endif - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,ivgtyp(isa,jsa) - -! soil type, it's in GFS surface file, hopefully will merge into gfsio soon - VarName='sotyp' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,buf) - VcoordName='sfc' - l=1 -!$omp parallel do private(i,j) - do j = jsta_2l, jend_2u - do i=1,im - if (buf(i,j) < spval) then - isltyp(i,j) = nint(buf(i,j)) - else - isltyp(i,j) = 0 !need to feed reasonable value to crtm - endif - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,isltyp(isa,jsa) - -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - smstav(i,j) = spval ! GFS does not have soil moisture availability -! smstot(i,j) = spval ! GFS does not have total soil moisture - sfcevp(i,j) = spval ! GFS does not have accumulated surface evaporation - acsnow(i,j) = spval ! GFS does not have averaged accumulated snow - acsnom(i,j) = spval ! GFS does not have snow melt -! sst(i,j) = spval ! GFS does not have sst???? - thz0(i,j) = ths(i,j) ! GFS does not have THZ0, use THS to substitute - qz0(i,j) = spval ! GFS does not output humidity at roughness length - uz0(i,j) = spval ! GFS does not output u at roughness length - vz0(i,j) = spval ! GFS does not output humidity at roughness length - enddo - enddo - do l=1,lm -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - EL_PBL(i,j,l) = spval ! GFS does not have mixing length - exch_h(i,j,l) = spval ! GFS does not output exchange coefficient - enddo - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,thz0(isa,jsa) - -! 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' -! l=1 -! if(debugprint)print*,'sample l',VarName,' = ',1,ptop(isa,jsa) - VarName='prescnvclt' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,ptop) - - -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - htop(i,j) = spval - if(ptop(i,j) <= 0.0) ptop(i,j) = spval - enddo - enddo - do j=jsta,jend - do i=1,im - if(ptop(i,j) < spval)then - do l=1,lm - if(ptop(i,j) <= pmid(i,j,l))then - htop(i,j) = l -! if(i==ii .and. j==jj)print*,'sample ptop,pmid pmid-1,pint= ', & -! ptop(i,j),pmid(i,j,l),pmid(i,j,l-1),pint(i,j,l),htop(i,j) - exit - end if - end do - end if - end do - end do - -! 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='prescnvclb' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,pbot) -! if(debugprint)print*,'sample l',VarName,VcoordName,' = ',1,pbot(isa,jsa) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - hbot(i,j) = spval - if(pbot(i,j) <= 0.0) pbot(i,j) = spval - enddo - enddo - do j=jsta,jend - do i=1,im -! if(.not.lb(i,j))print*,'false bitmask for pbot at ' -! + ,i,j,pbot(i,j) - if(pbot(i,j) < spval)then - do l=lm,1,-1 - if(pbot(i,j) >= pmid(i,j,l)) then - hbot(i,j) = l -! if(i==ii .and. j==jj)print*,'sample pbot,pmid= ', & -! pbot(i,j),pmid(i,j,l),hbot(i,j) - exit - end if - end do - end if - end do - end do - if(debugprint)print*,'sample hbot = ',hbot(isa,jsa) -! retrieve time averaged low cloud top pressure using nemsio - VarName='pres_avelct' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,ptopl) -! if(debugprint)print*,'sample l',VarName,' = ',1,ptopl(isa,jsa) - -! retrieve time averaged low cloud bottom pressure using nemsio - VarName='pres_avelcb' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,pbotl) -! if(debugprint)print*,'sample l',VarName,' = ',1,pbotl(isa,jsa) - -! retrieve time averaged low cloud top temperature using nemsio - VarName='tmp_avelct' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,Ttopl) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopl(isa,jsa) - -! retrieve time averaged middle cloud top pressure using nemsio - VarName='pres_avemct' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,ptopm) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptopm(isa,jsa) - -! retrieve time averaged middle cloud bottom pressure using nemsio - VarName='pres_avemcb' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,pbotm) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pbotm(isa,jsa) - -! retrieve time averaged middle cloud top temperature using nemsio - VarName='tmp_avemct' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,Ttopm) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopm(isa,jsa) - -! retrieve time averaged high cloud top pressure using nemsio ********* - VarName='pres_avehct' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,ptoph) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptoph(isa,jsa) - -! retrieve time averaged high cloud bottom pressure using nemsio - VarName='pres_avehcb' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,pboth) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pboth(isa,jsa) - -! retrieve time averaged high cloud top temperature using nemsio - VarName='tmp_avehct' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,Ttoph) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttoph(isa,jsa) - -! retrieve boundary layer cloud cover using nemsio - VarName='tcdc_avebndcl' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,pblcfr) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,pblcfr(isa,jsa) -! where (pblcfr /= spval)pblcfr=pblcfr/100. ! convert to fraction -!$omp parallel do private(i,j) - do j = jsta_2l, jend_2u - do i=1,im - if (pblcfr(i,j) < spval) pblcfr(i,j) = pblcfr(i,j) * 0.01 - enddo - enddo - -! retrieve cloud work function - VarName='cwork_aveclm' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,cldwork) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,cldwork(isa,jsa) - -! accumulated total (base+surface) runoff - VarName='watr_acc' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,runoff) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) runoff(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,runoff(isa,jsa) - -! retrieve shelter max temperature using nemsio - VarName='tmax_max2m' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,maxtshltr) - -! retrieve shelter min temperature using nemsio - VarName='tmin_min2m' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,mintshltr) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', & -! 1,mintshltr(im/2,(jsta+jend)/2) - -! retrieve shelter max RH -! VarName='rh02max' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,maxrhshltr) - -! retrieve shelter min temperature using nemsio -! VarName='rh02min' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,minrhshltr) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', & -! 1,mintshltr(im/2,(jsta+jend)/2) - -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - MAXRHSHLTR(i,j) = SPVAL - MINRHSHLTR(i,j) = SPVAL - enddo - enddo - -! retrieve ice thickness using nemsio - VarName='icetk' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,dzice) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,dzice(isa,jsa) - -! retrieve wilting point using nemsio - VarName='wilt' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,smcwlt) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) smcwlt(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,smcwlt(isa,jsa) - -! retrieve sunshine duration using nemsio - VarName='sunsd_acc' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,suntime) - -! retrieve field capacity using nemsio - VarName='fldcp' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,fieldcapa) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) fieldcapa(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,fieldcapa(isa,jsa) - -! retrieve time averaged surface visible beam downward solar flux - VarName='vbdsf_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avisbeamswin) - VcoordName='sfc' - l=1 - -! retrieve time averaged surface visible diffuse downward solar flux - VarName='vddsf_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avisdiffswin) - -! retrieve time averaged surface near IR beam downward solar flux - VarName='nbdsf_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,airbeamswin) - -! retrieve time averaged surface near IR diffuse downward solar flux - VarName='nddsf_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,airdiffswin) - -! retrieve time averaged surface clear sky outgoing LW - VarName='csulf' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,alwoutc) - -! retrieve time averaged TOA clear sky outgoing LW - VarName='csulftoa' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,alwtoac) - -! retrieve time averaged surface clear sky outgoing SW - VarName='csusf' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,aswoutc) - -! retrieve time averaged TOA clear sky outgoing LW - VarName='csusftoa' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,aswtoac) - -! retrieve time averaged surface clear sky incoming LW - VarName='csdlf' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,alwinc) - -! retrieve time averaged surface clear sky incoming SW - VarName='csdsf' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,aswinc) - -! retrieve shelter max specific humidity using nemsio - VarName='spfhmax_max2m' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,maxqshltr) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', -! 1,maxqshltr(isa,jsa) - -! retrieve shelter min temperature using nemsio - VarName='spfhmin_min2m' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,minqshltr) - -! retrieve storm runoff using nemsio - VarName='ssrun_acc' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,SSROFF) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) ssroff(i,j) = spval - enddo - enddo - -! retrieve direct soil evaporation - VarName='evbs_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgedir) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) avgedir(i,j) = spval - enddo - enddo - -! retrieve CANOPY WATER EVAP - VarName='evcw_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgecan) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) avgecan(i,j) = spval - enddo - enddo - -! retrieve PLANT TRANSPIRATION - VarName='trans_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgetrans) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) avgetrans(i,j) = spval - enddo - enddo - -! retrieve snow sublimation - VarName='sbsno_ave' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgesnow) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j)==1.0 .and. sice(i,j)==0.) avgesnow(i,j)=spval - enddo - enddo - -! retrive total soil moisture - VarName='soilm' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,smstot) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) smstot(i,j) = spval - enddo - enddo - -! retrieve snow phase change heat flux - VarName='snohf' - call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & - ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,snopcx) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) snopcx(i,j) = spval - enddo - enddo - -! GFS does not have deep convective cloud top and bottom fields - -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - HTOPD(i,j) = SPVAL - HBOTD(i,j) = SPVAL - HTOPS(i,j) = SPVAL - HBOTS(i,j) = SPVAL - CUPPT(i,j) = SPVAL - enddo - enddo - -! done with flux file, close it for now - Status=nf90_close(ncid2d) -! deallocate(tmp,recname,reclevtyp,reclev) - -! 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 == 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 -! - -! generate look up table for lifted parcel calculations - - THL = 210. - PLQ = 70000. - pt_TBL = 10000. ! this is for 100 hPa added by Moorthi - - CALL TABLE(PTBL,TTBL,PT_TBL, & - RDQ,RDTH,RDP,RDTHE,PL,THL,QS0,SQS,STHE,THE0) - - CALL TABLEQ(TTBLQ,RDPQ,RDTHEQ,PLQ,THL,STHEQ,THE0Q) - -! -! - IF(ME == 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 -! -!$omp parallel do private(l) - DO L = 1,LSM - ALSL(L) = LOG(SPL(L)) - END DO -! -!HC WRITE IGDS OUT FOR WEIGHTMAKER TO READ IN AS KGDSIN - if(me == 0)then - print*,'writing out igds' - igdout = 110 -! open(igdout,file='griddef.out',form='unformatted' -! + ,status='unknown') - if(maptype == 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 == 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 - ! Note: The calculation of the map scale factor at the standard - ! lat/lon and the PSMAPF - ! Get map factor at 60 degrees (N or S) for PS projection, which will - ! be needed to correctly define the DX and DY values in the GRIB GDS - if (TRUELAT1 < 0.) THEN - LAT = -60. - else - LAT = 60. - end if - - CALL MSFPS (LAT,TRUELAT1*0.001,PSMAPF) - - ELSE IF(MAPTYPE == 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 == 0 .OR. MAPTYPE == 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 -! -! - - RETURN - END - diff --git a/sorc/ncep_post.fd/INITPOST_GFS_NETCDF_PARA.f b/sorc/ncep_post.fd/INITPOST_GFS_NETCDF_PARA.f deleted file mode 100644 index aee28bb0b..000000000 --- a/sorc/ncep_post.fd/INITPOST_GFS_NETCDF_PARA.f +++ /dev/null @@ -1,2672 +0,0 @@ -!> @file -! . . . -!> SUBPROGRAM: INITPOST_GFS_NETCDF_PARA INITIALIZE POST FOR RUN -!! PRGRMMR: Wen Meng DATE: 2020-02-04 -!! -!! ABSTRACT: THIS ROUTINE INITIALIZES CONSTANTS AND -!! VARIABLES AT THE START OF GFS MODEL OR POST -!! PROCESSOR RUN. -!! -!! REVISION HISTORY -!! 2020-02-04 W Meng start from INITPOST_GFS_NETCDF.f -!! 2021-03-11 Bo Cui change local arrays to dimension (im,jsta:jend) -!! -!! USAGE: CALL INITPOST_GFS_NETCDF_PARA -!! 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 -!! - SUBROUTINE INITPOST_GFS_NETCDF_PARA(ncid3d) - - - use netcdf - use vrbls4d, only: dust, SALT, SUSO, SOOT, WASO, PP25, PP10 - use vrbls3d, only: t, q, uh, vh, pmid, pint, alpint, dpres, zint, zmid, o3, & - qqr, qqs, cwm, qqi, qqw, omga, rhomid, q2, cfr, rlwtt, rswtt, tcucn, & - tcucns, train, el_pbl, exch_h, vdifftt, vdiffmois, dconvmois, nradtt, & - o3vdiff, o3prod, o3tndy, mwpv, unknown, vdiffzacce, zgdrag,cnvctummixing, & - vdiffmacce, mgdrag, cnvctvmmixing, ncnvctcfrac, cnvctumflx, cnvctdmflx, & - cnvctzgdrag, sconvmois, cnvctmgdrag, cnvctdetmflx, duwt, duem, dusd, dudp, & - wh, qqg, ref_10cm - use vrbls2d, only: f, pd, fis, pblh, ustar, z0, ths, qs, twbs, qwbs, avgcprate, & - cprate, avgprec, prec, lspa, sno, si, cldefi, th10, q10, tshltr, pshltr, & - tshltr, albase, avgalbedo, avgtcdc, czen, czmean, mxsnal, landfrac, radot, sigt4, & - cfrach, cfracl, cfracm, avgcfrach, qshltr, avgcfracl, avgcfracm, cnvcfr, & - islope, cmc, grnflx, vegfrc, acfrcv, ncfrcv, acfrst, ncfrst, ssroff, & - bgroff, rlwin, rlwtoa, cldwork, alwin, alwout, alwtoa, rswin, rswinc, & - rswout, aswin, auvbin, auvbinc, aswout, aswtoa, sfcshx, sfclhx, subshx, & - snopcx, sfcux, sfcvx, sfcuxi, sfcvxi, sfcuvx, gtaux, gtauy, potevp, u10, v10, smstav, & - smstot, ivgtyp, isltyp, sfcevp, sfcexc, acsnow, acsnom, sst, thz0, qz0, & - uz0, vz0, ptop, htop, pbot, hbot, ptopl, pbotl, ttopl, ptopm, pbotm, ttopm, & - ptoph, pboth, pblcfr, ttoph, runoff, tecan, tetran, tedir, twa, maxtshltr, & - mintshltr, maxrhshltr, fdnsst, & - minrhshltr, dzice, smcwlt, suntime, fieldcapa, htopd, hbotd, htops, hbots, & - cuppt, dusmass, ducmass, dusmass25, ducmass25, aswintoa,rel_vort_maxhy1, & - maxqshltr, minqshltr, acond, sr, u10h, v10h,refd_max, w_up_max, w_dn_max, & - up_heli_max,up_heli_min,up_heli_max03,up_heli_min03,rel_vort_max01,u10max, v10max, & - avgedir,avgecan,paha,pahi,avgetrans,avgesnow,avgprec_cont,avgcprate_cont,rel_vort_max, & - avisbeamswin,avisdiffswin,airbeamswin,airdiffswin,refdm10c_max,wspd10max, & - alwoutc,alwtoac,aswoutc,aswtoac,alwinc,aswinc,avgpotevp,snoavg, & - ti,aod550,du_aod550,ss_aod550,su_aod550,oc_aod550,bc_aod550 - use soil, only: sldpth, sh2o, smc, stc - use masks, only: lmv, lmh, htm, vtm, gdlat, gdlon, dx, dy, hbm2, sm, sice - use physcons_post, only: grav => con_g, fv => con_fvirt, rgas => con_rd, & - eps => con_eps, epsm1 => con_epsm1 - use params_mod, only: erad, dtr, tfrz, h1, d608, rd, p1000, capa,pi - use lookup_mod, only: thl, plq, ptbl, ttbl, rdq, rdth, rdp, rdthe, pl, qs0, sqs, sthe, & - ttblq, rdpq, rdtheq, stheq, the0q, the0 - use ctlblk_mod, only: me, mpi_comm_comp, icnt, idsp, jsta, jend, ihrst, idat, sdat, ifhr, & - ifmin, filename, tprec, tclod, trdlw, trdsw, tsrfc, tmaxmin, td3d, restrt, sdat, & - jend_m, imin, imp_physics, dt, spval, pdtop, pt, qmin, nbin_du, nphs, dtq2, ardlw,& - ardsw, asrfc, avrain, avcnvc, theat, gdsdegr, spl, lsm, alsl, im, jm, im_jm, lm, & - jsta_2l, jend_2u, nsoil, lp1, icu_physics, ivegsrc, novegtype, nbin_ss, nbin_bc, & - nbin_oc, nbin_su, gocart_on, pt_tbl, hyb_sigp, filenameFlux, fileNameAER,rdaod - use gridspec_mod, only: maptype, gridtype, latstart, latlast, lonstart, lonlast, cenlon, & - dxval, dyval, truelat2, truelat1, psmapf, cenlat,lonstartv, lonlastv, cenlonv, & - latstartv, latlastv, cenlatv,latstart_r,latlast_r,lonstart_r,lonlast_r - - use upp_physics, only: fpvsnew -!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none -! -! type(nemsio_gfile) :: nfile,ffile,rfile - integer,parameter :: nvar2d=48 -! character(nemsio_charkind) :: name2d(nvar2d) - integer :: nvar3d, numDims -! character(nemsio_charkind), allocatable :: name3din(:), name3dout(:) -! character(nemsio_charkind) :: varname,levtype -! -! 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 -! -! 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_. - - real, parameter :: gravi = 1.0/grav - character(len=20) :: VarName, VcoordName - integer :: Status, fldsize, fldst, recn, recn_vvel - 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,IOOMG,IOALL -! logical, parameter :: debugprint = .true., zerout = .false. - logical, parameter :: debugprint = .false., zerout = .false. - logical :: convert_rad_to_deg=.false. - CHARACTER*32 varcharval -! CHARACTER*40 CONTRL,FILALL,FILMST,FILTMP,FILTKE,FILUNV,FILCLD,FILRAD,FILSFC - CHARACTER*4 RESTHR - CHARACTER FNAME*255,ENVAR*50 - INTEGER IDATE(8),JDATE(8),JPDS(200),JGDS(200),KPDS(200),KGDS(200) -! LOGICAL*1 LB(IM,JM) -! -! INCLUDE COMMON BLOCKS. -! -! DECLARE VARIABLES. -! -! REAL fhour - integer nfhour ! forecast hour from nems io file - integer fhzero !bucket - real dtp !physics time step - REAL RINC(5) - - REAL DUMMY(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, & - nframed2,iunitd3d,ierr,idum,iret,nrec,idrt - integer ncid3d,ncid2d,varid,nhcas - real TSTART,TLMH,TSPH,ES,FACT,soilayert,soilayerb,zhour,dum, & - tvll,pmll,tv, tx1, tx2 - - character*20,allocatable :: recname(:) - integer, allocatable :: reclev(:), kmsk(:,:) - real, allocatable :: glat1d(:), glon1d(:), qstl(:) - real, allocatable :: wrk1(:,:), wrk2(:,:) - real, allocatable :: p2d(:,:), t2d(:,:), q2d(:,:), & - qs2d(:,:), cw2d(:,:), cfr2d(:,:) - real(kind=4),allocatable :: vcoord4(:,:,:) - real, dimension(lm+1) :: ak5, bk5 - real*8, allocatable :: pm2d(:,:), pi2d(:,:) - real, allocatable :: tmp(:) - real :: buf(im,jsta_2l:jend_2u) - real :: buf3d(im,jsta_2l:jend_2u,lm) - -! real buf(im,jsta_2l:jend_2u),bufsoil(im,nsoil,jsta_2l:jend_2u) & -! ,buf3d(im,jsta_2l:jend_2u,lm),buf3d2(im,lp1,jsta_2l:jend_2u) - - real LAT - integer isa, jsa, latghf, jtem, idvc, idsl, nvcoord, ip1, nn, npass - - integer, parameter :: npass2=5, npass3=30 - real, parameter :: third=1.0/3.0 - INTEGER, DIMENSION(2) :: ij4min, ij4max - REAL :: omgmin, omgmax - real, allocatable :: d2d(:,:), u2d(:,:), v2d(:,:), omga2d(:,:) - REAL, ALLOCATABLE :: ps2d(:,:),psx2d(:,:),psy2d(:,:) - real, allocatable :: div3d(:,:,:) - real(kind=4),allocatable :: vcrd(:,:) - real :: dum_const - -!*********************************************************************** -! START INIT HERE. -! - WRITE(6,*)'INITPOST: ENTER INITPOST_GFS_NETCDF_PARA' - WRITE(6,*)'me=',me, & - 'jsta_2l=',jsta_2l,'jend_2u=', & - jend_2u,'im=',im -! - isa = im / 2 - jsa = (jsta+jend) / 2 - -!$omp parallel do private(i,j) - do j = jsta_2l, jend_2u - do i=1,im - buf(i,j) = spval - enddo - enddo - - Status=nf90_get_att(ncid3d,nf90_global,'ak',ak5) - if(Status/=0)then - print*,'ak not found; assigning missing value' - ak5=spval - else - if(me==0)print*,'ak5= ',ak5 - end if - Status=nf90_get_att(ncid3d,nf90_global,'idrt',idrt) - if(Status/=0)then - print*,'idrt not in netcdf file,reading grid' - Status=nf90_get_att(ncid3d,nf90_global,'grid',varcharval) - if(Status/=0)then - print*,'idrt and grid not in netcdf file, set default to latlon' - idrt=0 - MAPTYPE=0 - else - if(trim(varcharval)=='rotated_latlon')then - MAPTYPE=207 - idrt=207 - Status=nf90_get_att(ncid3d,nf90_global,'cen_lon',dum_const) - if(Status/=0)then - print*,'cen_lon not found; assigning missing value' - cenlon=spval - else - if(dum_const<0.)then - cenlon=nint((dum_const+360.)*gdsdegr) - else - cenlon=dum_const*gdsdegr - end if - end if - Status=nf90_get_att(ncid3d,nf90_global,'cen_lat',dum_const) - if(Status/=0)then - print*,'cen_lat not found; assigning missing value' - cenlat=spval - else - cenlat=dum_const*gdsdegr - end if - - Status=nf90_get_att(ncid3d,nf90_global,'lon1',dum_const) - if(Status/=0)then - print*,'lonstart_r not found; assigning missing value' - lonstart_r=spval - else - if(dum_const<0.)then - lonstart_r=nint((dum_const+360.)*gdsdegr) - else - lonstart_r=dum_const*gdsdegr - end if - end if - Status=nf90_get_att(ncid3d,nf90_global,'lat1',dum_const) - if(Status/=0)then - print*,'latstart_r not found; assigning missing value' - latstart_r=spval - else - latstart_r=dum_const*gdsdegr - end if - - Status=nf90_get_att(ncid3d,nf90_global,'lon2',dum_const) - if(Status/=0)then - print*,'lonlast_r not found; assigning missing value' - lonlast_r=spval - else - if(dum_const<0.)then - lonlast_r=nint((dum_const+360.)*gdsdegr) - else - lonlast_r=dum_const*gdsdegr - end if - end if - Status=nf90_get_att(ncid3d,nf90_global,'lat2',dum_const) - if(Status/=0)then - print*,'latlast_r not found; assigning missing value' - latlast_r=spval - else - latlast_r=dum_const*gdsdegr - end if - - Status=nf90_get_att(ncid3d,nf90_global,'dlon',dum_const) - if(Status/=0)then - print*,'dlmd not found; assigning missing value' - dxval=spval - else - dxval=dum_const*gdsdegr - end if - Status=nf90_get_att(ncid3d,nf90_global,'dlat',dum_const) - if(Status/=0)then - print*,'dphd not found; assigning missing value' - dyval=spval - else - dyval=dum_const*gdsdegr - end if - - print*,'lonstart,latstart,cenlon,cenlat,dyval,dxval', & - lonstart,latstart,cenlon,cenlat,dyval,dxval - -! Jili Dong add support for regular lat lon (2019/03/22) start - else if(trim(varcharval)=='latlon')then - MAPTYPE=0 - idrt=0 - - Status=nf90_get_att(ncid3d,nf90_global,'lon1',dum_const) - if(Status/=0)then - print*,'lonstart not found; assigning missing value' - lonstart=spval - else - if(dum_const<0.)then - lonstart=nint((dum_const+360.)*gdsdegr) - else - lonstart=dum_const*gdsdegr - end if - end if - Status=nf90_get_att(ncid3d,nf90_global,'lat1',dum_const) - if(Status/=0)then - print*,'latstart not found; assigning missing value' - latstart=spval - else - latstart=dum_const*gdsdegr - end if - - Status=nf90_get_att(ncid3d,nf90_global,'lon2',dum_const) - if(Status/=0)then - print*,'lonlast not found; assigning missing value' - lonlast=spval - else - if(dum_const<0.)then - lonlast=nint((dum_const+360.)*gdsdegr) - else - lonlast=dum_const*gdsdegr - end if - end if - Status=nf90_get_att(ncid3d,nf90_global,'lat2',dum_const) - if(Status/=0)then - print*,'latlast not found; assigning missing value' - latlast=spval - else - latlast=dum_const*gdsdegr - end if - - Status=nf90_get_att(ncid3d,nf90_global,'dlon',dum_const) - if(Status/=0)then - print*,'dlmd not found; assigning missing value' - dxval=spval - else - dxval=dum_const*gdsdegr - end if - Status=nf90_get_att(ncid3d,nf90_global,'dlat',dum_const) - if(Status/=0)then - print*,'dphd not found; assigning missing value' - dyval=spval - else - dyval=dum_const*gdsdegr - end if - - print*,'lonstart,latstart,dyval,dxval', & - lonstart,lonlast,latstart,latlast,dyval,dxval - -! Jili Dong add support for regular lat lon (2019/03/22) end - - else if(trim(varcharval)=='gaussian')then - MAPTYPE=4 - idrt=4 - else ! setting default maptype - MAPTYPE=0 - idrt=0 - end if - end if !end reading grid - end if !end reading idrt - if(me==0)print*,'idrt MAPTYPE= ',idrt,MAPTYPE -! STEP 1. READ MODEL OUTPUT FILE -! -! -!*** -! -! LMH and 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 - -! 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 - - Status=nf90_get_att(ncid3d,nf90_global,'nhcas',nhcas) - if(Status/=0)then - print*,'nhcas not in netcdf file, set default to nonhydro' - nhcas=0 - end if - if(me==0)print*,'nhcas= ',nhcas - if (nhcas == 0 ) then !non-hydrostatic case - nrec=15 - allocate (recname(nrec)) - recname=[character(len=20) :: 'ugrd','vgrd','spfh','tmp','o3mr', & - 'presnh','dzdt', 'clwmr','dpres', & - 'delz','icmr','rwmr', & - 'snmr','grle','cld_amt'] - else - nrec=8 - allocate (recname(nrec)) - recname=[character(len=20) :: 'ugrd','vgrd','tmp','spfh','o3mr', & - 'hypres', 'clwmr','dpres'] - endif - -! write(0,*)'nrec=',nrec - !allocate(recname(nrec),reclevtyp(nrec),reclev(nrec)) - allocate(glat1d(jm),glon1d(im)) - allocate(vcoord4(lm+1,3,2)) - -! hardwire idate for now -! idate=(/2017,08,07,00,0,0,0,0/) -! get cycle start time - Status=nf90_inq_varid(ncid3d,'time',varid) - if(Status/=0)then - print*,'time not in netcdf file, stopping' - STOP 1 - else - Status=nf90_get_att(ncid3d,varid,'units',varcharval) - if(Status/=0)then - print*,'time unit not available' - else - print*,'time unit read from netcdf file= ',varcharval -! assume use hours as unit -! idate_loc=index(varcharval,'since')+6 - read(varcharval,101)idate(1),idate(2),idate(3),idate(4),idate(5) - end if -! Status=nf90_inquire_dimension(ncid3d,varid,len=ntimes) -! allocate(fhours(ntimes)) -! status = nf90_inq_varid(ncid3d,varid,fhours) -! Status=nf90_get_var(ncid3d,varid,nfhour,start=(/1/), & -! count=(/1/)) -! if(Status/=0)then -! print*,'forecast hour not in netcdf file, stopping' -! STOP 1 -! end if - end if - 101 format(T13,i4,1x,i2,1x,i2,1x,i2,1x,i2) - print*,'idate= ',idate(1:5) -! get longitude - Status=nf90_inq_varid(ncid3d,'grid_xt',varid) - Status=nf90_inquire_variable(ncid3d,varid,ndims = numDims) - if(debugprint)print*,'number of dim for gdlon ',numDims - if(numDims==1)then - Status=nf90_get_var(ncid3d,varid,glon1d) - do j=jsta,jend - do i=1,im - gdlon(i,j) = real(glon1d(i),kind=4) - end do - end do - lonstart = nint(glon1d(1)*gdsdegr) - lonlast = nint(glon1d(im)*gdsdegr) - dxval = nint(abs(glon1d(1)-glon1d(2))*gdsdegr) - else if(numDims==2)then - Status=nf90_get_var(ncid3d,varid,dummy) - if(maxval(abs(dummy))<2.0*pi)convert_rad_to_deg=.true. - if(convert_rad_to_deg)then - do j=jsta,jend - do i=1,im - gdlon(i,j) = real(dummy(i,j),kind=4)*180./pi - end do - end do - else - do j=jsta,jend - do i=1,im - gdlon(i,j) = real(dummy(i,j),kind=4) - end do - end do - end if - if(convert_rad_to_deg)then - lonstart = nint(dummy(1,1)*gdsdegr)*180./pi - lonlast = nint(dummy(im,jm)*gdsdegr)*180./pi - dxval = nint(abs(dummy(1,1)-dummy(2,1))*gdsdegr)*180./pi - else - lonstart = nint(dummy(1,1)*gdsdegr) - lonlast = nint(dummy(im,jm)*gdsdegr) - dxval = nint(abs(dummy(1,1)-dummy(2,1))*gdsdegr) - end if - -! Jili Dong add support for regular lat lon (2019/03/22) start - if (MAPTYPE == 0) then - if(lonstart<0.)then - lonstart=lonstart+360.*gdsdegr - end if - if(lonlast<0.)then - lonlast=lonlast+360.*gdsdegr - end if - end if -! Jili Dong add support for regular lat lon (2019/03/22) end - - end if - print*,'lonstart,lonlast,dxval ',lonstart,lonlast,dxval -! get latitude - Status=nf90_inq_varid(ncid3d,'grid_yt',varid) - Status=nf90_inquire_variable(ncid3d,varid,ndims = numDims) - if(debugprint)print*,'number of dim for gdlat ',numDims - if(numDims==1)then - Status=nf90_get_var(ncid3d,varid,glat1d) - do j=jsta,jend - do i=1,im - gdlat(i,j) = real(glat1d(j),kind=4) - end do - end do - latstart = nint(glat1d(1)*gdsdegr) - latlast = nint(glat1d(jm)*gdsdegr) - dyval = nint(abs(glat1d(1)-glat1d(2))*gdsdegr) - else if(numDims==2)then - Status=nf90_get_var(ncid3d,varid,dummy) - if(maxval(abs(dummy))1000.)print*,'bad T ',t(i,j,l) - enddo - enddo - enddo - call read_netcdf_3d_para(ncid3d,im,jm,jsta,jsta_2l,jend,jend_2u, & - spval,recname(11),qqi(1,jsta_2l,1),lm) - call read_netcdf_3d_para(ncid3d,im,jm,jsta,jsta_2l,jend,jend_2u, & - spval,recname(12),qqr(1,jsta_2l,1),lm) - call read_netcdf_3d_para(ncid3d,im,jm,jsta,jsta_2l,jend,jend_2u, & - spval,recname(13),qqs(1,jsta_2l,1),lm) - call read_netcdf_3d_para(ncid3d,im,jm,jsta,jsta_2l,jend,jend_2u, & - spval,recname(14),qqg(1,jsta_2l,1),lm) - call read_netcdf_3d_para(ncid3d,im,jm,jsta,jsta_2l,jend,jend_2u, & - spval,recname(15),cfr(1,jsta_2l,1),lm) - -! calculate CWM from FV3 output - do l=1,lm - do j=jsta,jend - do i=1,im - cwm(i,j,l)=qqg(i,j,l)+qqs(i,j,l)+qqr(i,j,l)+qqi(i,j,l)+qqw(i,j,l) - enddo - enddo - if(debugprint)print*,'sample l,t,q,u,v,w,= ',isa,jsa,l & - ,t(isa,jsa,l),q(isa,jsa,l),uh(isa,jsa,l),vh(isa,jsa,l) & - ,wh(isa,jsa,l) - if(debugprint)print*,'sample l cwm for FV3',l, & - cwm(isa,jsa,l) - end do - -! surface pressure - VarName='pressfc' - call read_netcdf_2d_para(ncid3d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,pint(1,jsta_2l,lp1)) - do j=jsta,jend - do i=1,im -! if(pint(i,j,lp1)>1.0E6 .or. pint(1,jsta_2l,lp1)<50000.) & -! print*,'bad psfc ',i,j,pint(i,j,lp1) - end do - end do - if(debugprint)print*,'sample ',VarName,' =',pint(isa,jsa,lp1) - - pt = ak5(1) - - do j=jsta,jend - do i=1,im - pint(i,j,1)= pt - end do - end do - - do l=2,lp1 - do j=jsta,jend - do i=1,im - pint(i,j,l) = pint(i,j,l-1) + dpres(i,j,l-1) - enddo - enddo -! if (me == 0) print*,'sample model pint,pmid' ,ii,jj,l & -! ,pint(ii,jj,l),pmid(ii,jj,l) - end do - -!compute pmid from averaged two layer pint - do l=lm,1,-1 - do j=jsta,jend - do i=1,im - pmid(i,j,l) = 0.5*(pint(i,j,l)+pint(i,j,l+1)) - enddo - enddo - enddo - -! surface height from FV3 -! dong set missing value for zint -! zint=spval - VarName='hgtsfc' - call read_netcdf_2d_para(ncid3d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,zint(1,jsta_2l,lp1)) - if(debugprint)print*,'sample ',VarName,' =',zint(isa,jsa,lp1) - do j=jsta,jend - do i=1,im - if (zint(i,j,lp1) /= spval) then - fis(i,j) = zint(i,j,lp1) * grav - else - fis(i,j) = spval - endif - enddo - enddo - - do l=lm,1,-1 - do j=jsta,jend - do i=1,im - if(zint(i,j,l+1)/=spval .and. buf3d(i,j,l)/=spval)then -!make sure delz is positive - zint(i,j,l)=zint(i,j,l+1)+abs(buf3d(i,j,l)) -! if(zint(i,j,l)>1.0E6)print*,'bad H ',i,j,l,zint(i,j,l) - else - zint(i,j,l)=spval - end if - end do - end do - print*,'sample zint= ',isa,jsa,l,zint(isa,jsa,l) - end do - - do l=lp1,1,-1 - do j=jsta,jend - do i=1,im - alpint(i,j,l)=log(pint(i,j,l)) - end do - end do - end do - - do l=lm,1,-1 - do j=jsta,jend - do i=1,im - if(zint(i,j,l+1)/=spval .and. zint(i,j,l)/=spval & - .and. pmid(i,j,l)/=spval)then - zmid(i,j,l)=zint(i,j,l+1)+(zint(i,j,l)-zint(i,j,l+1))* & - (log(pmid(i,j,l))-alpint(i,j,l+1))/ & - (alpint(i,j,l)-alpint(i,j,l+1)) - if(zmid(i,j,l)>1.0E6)print*,'bad Hmid ',i,j,l,zmid(i,j,l) - else - zmid(i,j,l)=spval - endif - end do - end do - end do - - - pt = ak5(1) - -! - - deallocate (vcoord4) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - -! - -! done with 3d file, close it for now - Status=nf90_close(ncid3d) - deallocate(recname) - -! open flux file - Status = nf90_open(trim(fileNameFlux),ior(nf90_nowrite, nf90_mpiio), & - ncid2d,comm=mpi_comm_world,info=mpi_info_null) - if ( Status /= 0 ) then - print*,'error opening ',fileNameFlux, ' Status = ', Status - print*,'skip reading of flux file' - endif - -! IVEGSRC=1 for IGBP, 0 for USGS, 2 for UMD - VarName='IVEGSRC' - Status=nf90_get_att(ncid2d,nf90_global,'IVEGSRC',IVEGSRC) - if (Status /= 0) then - print*,VarName,' not found-Assigned 1 for IGBP as default' - IVEGSRC=1 - end if - if (me == 0) print*,'IVEGSRC= ',IVEGSRC - -! set novegtype based on vegetation classification - if(ivegsrc==2)then - novegtype=13 - else if(ivegsrc==1)then - novegtype=20 - else if(ivegsrc==0)then - novegtype=24 - end if - if (me == 0) print*,'novegtype= ',novegtype - - Status=nf90_get_att(ncid2d,nf90_global,'imp_physics',imp_physics) - if (Status /= 0) then - print*,VarName,' not found-Assigned 11 GFDL as default' - imp_physics=11 - end if - if (me == 0) print*,'MP_PHYSICS= ',imp_physics -! - Status=nf90_get_att(ncid2d,nf90_global,'fhzero',fhzero) - if (Status /= 0) then - print*,VarName,' not found-Assigned 3 hours as default' - fhzero=3 - end if - if (me == 0) print*,'fhzero= ',fhzero -! - Status=nf90_get_att(ncid2d,nf90_global,'dtp',dtp) - if (Status /= 0) then - print*,VarName,' not found-Assigned 90s as default' - dtp=90 - end if - if (me == 0) print*,'dtp= ',dtp -! Initializes constants for Ferrier microphysics - if(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95) then - CALL MICROINIT(imp_physics) - end if - - tprec = float(fhzero) - if(ifhr>240)tprec=12. - tclod = tprec - trdlw = tprec - trdsw = tprec - tsrfc = tprec - tmaxmin = tprec - td3d = tprec - print*,'tprec = ',tprec - - -!Set REF_10CM as missning since gfs doesn't ouput it - do l=1,lm - do j=jsta,jend - do i=1,im - REF_10CM(i,j,l)=spval - enddo - enddo - enddo - - VarName='land' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,sm) - if(debugprint)print*,'sample ',VarName,' =',sm(im/2,(jsta+jend)/2) - -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= spval) sm(i,j) = 1.0 - sm(i,j) - enddo - enddo - -! sea ice mask - - VarName = 'icec' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,sice) - if(debugprint)print*,'sample ',VarName,' = ',sice(isa,jsa) - -! where(sice /=spval .and. sice >=1.0)sm=0.0 !sea ice has sea -! mask=0 -! GFS flux files have land points with non-zero sea ice, per Iredell, -! these -! points have sea ice changed to zero, i.e., trust land mask more than -! sea ice -! where(sm/=spval .and. sm==0.0)sice=0.0 !specify sea ice=0 at land - -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= spval .and. sm(i,j) == 0.0) sice(i,j) = 0.0 - enddo - enddo - - -! PBL height using nemsio - VarName = 'hpbl' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,pblh) - if(debugprint)print*,'sample ',VarName,' = ',pblh(isa,jsa) - -! frictional velocity using nemsio - VarName='fricv' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,ustar) -! if(debugprint)print*,'sample ',VarName,' = ',ustar(isa,jsa) - -! roughness length using getgb - VarName='sfcr' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,z0) -! if(debugprint)print*,'sample ',VarName,' = ',z0(isa,jsa) - -! sfc exchange coeff - VarName='sfexc' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,SFCEXC) - -! aerodynamic conductance - VarName='acond' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,acond) - if(debugprint)print*,'sample ',VarName,' = ',acond(isa,jsa) - -! mid day avg albedo - VarName='albdo_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgalbedo) - if(debugprint)print*,'sample ',VarName,' = ',avgalbedo(isa,jsa) - do j=jsta,jend - do i=1,im - if (avgalbedo(i,j) /= spval) avgalbedo(i,j) = avgalbedo(i,j) * 0.01 - enddo - enddo - -! surface potential T using getgb - VarName='tmpsfc' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,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 - QS(i,j) = SPVAL ! GFS does not have surface specific humidity - twbs(i,j) = SPVAL ! GFS does not have inst sensible heat flux - qwbs(i,j) = SPVAL ! GFS does not have inst latent heat flux -!assign sst - if (sm(i,j) /= 0.0) then - sst(i,j) = ths(i,j) * (pint(i,j,lp1)/p1000)**capa - else - sst(i,j) = spval - endif - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',ths(isa,jsa) - -! foundation temperature - VarName='tref' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,fdnsst) - if(debugprint)print*,'sample ',VarName,' = ',fdnsst(isa,jsa) - -! GFS does not have time step and physics time step, make up ones since they -! are not really used anyway -! NPHS=1. -! DT=90. -! DTQ2 = DT * NPHS !MEB need to get physics DT - DTQ2 = DTP !MEB need to get physics DT - NPHS=1 - DT = DTQ2/NPHS !MEB need to get DT - TSPH = 3600./DT - -! convective precip in m per physics time step using getgb -! read 3 hour bucket - VarName='cpratb_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgcprate) -! where(avgcprate /= spval)avgcprate=avgcprate*dtq2/1000. ! convert to m -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgcprate(i,j) /= spval) avgcprate(i,j) = avgcprate(i,j) * (dtq2*0.001) - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',avgcprate(isa,jsa) - -! print*,'maxval CPRATE: ', maxval(CPRATE) - -! read continuous bucket - VarName='cprat_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgcprate_cont) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgcprate_cont(i,j) /= spval) avgcprate_cont(i,j) = & - avgcprate_cont(i,j) * (dtq2*0.001) - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',avgcprate_cont(isa,jsa) - -! print*,'maxval CPRATE: ', maxval(CPRATE) - -! precip rate in m per physics time step using getgb - VarName='prateb_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgprec) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if(avgprec(i,j) /= spval)avgprec(i,j)=avgprec(i,j)*(dtq2*0.001) - enddo - enddo - - if(debugprint)print*,'sample ',VarName,' = ',avgprec(isa,jsa) - -! prec = avgprec !set avg cprate to inst one to derive other fields - - VarName='prate_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgprec_cont) -! where(avgprec /= spval)avgprec=avgprec*dtq2/1000. ! convert to m -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgprec_cont(i,j) /=spval)avgprec_cont(i,j)=avgprec_cont(i,j) & - * (dtq2*0.001) - enddo - enddo - - if(debugprint)print*,'sample ',VarName,' = ',avgprec_cont(isa,jsa) -! precip rate in m per physics time step - VarName='tprcp' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,prec) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (prec(i,j) /= spval) prec(i,j)=prec(i,j)* (dtq2*0.001) & - * 1000. / dtp - enddo - enddo - -! convective precip rate in m per physics time step - VarName='cnvprcp' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,cprate) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (cprate(i,j) /= spval) then - cprate(i,j) = max(0.,cprate(i,j)) * (dtq2*0.001) & - * 1000. / dtp - else - cprate(i,j) = 0. - endif - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',cprate(isa,jsa) - -! GFS does not have accumulated total, gridscale, and convective precip, will use inst precip to derive in SURFCE.f - -! max hourly 1-km agl reflectivity -! VarName='refdmax' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,refd_max) -! if(debugprint)print*,'sample ',VarName,' = ',refd_max(isa,jsa) -! max hourly -10C reflectivity -! VarName='refdmax263k' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,refdm10c_max) -! if(debugprint)print*,'sample ',VarName,' = ',refdm10c_max(isa,jsa) - -! max hourly u comp of 10m agl wind -! VarName='u10max' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,u10max) -! if(debugprint)print*,'sample ',VarName,' = ',u10max(isa,jsa) -! max hourly v comp of 10m agl wind -! VarName='v10max' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,v10max) -! if(debugprint)print*,'sample ',VarName,' = ',v10max(isa,jsa) -! max hourly 10m agl wind speed -! VarName='spd10max' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,wspd10max) -! if(debugprint)print*,'sample ',VarName,' = ',wspd10max(isa,jsa) - - -! 2m T using nemsio - VarName='tmp2m' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,tshltr) - if(debugprint)print*,'sample ',VarName,' = ',tshltr(isa,jsa) - -! inst snow water eqivalent using nemsio - VarName='weasd' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,sno) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j)==0.) sno(i,j) = spval - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',sno(isa,jsa) - -! ave snow cover - VarName='snowc_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,snoavg) -! snow cover is multipled by 100 in SURFCE before writing it out - do j=jsta,jend - do i=1,im - if (sm(i,j)==1.0 .and. sice(i,j)==0.) snoavg(i,j)=spval - if(snoavg(i,j)/=spval)snoavg(i,j)=snoavg(i,j)/100. - end do - end do - -! snow depth in mm using nemsio - VarName='snod' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,si) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j)==1.0 .and. sice(i,j)==0.) si(i,j)=spval - if (si(i,j) /= spval) si(i,j) = si(i,j) * 1000.0 - CLDEFI(i,j) = SPVAL ! GFS does not have convective cloud efficiency - lspa(i,j) = spval ! GFS does not have similated precip - TH10(i,j) = SPVAL ! GFS does not have 10 m theta - TH10(i,j) = SPVAL ! GFS does not have 10 m theta - Q10(i,j) = SPVAL ! GFS does not have 10 m humidity - ALBASE(i,j) = SPVAL ! GFS does not have snow free albedo - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',si(isa,jsa) - -! 2m T using nemsio - VarName='tmp2m' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,tshltr) - if(debugprint)print*,'sample ',VarName,' = ',tshltr(isa,jsa) - -! GFS does not have 2m pres, estimate it, also convert t to theta - Do j=jsta,jend - Do i=1,im - 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 -! if (j == jm/2 .and. mod(i,50) == 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 nemsio - VarName='spfh2m' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,qshltr) - if(debugprint)print*,'sample ',VarName,' = ',qshltr(isa,jsa) - -! mid day avg albedo in fraction using nemsio -! VarName='albdosfc' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,avgalbedo) -!! where(avgalbedo /= spval)avgalbedo=avgalbedo/100. ! convert to fraction -!!$omp parallel do private(i,j) -! do j=jsta,jend -! do i=1,im -! if (avgalbedo(i,j) /= spval) avgalbedo(i,j) = avgalbedo(i,j) * 0.01 -! enddo -! enddo -! if(debugprint)print*,'sample ',VarName,' = ',avgalbedo(isa,jsa) - -! time averaged column cloud fractionusing nemsio - VarName='tcdc_aveclm' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgtcdc) -! where(avgtcdc /= spval)avgtcdc=avgtcdc/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgtcdc(i,j) /= spval) avgtcdc(i,j) = avgtcdc(i,j) * 0.01 - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',avgtcdc(isa,jsa) - -! GFS probably does not use zenith angle -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - Czen(i,j) = spval - CZMEAN(i,j) = SPVAL - enddo - enddo - -! maximum snow albedo in fraction using nemsio - VarName='snoalb' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,mxsnal) -! where(mxsnal /= spval)mxsnal=mxsnal/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (mxsnal(i,j) /= spval) mxsnal(i,j) = mxsnal(i,j) * 0.01 - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',mxsnal(isa,jsa) - -! land fraction - VarName='lfrac' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,landfrac) - -! 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) * T(I,J,LM) - Sigt4(i,j) = 5.67E-8 * TLMH * TLMH - End do - End do - -! TG is not used, skip it for now - -! GFS does not have inst cloud fraction for high, middle, and low cloud -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - cfrach(i,j) = spval - cfracl(i,j) = spval - cfracm(i,j) = spval - enddo - enddo - -! ave high cloud fraction using nemsio - VarName='tcdc_avehcl' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgcfrach) -! where(avgcfrach /= spval)avgcfrach=avgcfrach/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgcfrach(i,j) /= spval) avgcfrach(i,j) = avgcfrach(i,j) * 0.01 - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',avgcfrach(isa,jsa) - -! ave low cloud fraction using nemsio - VarName='tcdc_avelcl' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgcfracl) -! where(avgcfracl /= spval)avgcfracl=avgcfracl/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgcfracl(i,j) /= spval) avgcfracl(i,j) = avgcfracl(i,j) * 0.01 - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',avgcfracl(isa,jsa) - -! ave middle cloud fraction using nemsio - VarName='tcdc_avemcl' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgcfracm) -! where(avgcfracm /= spval)avgcfracm=avgcfracm/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgcfracm(i,j) /= spval) avgcfracm(i,j) = avgcfracm(i,j) * 0.01 - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',avgcfracm(isa,jsa) - -! inst convective cloud fraction using nemsio - VarName='tcdccnvcl' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,cnvcfr) -! where(cnvcfr /= spval)cnvcfr=cnvcfr/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (cnvcfr(i,j) /= spval) cnvcfr (i,j)= cnvcfr(i,j) * 0.01 - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',cnvcfr(isa,jsa) - -! slope type using nemsio - VarName='sltyp' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,buf) -!$omp parallel do private(i,j) - do j = jsta_2l, jend_2u - do i=1,im - if (buf(i,j) < spval) then - islope(i,j) = nint(buf(i,j)) - else - islope(i,j) = 0 - endif - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',islope(isa,jsa) - -! plant canopy sfc wtr in m - VarName='cnwat' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,cmc) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (cmc(i,j) /= spval) cmc(i,j) = cmc(i,j) * 0.001 - if (sm(i,j) /= 0.0) cmc(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',cmc(isa,jsa) - -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - grnflx(i,j) = spval ! GFS does not have inst ground heat flux - enddo - enddo - -! frozen precip fraction using nemsio - VarName='cpofp' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,sr) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if(sr(i,j) /= spval) then -!set range within (0,1) - sr(i,j)=min(1.,max(0.,sr(i,j))) - endif - enddo - enddo - -! sea ice skin temperature - VarName='tisfc' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,ti) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sice(i,j) == spval .or. sice(i,j) == 0.) ti(i,j)=spval - enddo - enddo - -! vegetation fraction in fraction. using nemsio - VarName='veg' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,vegfrc) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (vegfrc(i,j) /= spval) then - vegfrc(i,j) = vegfrc(i,j) * 0.01 - else - vegfrc(i,j) = 0.0 - endif - enddo - enddo -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) vegfrc(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample ',VarName,' = ',vegfrc(isa,jsa) - -! 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='soill1' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,sh2o(1,jsta_2l,1)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) sh2o(i,j,1) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,sh2o(isa,jsa,1) - - VarName='soill2' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,sh2o(1,jsta_2l,2)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) sh2o(i,j,2) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,sh2o(isa,jsa,2) - - VarName='soill3' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,sh2o(1,jsta_2l,3)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) sh2o(i,j,3) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,sh2o(isa,jsa,3) - - VarName='soill4' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,sh2o(1,jsta_2l,4)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) sh2o(i,j,4) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,sh2o(isa,jsa,4) - -! volumetric soil moisture using nemsio - VarName='soilw1' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,smc(1,jsta_2l,1)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) smc(i,j,1) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,1) - - VarName='soilw2' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,smc(1,jsta_2l,2)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) smc(i,j,2) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,2) - - VarName='soilw3' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,smc(1,jsta_2l,3)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) smc(i,j,3) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,3) - - VarName='soilw4' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,smc(1,jsta_2l,4)) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) smc(i,j,4) = spval - enddo - enddo - if(debugprint)print*,'sample l',VarName,' = ',1,smc(isa,jsa,4) - -! soil temperature using nemsio - VarName='soilt1' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,stc(1,jsta_2l,1)) -! mask open water areas, combine with sea ice tmp -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,1) = spval - !if (sm(i,j) /= 0.0) stc(i,j,1) = spval - enddo - enddo - if(debugprint)print*,'sample l','stc',' = ',1,stc(isa,jsa,1) - - VarName='soilt2' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,stc(1,jsta_2l,2)) -! mask open water areas, combine with sea ice tmp -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,2) = spval - !if (sm(i,j) /= 0.0) stc(i,j,2) = spval - enddo - enddo - if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,2) - - VarName='soilt3' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,stc(1,jsta_2l,3)) -! mask open water areas, combine with sea ice tmp -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,3) = spval - !if (sm(i,j) /= 0.0) stc(i,j,3) = spval - enddo - enddo - if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,3) - - VarName='soilt4' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,stc(1,jsta_2l,4)) -! mask open water areas, combine with sea ice tmp -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) stc(i,j,4) = spval - !if (sm(i,j) /= 0.0) stc(i,j,4) = spval - enddo - enddo - if(debugprint)print*,'sample stc = ',1,stc(isa,jsa,4) - -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - acfrcv(i,j) = spval ! GFS does not output time averaged convective and strat cloud fraction, set acfrcv to spval, ncfrcv to 1 - ncfrcv(i,j) = 1.0 - acfrst(i,j) = spval ! GFS does not output time averaged cloud fraction, set acfrst to spval, ncfrst to 1 - ncfrst(i,j) = 1.0 - bgroff(i,j) = spval ! GFS does not have UNDERGROUND RUNOFF - rlwtoa(i,j) = spval ! GFS does not have inst model top outgoing longwave - enddo - enddo -! trdlw(i,j) = 6.0 - ardlw = 1.0 ! GFS incoming sfc longwave has been averaged over 6 hr bucket, set ARDLW to 1 - -! time averaged incoming sfc longwave - VarName='dlwrf_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,alwin) - -! inst incoming sfc longwave - VarName='dlwrf' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,rlwin) - -! time averaged outgoing sfc longwave - VarName='ulwrf_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,alwout) -! inst outgoing sfc longwave - VarName='ulwrf' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,radot) - -! where(alwout /= spval) alwout=-alwout ! CLDRAD puts a minus sign before gribbing -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (alwout(i,j) /= spval) alwout(i,j) = -alwout(i,j) - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,alwout(isa,jsa) - -! time averaged outgoing model top longwave using gfsio - VarName='ulwrf_avetoa' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,alwtoa) -! if(debugprint)print*,'sample l',VarName,' = ',1,alwtoa(isa,jsa) - -! GFS incoming sfc longwave has been averaged, set ARDLW to 1 - ardsw=1.0 -! trdsw=6.0 - -! time averaged incoming sfc shortwave - VarName='dswrf_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,aswin) -! if(debugprint)print*,'sample l',VarName,' = ',1,aswin(isa,jsa) - -! inst incoming sfc shortwave - VarName='dswrf' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,rswin) - -! inst incoming clear sky sfc shortwave - VarName='csdlf' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,rswinc) - -! time averaged incoming sfc uv-b using getgb - VarName='duvb_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,auvbin) -! if(debugprint)print*,'sample l',VarName,' = ',1,auvbin(isa,jsa) - -! time averaged incoming sfc clear sky uv-b using getgb - VarName='cduvb_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,auvbinc) -! if(debugprint)print*,'sample l',VarName,' = ',1,auvbinc(isa,jsa) - -! time averaged outgoing sfc shortwave using gfsio - VarName='uswrf_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,aswout) -! where(aswout /= spval) aswout=-aswout ! CLDRAD puts a minus sign before gribbing -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (aswout(i,j) /= spval) aswout(i,j) = -aswout(i,j) - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,aswout(isa,jsa) - -! inst outgoing sfc shortwave using gfsio - VarName='uswrf' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,rswout) - -! time averaged model top incoming shortwave - VarName='dswrf_avetoa' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,aswintoa) -! if(debugprint)print*,'sample l',VarName,' = ',1,aswintoa(isa,jsa) - -! time averaged model top outgoing shortwave - VarName='uswrf_avetoa' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,aswtoa) -! if(debugprint)print*,'sample l',VarName,' = ',1,aswtoa(isa,jsa) - -! time averaged surface sensible heat flux, multiplied by -1 because wrf model flux -! has reversed sign convention using gfsio - VarName='shtfl_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,sfcshx) -! where (sfcshx /= spval)sfcshx=-sfcshx -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sfcshx(i,j) /= spval) sfcshx(i,j) = -sfcshx(i,j) - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,sfcshx(isa,jsa) - -! inst surface sensible heat flux - VarName='shtfl' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,twbs) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (twbs(i,j) /= spval) twbs(i,j) = -twbs(i,j) - enddo - enddo - -! 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 - VarName='lhtfl_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,sfclhx) -! where (sfclhx /= spval)sfclhx=-sfclhx -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sfclhx(i,j) /= spval) sfclhx(i,j) = -sfclhx(i,j) - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,sfclhx(isa,jsa) - -! inst surface latent heat flux - VarName='lhtfl' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,qwbs) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (qwbs(i,j) /= spval) qwbs(i,j) = -qwbs(i,j) - enddo - enddo - - if(me==0)print*,'rdaod= ',rdaod -! inst aod550 optical depth - if(rdaod) then - VarName='aod550' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,aod550) - - VarName='du_aod550' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,du_aod550) - - VarName='ss_aod550' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,ss_aod550) - - VarName='su_aod550' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,su_aod550) - - VarName='oc_aod550' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,oc_aod550) - - VarName='bc_aod550' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,bc_aod550) - endif !end if rdaod - - -! time averaged ground heat flux using nemsio - VarName='gflux_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,subshx) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) subshx(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,subshx(isa,jsa) - -! inst ground heat flux using nemsio - VarName='gflux' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,grnflx) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) grnflx(i,j) = spval - enddo - enddo - -! time averaged zonal momentum flux using gfsio - VarName='uflx_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,sfcux) -! if(debugprint)print*,'sample l',VarName,' = ',1,sfcux(isa,jsa) - -! time averaged meridional momentum flux using nemsio - VarName='vflx_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,sfcvx) -! if(debugprint)print*,'sample l',VarName,' = ',1,sfcvx(isa,jsa) - -! dong read in inst surface flux -! inst zonal momentum flux using gfsio -! VarName='uflx' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sfcuxi) -! if(debugprint)print*,'sample l',VarName,' = ',1,sfcuxi(isa,jsa) - -! inst meridional momentum flux using nemsio -! VarName='vflx' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,sfcvxi) -! if(debugprint)print*,'sample l',VarName,' = ',1,sfcvxi(isa,jsa) - - -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - sfcuvx(i,j) = spval ! GFS does not use total momentum flux - enddo - enddo - -! time averaged zonal gravity wave stress using nemsio - VarName='u-gwd_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,gtaux) -! if(debugprint)print*,'sample l',VarName,' = ',1,gtaux(isa,jsa) - -! time averaged meridional gravity wave stress using getgb - VarName='v-gwd_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,gtauy) -! if(debugprint)print*,'sample l',VarName,' = ',1,gtauy(isa,jsa) - -! time averaged accumulated potential evaporation - VarName='pevpr_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgpotevp) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) avgpotevp(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,potevp(isa,jsa) - -! inst potential evaporation - VarName='pevpr' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,potevp) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) == 1.0 .and. sice(i,j) ==0.) potevp(i,j) = spval - enddo - enddo - - do l=1,lm -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im -! GFS does not have temperature tendency due to long wave radiation - rlwtt(i,j,l) = spval -! GFS does not have temperature tendency due to short wave radiation - rswtt(i,j,l) = spval -! GFS does not have temperature tendency due to latent heating from convection - tcucn(i,j,l) = spval - tcucns(i,j,l) = spval -! GFS does not have temperature tendency due to latent heating from grid scale - train(i,j,l) = spval - enddo - enddo - enddo - -! set avrain to 1 - avrain=1.0 - avcnvc=1.0 - theat=6.0 ! just in case GFS decides to output T tendency - -! 10 m u using nemsio - VarName='ugrd10m' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,u10) - - do j=jsta,jend - do i=1,im - u10h(i,j)=u10(i,j) - end do - end do -! if(debugprint)print*,'sample l',VarName,' = ',1,u10(isa,jsa) - -! 10 m v using gfsio - VarName='vgrd10m' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,v10) - - do j=jsta,jend - do i=1,im - v10h(i,j)=v10(i,j) - end do - end do -! if(debugprint)print*,'sample l',VarName,' = ',1,v10(isa,jsa) - -! vegetation type, it's in GFS surface file, hopefully will merge into gfsio soon - VarName='vtype' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,buf) -! where (buf /= spval) -! ivgtyp=nint(buf) -! elsewhere -! ivgtyp=0 !need to feed reasonable value to crtm -! end where -!$omp parallel do private(i,j) - do j = jsta_2l, jend_2u - do i=1,im - if (buf(i,j) < spval) then - ivgtyp(i,j) = nint(buf(i,j)) - else - ivgtyp(i,j) = 0 - endif - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,ivgtyp(isa,jsa) - -! soil type, it's in GFS surface file, hopefully will merge into gfsio soon - VarName='sotyp' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,buf) -!$omp parallel do private(i,j) - do j = jsta_2l, jend_2u - do i=1,im - if (buf(i,j) < spval) then - isltyp(i,j) = nint(buf(i,j)) - else - isltyp(i,j) = 0 !need to feed reasonable value to crtm - endif - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,isltyp(isa,jsa) - -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - smstav(i,j) = spval ! GFS does not have soil moisture availability -! smstot(i,j) = spval ! GFS does not have total soil moisture - sfcevp(i,j) = spval ! GFS does not have accumulated surface evaporation - acsnow(i,j) = spval ! GFS does not have averaged accumulated snow - acsnom(i,j) = spval ! GFS does not have snow melt -! sst(i,j) = spval ! GFS does not have sst???? - thz0(i,j) = ths(i,j) ! GFS does not have THZ0, use THS to substitute - qz0(i,j) = spval ! GFS does not output humidity at roughness length - uz0(i,j) = spval ! GFS does not output u at roughness length - vz0(i,j) = spval ! GFS does not output humidity at roughness length - enddo - enddo - do l=1,lm -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - EL_PBL(i,j,l) = spval ! GFS does not have mixing length - exch_h(i,j,l) = spval ! GFS does not output exchange coefficient - enddo - enddo - enddo -! if(debugprint)print*,'sample l',VarName,' = ',1,thz0(isa,jsa) - -! 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' -! l=1 -! if(debugprint)print*,'sample l',VarName,' = ',1,ptop(isa,jsa) - VarName='prescnvclt' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,ptop) - - -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - htop(i,j) = spval - if(ptop(i,j) <= 0.0) ptop(i,j) = spval - enddo - enddo - do j=jsta,jend - do i=1,im - if(ptop(i,j) < spval)then - do l=1,lm - if(ptop(i,j) <= pmid(i,j,l))then - htop(i,j) = l -! if(i==ii .and. j==jj)print*,'sample ptop,pmid pmid-1,pint= ', & -! ptop(i,j),pmid(i,j,l),pmid(i,j,l-1),pint(i,j,l),htop(i,j) - exit - end if - end do - end if - end do - end do - -! 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='prescnvclb' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,pbot) -! if(debugprint)print*,'sample l',VarName,VcoordName,' = ',1,pbot(isa,jsa) -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - hbot(i,j) = spval - if(pbot(i,j) <= 0.0) pbot(i,j) = spval - enddo - enddo - do j=jsta,jend - do i=1,im -! if(.not.lb(i,j))print*,'false bitmask for pbot at ' -! + ,i,j,pbot(i,j) - if(pbot(i,j) < spval)then - do l=lm,1,-1 - if(pbot(i,j) >= pmid(i,j,l)) then - hbot(i,j) = l -! if(i==ii .and. j==jj)print*,'sample pbot,pmid= ', & -! pbot(i,j),pmid(i,j,l),hbot(i,j) - exit - end if - end do - end if - end do - end do - if(debugprint)print*,'sample hbot = ',hbot(isa,jsa) -! retrieve time averaged low cloud top pressure using nemsio - VarName='pres_avelct' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,ptopl) -! if(debugprint)print*,'sample l',VarName,' = ',1,ptopl(isa,jsa) - -! retrieve time averaged low cloud bottom pressure using nemsio - VarName='pres_avelcb' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,pbotl) -! if(debugprint)print*,'sample l',VarName,' = ',1,pbotl(isa,jsa) - -! retrieve time averaged low cloud top temperature using nemsio - VarName='tmp_avelct' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,Ttopl) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopl(isa,jsa) - -! retrieve time averaged middle cloud top pressure using nemsio - VarName='pres_avemct' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,ptopm) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptopm(isa,jsa) - -! retrieve time averaged middle cloud bottom pressure using nemsio - VarName='pres_avemcb' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,pbotm) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pbotm(isa,jsa) - -! retrieve time averaged middle cloud top temperature using nemsio - VarName='tmp_avemct' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,Ttopm) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopm(isa,jsa) - -! retrieve time averaged high cloud top pressure using nemsio ********* - VarName='pres_avehct' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,ptoph) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptoph(isa,jsa) - -! retrieve time averaged high cloud bottom pressure using nemsio - VarName='pres_avehcb' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,pboth) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pboth(isa,jsa) - -! retrieve time averaged high cloud top temperature using nemsio - VarName='tmp_avehct' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,Ttoph) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttoph(isa,jsa) - -! retrieve boundary layer cloud cover using nemsio - VarName='tcdc_avebndcl' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,pblcfr) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,pblcfr(isa,jsa) -! where (pblcfr /= spval)pblcfr=pblcfr/100. ! convert to fraction -!$omp parallel do private(i,j) - do j = jsta_2l, jend_2u - do i=1,im - if (pblcfr(i,j) < spval) pblcfr(i,j) = pblcfr(i,j) * 0.01 - enddo - enddo - -! retrieve cloud work function - VarName='cwork_aveclm' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,cldwork) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,cldwork(isa,jsa) - -! accumulated total (base+surface) runoff - VarName='watr_acc' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,runoff) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) runoff(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,runoff(isa,jsa) - -! accumulated evaporation of intercepted water - VarName='ecan_acc' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,tecan) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) tecan(i,j) = spval - enddo - enddo - -! accumulated plant transpiration - VarName='etran_acc' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,tetran) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) tetran(i,j) = spval - enddo - enddo - -! accumulated soil surface evaporation - VarName='edir_acc' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,tedir) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) tedir(i,j) = spval - enddo - enddo - -! total water storage in aquifer - VarName='wa_acc' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,twa) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) twa(i,j) = spval - enddo - enddo - -! retrieve shelter max temperature using nemsio - VarName='tmax_max2m' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,maxtshltr) - -! retrieve shelter min temperature using nemsio - VarName='tmin_min2m' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,mintshltr) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', & -! 1,mintshltr(im/2,(jsta+jend)/2) - -! retrieve shelter max RH -! VarName='rh02max' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,maxrhshltr) - -! retrieve shelter min temperature using nemsio -! VarName='rh02min' -! call read_netcdf_2d_scatter(me,ncid2d,1,im,jm,jsta,jsta_2l & -! ,jend_2u,MPI_COMM_COMP,icnt,idsp,spval,VarName,minrhshltr) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', & -! 1,mintshltr(im/2,(jsta+jend)/2) - -!$omp parallel do private(i,j) - do j=jsta_2l,jend_2u - do i=1,im - MAXRHSHLTR(i,j) = SPVAL - MINRHSHLTR(i,j) = SPVAL - enddo - enddo - -! retrieve ice thickness using nemsio - VarName='icetk' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,dzice) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,dzice(isa,jsa) - -! retrieve wilting point using nemsio - VarName='wilt' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,smcwlt) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) smcwlt(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,smcwlt(isa,jsa) - -! retrieve sunshine duration using nemsio - VarName='sunsd_acc' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,suntime) - -! retrieve field capacity using nemsio - VarName='fldcp' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,fieldcapa) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) fieldcapa(i,j) = spval - enddo - enddo -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,fieldcapa(isa,jsa) - -! retrieve time averaged surface visible beam downward solar flux - VarName='vbdsf_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avisbeamswin) - l=1 - -! retrieve time averaged surface visible diffuse downward solar flux - VarName='vddsf_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avisdiffswin) - -! retrieve time averaged surface near IR beam downward solar flux - VarName='nbdsf_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,airbeamswin) - -! retrieve time averaged surface near IR diffuse downward solar flux - VarName='nddsf_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,airdiffswin) - -! retrieve time averaged surface clear sky outgoing LW - VarName='csulf' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,alwoutc) - -! retrieve time averaged TOA clear sky outgoing LW - VarName='csulftoa' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,alwtoac) - -! retrieve time averaged surface clear sky outgoing SW - VarName='csusf' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,aswoutc) - -! retrieve time averaged TOA clear sky outgoing LW - VarName='csusftoa' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,aswtoac) - -! retrieve time averaged surface clear sky incoming LW - VarName='csdlf' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,alwinc) - -! retrieve time averaged surface clear sky incoming SW - VarName='csdsf' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,aswinc) - -! retrieve shelter max specific humidity using nemsio - VarName='spfhmax_max2m' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,maxqshltr) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', -! 1,maxqshltr(isa,jsa) - -! retrieve shelter min temperature using nemsio - VarName='spfhmin_min2m' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,minqshltr) - -! retrieve storm runoff using nemsio - VarName='ssrun_acc' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,SSROFF) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) ssroff(i,j) = spval - enddo - enddo - -! retrieve direct soil evaporation - VarName='evbs_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgedir) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) avgedir(i,j) = spval - enddo - enddo - -! retrieve CANOPY WATER EVAP - VarName='evcw_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgecan) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) avgecan(i,j) = spval - enddo - enddo - -! retrieve AVERAGED PRECIP ADVECTED HEAT FLUX - VarName='pah_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,paha) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) paha(i,j) = spval - enddo - enddo - -! retrieve nstantaneous PRECIP ADVECTED HEAT FLUX - VarName='pahi' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,pahi) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) pahi(i,j) = spval - enddo - enddo - -! retrieve PLANT TRANSPIRATION - VarName='trans_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgetrans) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) avgetrans(i,j) = spval - enddo - enddo - -! retrieve snow sublimation - VarName='sbsno_ave' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgesnow) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j)==1.0 .and. sice(i,j)==0.) avgesnow(i,j)=spval - enddo - enddo - -! retrive total soil moisture - VarName='soilm' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,smstot) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) smstot(i,j) = spval - enddo - enddo - -! retrieve snow phase change heat flux - VarName='snohf' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,snopcx) -! mask water areas -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (sm(i,j) /= 0.0) snopcx(i,j) = spval - enddo - enddo - -! GFS does not have deep convective cloud top and bottom fields - -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - HTOPD(i,j) = SPVAL - HBOTD(i,j) = SPVAL - HTOPS(i,j) = SPVAL - HBOTS(i,j) = SPVAL - CUPPT(i,j) = SPVAL - enddo - enddo - -! done with flux file, close it for now - Status=nf90_close(ncid2d) -! deallocate(tmp,recname,reclevtyp,reclev) - -! 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 == 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 -! - -! generate look up table for lifted parcel calculations - - THL = 210. - PLQ = 70000. - pt_TBL = 10000. ! this is for 100 hPa added by Moorthi - - CALL TABLE(PTBL,TTBL,PT_TBL, & - RDQ,RDTH,RDP,RDTHE,PL,THL,QS0,SQS,STHE,THE0) - - CALL TABLEQ(TTBLQ,RDPQ,RDTHEQ,PLQ,THL,STHEQ,THE0Q) - -! -! - IF(ME == 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 -! -!$omp parallel do private(l) - DO L = 1,LSM - ALSL(L) = LOG(SPL(L)) - END DO -! -!HC WRITE IGDS OUT FOR WEIGHTMAKER TO READ IN AS KGDSIN - if(me == 0)then - print*,'writing out igds' - igdout = 110 -! open(igdout,file='griddef.out',form='unformatted' -! + ,status='unknown') - if(maptype == 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 == 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 - ! Note: The calculation of the map scale factor at the standard - ! lat/lon and the PSMAPF - ! Get map factor at 60 degrees (N or S) for PS projection, which will - ! be needed to correctly define the DX and DY values in the GRIB GDS - if (TRUELAT1 < 0.) THEN - LAT = -60. - else - LAT = 60. - end if - - CALL MSFPS (LAT,TRUELAT1*0.001,PSMAPF) - - ELSE IF(MAPTYPE == 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 == 0 .OR. MAPTYPE == 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 -! -! - - RETURN - END - - subroutine read_netcdf_3d_para(ncid,im,jm,jsta,jsta_2l,jend,jend_2u, & - spval,varname,buf,lm) - - use netcdf - use ctlblk_mod, only : me - use params_mod, only : small - implicit none - INCLUDE "mpif.h" - - character(len=20),intent(in) :: varname - real,intent(in) :: spval - integer,intent(in) :: ncid,im,jm,lm,jsta_2l,jend_2u,jsta,jend - real,intent(out) :: buf(im,jsta_2l:jend_2u,lm) - integer :: varid,iret,jj,i,j,l,kk - integer :: start(3), count(3), stride(3) - real,parameter :: spval_netcdf=9.99e+20 - real :: fill_value - - iret = nf90_inq_varid(ncid,trim(varname),varid) - if (iret /= 0) then - if (me == 0) print*,VarName," not found -Assigned missing values" -!$omp parallel do private(i,j,l) - do l=1,lm - do j=jsta,jend - do i=1,im - buf(i,j,l)=spval - enddo - enddo - enddo - else - iret = nf90_get_att(ncid,varid,"_FillValue",fill_value) - if (iret /= 0) fill_value = spval_netcdf - start = (/1,jsta,1/) - jj=jend-jsta+1 - count = (/im,jj,lm/) - iret = nf90_get_var(ncid,varid,buf(1:im,jsta:jend,1:lm),start=start,count=count) - do l=1,lm - do j=jsta,jend - do i=1,im - if(abs(buf(i,j,l)-fill_value) con_g, fv => con_fvirt, rgas => con_rd, & @@ -80,7 +84,7 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) ardsw, asrfc, avrain, avcnvc, theat, gdsdegr, spl, lsm, alsl, im, jm, im_jm, lm, & jsta_2l, jend_2u, nsoil, lp1, icu_physics, ivegsrc, novegtype, nbin_ss, nbin_bc, & nbin_oc, nbin_su, gocart_on, pt_tbl, hyb_sigp, filenameFlux, fileNameAER, & - iSF_SURFACE_PHYSICS,rdaod, aqfcmaq_on + iSF_SURFACE_PHYSICS,rdaod, aqfcmaq_on, modelname use gridspec_mod, only: maptype, gridtype, latstart, latlast, lonstart, lonlast, cenlon, & dxval, dyval, truelat2, truelat1, psmapf, cenlat,lonstartv, lonlastv, cenlonv, & latstartv, latlastv, cenlatv,latstart_r,latlast_r,lonstart_r,lonlast_r, STANDLON @@ -160,7 +164,6 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) real, allocatable :: wrk1(:,:), wrk2(:,:) real, allocatable :: p2d(:,:), t2d(:,:), q2d(:,:), & qs2d(:,:), cw2d(:,:), cfr2d(:,:) - real(kind=4),allocatable :: vcoord4(:,:,:) real, dimension(lm+1) :: ak5, bk5 real*8, allocatable :: pm2d(:,:), pi2d(:,:) real, allocatable :: tmp(:) @@ -382,7 +385,7 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) !$omp parallel do private(i,j) do j = jsta_2l, jend_2u - do i=1,im + do i= 1, im buf(i,j) = spval enddo enddo @@ -617,9 +620,13 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) end if STANDLON = cenlon - print*,'lonstart,latstart,cenlon,cenlat,truelat1,truelat2,stadlon,dyval,dxval', & + print*,'lonstart,latstart,cenlon,cenlat,truelat1,truelat2, & + stadlon,dyval,dxval', & lonstart,latstart,cenlon,cenlat,truelat1,truelat2,standlon,dyval,dxval + else if(trim(varcharval)=='gaussian')then + MAPTYPE=4 + idrt=4 else ! setting default maptype MAPTYPE=0 idrt=0 @@ -677,7 +684,6 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) ! write(0,*)'nrec=',nrec !allocate(recname(nrec),reclevtyp(nrec),reclev(nrec)) allocate(glat1d(jm),glon1d(im)) - allocate(vcoord4(lm+1,3,2)) ! hardwire idate for now ! idate=(/2017,08,07,00,0,0,0,0/) @@ -712,7 +718,7 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) ! Jili Dong check output format for coordinate reading Status=nf90_inq_varid(ncid3d,'grid_xt',varid) Status=nf90_inquire_variable(ncid3d,varid,ndims = numDims) - if(numDims==1) then + if(numDims==1.and.modelname=="FV3R") then read_lonlat=.true. else read_lonlat=.false. @@ -836,7 +842,7 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) if(debugprint)print*,'me sample gdlon gdlat= ' & ,me,gdlon(isa,jsa),gdlat(isa,jsa) -! Specigy grid staggering type +! Specify grid staggering type gridtype = 'A' maptype=idrt if (me == 0) print *, 'maptype and gridtype is ', & @@ -859,21 +865,11 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) end if end if -! -! if (hyb_sigp) then -! do l=1,lm+1 -! ak5(l) = vcoord4(l,1,1) -! bk5(l) = vcoord4(l,2,1) -! enddo -! endif - -!--Fanglin Yang: nemsio file created from FV3 does not have vcoord. -! if ( minval(ak5) <0 .or. minval(bk5) <0 ) then deallocate(glat1d,glon1d) print*,'idate = ',(idate(i),i=1,7) - print*,'nfhour = ',nfhour +! print*,'nfhour = ',nfhour ! sample print point ii = im/2 @@ -1013,8 +1009,6 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) spval,recname(13),qqs(1,jsta_2l,1),lm) call read_netcdf_3d_para(ncid3d,im,jm,jsta,jsta_2l,jend,jend_2u, & spval,recname(14),qqg(1,jsta_2l,1),lm) -!wm call read_netcdf_3d_para(ncid3d,im,jm,jsta,jsta_2l,jend,jend_2u, & -!wm spval,recname(15),cfr(1,jsta_2l,1),lm) ! calculate CWM from FV3 output do l=1,lm @@ -1030,6 +1024,21 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) cwm(isa,jsa,l) end do +! instantaneous 3D cloud fraction + if ( imp_physics==11) then !GFDL MP + VarName='cld_amt' + call read_netcdf_3d_para(ncid3d,im,jm,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,cfr(1,jsta_2l,1),lm) + else + VarName='cldfra' + call read_netcdf_3d_para(ncid2d,im,jm,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,cfr(1,jsta_2l,1),lm) + endif +! do l=1,lm +! if(debugprint)print*,'sample ',VarName,'isa,jsa,l =' & +! ,cfr(isa,jsa,l),isa,jsa,l +! enddo + !============================= ! For AQF Chemical species !============================= @@ -1515,12 +1524,13 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) endif ! -- aqfcmaq_on !============================ +! read for regional FV3 + if (modelname == 'FV3R') then ! max hourly updraft velocity VarName='upvvelmax' call read_netcdf_2d_para(ncid3d,im,jsta,jsta_2l,jend,jend_2u, & spval,VarName,w_up_max) if(debugprint)print*,'sample ',VarName,' = ',w_up_max(isa,jsa) - ! max hourly downdraft velocity VarName='dnvvelmax' call read_netcdf_2d_para(ncid3d,im,jsta,jsta_2l,jend,jend_2u, & @@ -1562,6 +1572,8 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) call read_netcdf_2d_para(ncid3d,im,jsta,jsta_2l,jend,jend_2u, & spval,VarName,rel_vort_maxhy1) if(debugprint)print*,'sample ',VarName,' =',rel_vort_maxhy1(isa,jsa) + endif + ! surface pressure VarName='pressfc' call read_netcdf_2d_para(ncid3d,im,jsta,jsta_2l,jend,jend_2u, & @@ -1666,10 +1678,6 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) end do - pt = ak5(1) - - - deallocate (vcoord4) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! @@ -1697,23 +1705,16 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) end if if (me == 0) print*,'novegtype= ',novegtype - Status=nf90_get_att(ncid2d,nf90_global,'imp_physics',imp_physics) - if (Status /= 0) then - print*,VarName,' not found-Assigned 11 GFDL as default' - imp_physics=11 - end if - if (me == 0) print*,'MP_PHYSICS= ',imp_physics -! Status=nf90_get_att(ncid2d,nf90_global,'fhzero',fhzero) if (Status /= 0) then - print*,VarName,' not found-Assigned 3 hours as default' + print*,'fhzero not found-Assigned 3 hours as default' fhzero=3 end if if (me == 0) print*,'fhzero= ',fhzero ! Status=nf90_get_att(ncid2d,nf90_global,'dtp',dtp) if (Status /= 0) then - print*,VarName,' not found-Assigned 90s as default' + print*,'dtp not found-Assigned 90s as default' dtp=90 end if if (me == 0) print*,'dtp= ',dtp @@ -1733,15 +1734,6 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) print*,'tprec = ',tprec -! instantaneous 3D cloud fraction - VarName='cldfra' - call read_netcdf_3d_para(ncid2d,im,jm,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,cfr(1,jsta_2l,1),lm) -! do l=1,lm -! if(debugprint)print*,'sample ',VarName,'isa,jsa,l =' & -! ,cfr(isa,jsa,l),isa,jsa,l -! enddo - VarName='refl_10cm' call read_netcdf_3d_para(ncid2d,im,jm,jsta,jsta_2l,jend,jend_2u, & spval,VarName,REF_10CM(1,jsta_2l,1),lm) @@ -1839,12 +1831,13 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) VarName='albdo_ave' call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & spval,VarName,avgalbedo) - if(debugprint)print*,'sample ',VarName,' = ',avgalbedo(isa,jsa) +!$omp parallel do private(i,j) do j=jsta,jend do i=1,im if (avgalbedo(i,j) /= spval) avgalbedo(i,j) = avgalbedo(i,j) * 0.01 enddo enddo + if(debugprint)print*,'sample ',VarName,' = ',avgalbedo(isa,jsa) ! surface potential T using getgb VarName='tmpsfc' @@ -1877,6 +1870,12 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) enddo if(debugprint)print*,'sample ',VarName,' = ',ths(isa,jsa) +! foundation temperature + VarName='tref' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,fdnsst) + if(debugprint)print*,'sample ',VarName,' = ',fdnsst(isa,jsa) + ! GFS does not have time step and physics time step, make up ones since they ! are not really used anyway @@ -1960,11 +1959,18 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) enddo ! convective precip rate in m per physics time step -! VarName='cnvprcp' + VarName='cnvprcp' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,cprate) !set cprate as 0. do j=jsta,jend do i=1,im - cprate(i,j) = 0. + if (cprate(i,j) /= spval) then + cprate(i,j) = max(0.,cprate(i,j)) * (dtq2*0.001) & + * 1000. / dtp + else + cprate(i,j) = 0. + endif enddo enddo @@ -2002,13 +2008,6 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) spval,VarName,wspd10max) if(debugprint)print*,'sample ',VarName,' = ',wspd10max(isa,jsa) - -! 2m T using nemsio - VarName='tmp2m' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,tshltr) - if(debugprint)print*,'sample ',VarName,' = ',tshltr(isa,jsa) - ! inst snow water eqivalent using nemsio VarName='weasd' call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & @@ -2076,19 +2075,6 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) spval,VarName,qshltr) if(debugprint)print*,'sample ',VarName,' = ',qshltr(isa,jsa) -! mid day avg albedo in fraction using nemsio - VarName='albdosfc' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,avgalbedo) -! where(avgalbedo /= spval)avgalbedo=avgalbedo/100. ! convert to fraction -!$omp parallel do private(i,j) - do j=jsta,jend - do i=1,im - if (avgalbedo(i,j) /= spval) avgalbedo(i,j) = avgalbedo(i,j) * 0.01 - enddo - enddo - if(debugprint)print*,'sample ',VarName,' = ',avgalbedo(isa,jsa) - ! time averaged column cloud fractionusing nemsio VarName='tcdc_aveclm' call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & @@ -2115,14 +2101,11 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) VarName='snoalb' call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & spval,VarName,mxsnal) -! where(mxsnal /= spval)mxsnal=mxsnal/100. ! convert to fraction -!!$omp parallel do private(i,j) -! do j=jsta,jend -! do i=1,im -! if (mxsnal(i,j) /= spval) mxsnal(i,j) = mxsnal(i,j) * 0.01 -! enddo -! enddo -! if(debugprint)print*,'sample ',VarName,' = ',mxsnal(isa,jsa) + +! land fraction + VarName='lfrac' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,landfrac) ! GFS probably does not use sigt4, set it to sig*t^4 !$omp parallel do private(i,j,tlmh) @@ -2652,12 +2635,6 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & spval,VarName,rswin) -! time averaged model top incoming shortwave - VarName='dswrf_avetoa' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,aswintoa) -! if(debugprint)print*,'sample l',VarName,' = ',1,aswintoa(isa,jsa) - ! inst incoming clear sky sfc shortwave ! FV3 do not output instant incoming clear sky sfc shortwave !$omp parallel do private(i,j) @@ -2697,6 +2674,12 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & spval,VarName,rswout) +! time averaged model top incoming shortwave + VarName='dswrf_avetoa' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,aswintoa) +! if(debugprint)print*,'sample l',VarName,' = ',1,aswintoa(isa,jsa) + ! time averaged model top outgoing shortwave VarName='uswrf_avetoa' call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & @@ -3054,63 +3037,63 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) end do if(debugprint)print*,'sample hbot = ',hbot(isa,jsa) ! retrieve time averaged low cloud top pressure using nemsio - VarName='pres_ave' - VcoordName='low cld top' - l=1 + VarName='pres_avelct' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,ptopl) ! if(debugprint)print*,'sample l',VarName,' = ',1,ptopl(isa,jsa) ! retrieve time averaged low cloud bottom pressure using nemsio - VarName='pres_ave' - VcoordName='low cld bot' - l=1 + VarName='pres_avelcb' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,pbotl) ! if(debugprint)print*,'sample l',VarName,' = ',1,pbotl(isa,jsa) ! retrieve time averaged low cloud top temperature using nemsio - VarName='tmp_ave' - VcoordName='low cld top' - l=1 + VarName='tmp_avelct' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,Ttopl) ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopl(isa,jsa) ! retrieve time averaged middle cloud top pressure using nemsio - VarName='pres_ave' - VcoordName='mid cld top' - l=1 + VarName='pres_avemct' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,ptopm) ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptopm(isa,jsa) ! retrieve time averaged middle cloud bottom pressure using nemsio - VarName='pres_ave' - VcoordName='mid cld bot' - l=1 + VarName='pres_avemcb' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,pbotm) ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pbotm(isa,jsa) ! retrieve time averaged middle cloud top temperature using nemsio - VarName='tmp_ave' - VcoordName='mid cld top' - l=1 + VarName='tmp_avemct' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,Ttopm) ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttopm(isa,jsa) ! retrieve time averaged high cloud top pressure using nemsio ********* - VarName='pres_ave' - VcoordName='high cld top' - l=1 + VarName='pres_avehct' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,ptoph) ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,ptoph(isa,jsa) ! retrieve time averaged high cloud bottom pressure using nemsio - VarName='pres_ave' - VcoordName='high cld bot' - l=1 + VarName='pres_avehcb' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,pboth) ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,pboth(isa,jsa) ! retrieve time averaged high cloud top temperature using nemsio - VarName='tmp_ave' - VcoordName='high cld top' - l=1 + VarName='tmp_avehct' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,Ttoph) ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ',1,Ttoph(isa,jsa) ! retrieve boundary layer cloud cover using nemsio VarName='tcdc_avebndcl' - VcoordName='bndary-layer cld' - l=1 + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,pblcfr) ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,pblcfr(isa,jsa) ! where (pblcfr /= spval)pblcfr=pblcfr/100. ! convert to fraction !$omp parallel do private(i,j) @@ -3121,7 +3104,7 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) enddo ! retrieve cloud work function - VarName='cwork_ave' + VarName='cwork_aveclm' call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & spval,VarName,cldwork) ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,cldwork(isa,jsa) @@ -3137,15 +3120,65 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) if (sm(i,j) /= 0.0) runoff(i,j) = spval enddo enddo + +! total water storage in aquifer + VarName='wa_acc' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,twa) +! mask water areas +!$omp parallel do private(i,j) + do j=jsta,jend + do i=1,im + if (sm(i,j) /= 0.0) twa(i,j) = spval + enddo + enddo ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', 1,runoff(isa,jsa) - + +! accumulated evaporation of intercepted water + VarName='ecan_acc' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,tecan) +! mask water areas +!$omp parallel do private(i,j) + do j=jsta,jend + do i=1,im + if (sm(i,j) /= 0.0) tecan(i,j) = spval + enddo + enddo + +! accumulated plant transpiration + VarName='etran_acc' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,tetran) +! mask water areas +!$omp parallel do private(i,j) + do j=jsta,jend + do i=1,im + if (sm(i,j) /= 0.0) tetran(i,j) = spval + enddo + enddo + +! accumulated soil surface evaporation + VarName='edir_acc' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,tedir) +! mask water areas +!$omp parallel do private(i,j) + do j=jsta,jend + do i=1,im + if (sm(i,j) /= 0.0) tedir(i,j) = spval + enddo + enddo + ! retrieve shelter max temperature using nemsio VarName='t02max' + if(modelname=='GFS') VarName='tmax_max2m' call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & spval,VarName,maxtshltr) ! retrieve shelter min temperature using nemsio VarName='t02min' + if(modelname=='GFS') VarName='tmin_min2m' call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & spval,VarName,mintshltr) ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', & @@ -3162,15 +3195,19 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) spval,VarName,minrhshltr) ! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', & ! 1,mintshltr(im/2,(jsta+jend)/2) + +! retrieve shelter max specific humidity using nemsio + VarName='spfhmax_max2m' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,maxqshltr) +! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', +! 1,maxqshltr(isa,jsa) + +! retrieve shelter min temperature using nemsio + VarName='spfhmin_min2m' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,minqshltr) -!!$omp parallel do private(i,j) -! do j=jsta_2l,jend_2u -! do i=1,im -! MAXRHSHLTR(i,j) = SPVAL -! MINRHSHLTR(i,j) = SPVAL -! enddo -! enddo - ! retrieve ice thickness using nemsio VarName='icetk' call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & @@ -3260,18 +3297,6 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & spval,VarName,aswinc) -! retrieve shelter max specific humidity using nemsio - VarName='spfhmax_max2m' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,maxqshltr) -! if(debugprint)print*,'sample l',VcoordName,VarName,' = ', -! 1,maxqshltr(isa,jsa) - -! retrieve shelter min temperature using nemsio - VarName='spfhmin_min2m' - call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & - spval,VarName,minqshltr) - ! retrieve storm runoff using nemsio VarName='ssrun_acc' call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & @@ -3308,6 +3333,30 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) enddo enddo +! retrieve AVERAGED PRECIP ADVECTED HEAT FLUX + VarName='pah_ave' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,paha) +! mask water areas +!$omp parallel do private(i,j) + do j=jsta,jend + do i=1,im + if (sm(i,j) /= 0.0) paha(i,j) = spval + enddo + enddo + +! retrieve instantaneous PRECIP ADVECTED HEAT FLUX + VarName='pahi' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,pahi) +! mask water areas +!$omp parallel do private(i,j) + do j=jsta,jend + do i=1,im + if (sm(i,j) /= 0.0) pahi(i,j) = spval + enddo + enddo + ! retrieve PLANT TRANSPIRATION VarName='trans_ave' call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & @@ -3355,6 +3404,11 @@ SUBROUTINE INITPOST_NETCDF(ncid2d,ncid3d) if (sm(i,j) /= 0.0) snopcx(i,j) = spval enddo enddo + +! retrieve pwater + VarName='pwat' + call read_netcdf_2d_para(ncid2d,im,jsta,jsta_2l,jend,jend_2u, & + spval,VarName,pwat) ! GFS does not have deep convective cloud top and bottom fields @@ -3613,3 +3667,93 @@ subroutine read_netcdf_2d_scatter(me,ncid,ifhr,im,jm,jsta,jsta_2l & ,buf(1,jsta),icnt(me),mpi_real,0,MPI_COMM_COMP,iret) end subroutine read_netcdf_2d_scatter + + subroutine read_netcdf_3d_para(ncid,im,jm,jsta,jsta_2l,jend,jend_2u, & + spval,varname,buf,lm) + + use netcdf + use ctlblk_mod, only : me + use params_mod, only : small + implicit none + INCLUDE "mpif.h" + + character(len=20),intent(in) :: varname + real,intent(in) :: spval + integer,intent(in) :: ncid,im,jm,lm,jsta_2l,jend_2u,jsta,jend + real,intent(out) :: buf(im,jsta_2l:jend_2u,lm) + integer :: varid,iret,jj,i,j,l,kk + integer :: start(3), count(3), stride(3) + real,parameter :: spval_netcdf=9.99e+20 + real :: fill_value + + iret = nf90_inq_varid(ncid,trim(varname),varid) + if (iret /= 0) then + if (me == 0) print*,VarName," not found -Assigned missing values" +!$omp parallel do private(i,j,l) + do l=1,lm + do j=jsta,jend + do i=1,im + buf(i,j,l)=spval + enddo + enddo + enddo + else + iret = nf90_get_att(ncid,varid,"_FillValue",fill_value) + if (iret /= 0) fill_value = spval_netcdf + start = (/1,jsta,1/) + jj=jend-jsta+1 + count = (/im,jj,lm/) + iret = nf90_get_var(ncid,varid,buf(1:im,jsta:jend,1:lm),start=start,count=count) + do l=1,lm + do j=jsta,jend + do i=1,im + if(abs(buf(i,j,l)-fill_value)