module gsi_nemsio_mod !$$$ module documentation block ! . . . . ! module: gsi_nemsio_mod ! prgmmr: ! ! abstract: ! ! program history log: ! 2009-08-04 lueken - added module doc block ! 2014-06-30 wu - remove debugging printout ! 2015_05_13 wu - output error flag of nemsio_open ! 2015-06-10 s.liu - add gsi_nemsio_read_fraction to handle NMMB f_rain and f_ice ! 2015-06-10 s.liu - add gsi_nemsio_write_fraction to handle NMMB f_rain and f_ice ! 2016-02-05 s.liu - add fraction2variable and variable2fraction to handle NMMB f_rain and f_ice ! ! subroutines included: ! sub gsi_nemsio_open ! sub gsi_nemsio_update ! sub gsi_nemsio_close ! sub gsi_nemsio_read ! sub gsi_nemsio_read_fraction ! sub gsi_nemsio_write ! sub gsi_nemsio_write_fraction ! sub fraction2variable ! sub variable2fraction ! ! variable definitions: ! ! attributes: ! langauge: f90 ! machine: ! !$$$ end documentation block use kinds, only: r_kind,i_kind,r_single use nemsio_module, only: nemsio_gfile use gridmod, only: nlon_regional,nlat_regional implicit none type(nemsio_gfile) :: gfile save gfile real(r_single),allocatable::work_saved(:) ! set default to private private ! set subroutines to public public :: gsi_nemsio_open public :: gsi_nemsio_update public :: gsi_nemsio_close public :: gsi_nemsio_read public :: gsi_nemsio_read_fraction public :: gsi_nemsio_read_fractionnew public :: gsi_nemsio_write public :: gsi_nemsio_write_fraction public :: gsi_nemsio_write_fractionnew contains subroutine gsi_nemsio_open(file_name,iostatus,message,mype,mype_io,ierr) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_open ! pgrmmr: ! ! abstract: ! ! program history log: ! 2009-08-04 lueken - added subprogram doc block ! ! input argument list: ! file_name ! iostatus ! message ! mype - mpi task id ! mype_io ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use nemsio_module, only: nemsio_init,nemsio_open implicit none character(*) ,intent(in ) :: file_name ! input file name character(*) ,intent(in ) :: iostatus ! 'READ' for read only, 'rdwr' for read/write character(*) ,intent(in ) :: message ! info to appear in write statement on status of file open integer(i_kind),intent(in ) :: mype,mype_io integer(i_kind),intent(out ) :: ierr integer(i_kind) iret if(mype==mype_io) then call nemsio_init(iret=iret) if(iret/=0) then write(6,*)trim(message),' problem with nemsio_init, Status = ',iret call stop2(74) end if ierr=0 call nemsio_open(gfile,file_name,trim(iostatus),iret=iret) if(iret/=0) then write(6,*)trim(message),' problem opening file',trim(file_name),', Status = ',iret ierr=1 return end if end if allocate(work_saved(nlon_regional*nlat_regional)) end subroutine gsi_nemsio_open subroutine gsi_nemsio_update(file_name,message,mype,mype_io) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_update ! pgrmmr: ! ! abstract: ! ! program history log: ! 2009-08-04 lueken - added subprogram doc block ! ! input argument list: ! file_name ! message ! mype - mpi task id ! mype_io ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use nemsio_module, only: nemsio_init,nemsio_open,nemsio_getfilehead,nemsio_close,nemsio_setheadvar use nemsio_module, only: nemsio_getheadvar use constants, only: zero implicit none character(*) ,intent(in ) :: file_name ! input file name character(*) ,intent(in ) :: message ! info to appear in write statement on status of file open integer(i_kind),intent(in ) :: mype,mype_io integer(i_kind) iret,nrec integer(i_kind) idate(7),jdate(7),nfhour,nfminute,nfsecondn,nfday,ihrst,idat(3) integer(i_kind),dimension(8):: ida,jda real(r_kind),dimension(5):: fha integer(i_kind) im,jm,lm,nfsecondd,nframe,ntrac,nsoil,nmeta,ntimestep logical extrameta character(4) gdatatype,modelname character(32) gtype if(mype==mype_io) then call nemsio_init(iret=iret) if(iret/=0) then write(6,*)trim(message),' problem with nemsio_init, Status = ',iret call stop2(74) end if call nemsio_open(gfile,file_name,'RDWR',iret=iret) if(iret/=0) then write(6,*)trim(message),' problem opening file',trim(file_name),', Status = ',iret call stop2(74) end if call nemsio_getheadvar(gfile,'idat',idat,iret) write(6,*)' check old idat after getheadvar, idat,iret=',idat,iret call nemsio_getheadvar(gfile,'ihrst',ihrst,iret) write(6,*)' check old ihrst after getheadvar, ihrst,iret=',ihrst,iret call nemsio_getheadvar(gfile,'ntimestep',ntimestep,iret) write(6,*)' check old ntimestep after getheadvar, ntimestep,iret=',ntimestep,iret call nemsio_getfilehead(gfile,iret=iret,nrec=nrec,dimx=im,dimy=jm, & dimz=lm,idate=idate,gdatatype=gdatatype,gtype=gtype,modelname=modelname, & nfhour=nfhour,nfminute=nfminute,nfsecondn=nfsecondn,nfsecondd=nfsecondd, & nfday=nfday, & nframe=nframe,ntrac=ntrac,nsoil=nsoil,extrameta=extrameta,nmeta=nmeta) write(6,*)' in gsi_nemsio_update, guess yr,mn,dy,hr,fhr=',idate(1:4),nfhour fha=zero ; ida=0 ; jda=0 fha(2)=nfhour fha(3)=nfminute ida(1)=idate(1) ! year ida(2)=idate(2) ! month ida(3)=idate(3) ! day ida(4)=0 ! time zone ida(5)=idate(4) ! hour ida(6)=idate(5) ! minute call w3movdat(fha,ida,jda) jdate(1)=jda(1) ! new year jdate(2)=jda(2) ! new month jdate(3)=jda(3) ! new day jdate(4)=jda(5) ! new hour jdate(5)=jda(6) ! new minute jdate(6)=0 ! new scaled seconds jdate(7)=idate(7) ! new seconds multiplier nfhour=0 ! new forecast hour nfminute=0 nfsecondn=0 ntimestep=0 call nemsio_setheadvar(gfile,'idate',jdate,iret) write(6,*)' after setheadvar, jdate,iret=',jdate,iret call nemsio_setheadvar(gfile,'nfhour',nfhour,iret) write(6,*)' after setheadvar, nfhour,iret=',nfhour,iret call nemsio_setheadvar(gfile,'nfminute',nfminute,iret) write(6,*)' after setheadvar, nfminute,iret=',nfminute,iret call nemsio_setheadvar(gfile,'nfsecondn',nfsecondn,iret) write(6,*)' after setheadvar, nfsecondn,iret=',nfsecondn,iret ! idat(3)=jdate(1) ! forecast starting year idat(2)=jdate(2) ! forecast starting month idat(1)=jdate(3) ! forecast starting day ihrst=jdate(4) ! forecast starting hour (0-23) call nemsio_setheadvar(gfile,'idat',idat,iret) write(6,*)' after setheadvar, idat,iret=',idat,iret call nemsio_setheadvar(gfile,'ihrst',ihrst,iret) write(6,*)' after setheadvar, ihrst,iret=',ihrst,iret call nemsio_setheadvar(gfile,'ntimestep',ntimestep,iret) write(6,*)' after setheadvar, ntimestep,iret=',ntimestep,iret ! Following is diagnostic to check if date updated: call nemsio_getfilehead(gfile,iret=iret,nrec=nrec,dimx=im,dimy=jm, & dimz=lm,idate=idate,gdatatype=gdatatype,gtype=gtype,modelname=modelname, & nfhour=nfhour,nfminute=nfminute,nfsecondn=nfsecondn,nfsecondd=nfsecondd, & nfday=nfday, & nframe=nframe,ntrac=ntrac,nsoil=nsoil,extrameta=extrameta,nmeta=nmeta) write(6,*)' in gsi_nemsio_update, analysis yr,mn,dy,hr,fhr=',idate(1:4),nfhour call nemsio_getheadvar(gfile,'idat',idat,iret) write(6,*)' check new idat after getheadvar, idat,iret=',idat,iret call nemsio_getheadvar(gfile,'ihrst',ihrst,iret) write(6,*)' check new ihrst after getheadvar, ihrst,iret=',ihrst,iret call nemsio_getheadvar(gfile,'ntimestep',ntimestep,iret) write(6,*)' check new ntimestep after getheadvar, ntimestep,iret=',ntimestep,iret call nemsio_close(gfile,iret=iret) if(iret/=0) then write(6,*)trim(message),' problem closing file',trim(file_name),', Status = ',iret call stop2(74) end if end if end subroutine gsi_nemsio_update subroutine gsi_nemsio_close(file_name,message,mype,mype_io) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_close ! pgrmmr: ! ! abstract: ! ! program history log: ! 2009-08-04 lueken - added subprogram doc block ! ! input argument list: ! file_name ! message ! mype - mpi task id ! mype_io ! ! output argument list: ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use nemsio_module, only: nemsio_close implicit none character(*) ,intent(in ) :: file_name ! input file name character(*) ,intent(in ) :: message ! info to appear in write statement on status of file open integer(i_kind),intent(in ) :: mype,mype_io integer(i_kind) iret if(mype==mype_io) then call nemsio_close(gfile,iret=iret) if(iret/=0) then write(6,*)trim(message),' problem closing file',trim(file_name),', Status = ',iret call stop2(74) end if end if deallocate(work_saved) end subroutine gsi_nemsio_close subroutine gsi_nemsio_read(varname,vartype,gridtype,lev,var,mype,mype_io,good_var) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_read ! pgrmmr: parrish ! ! abstract: intermediate level routine to read nmmb model fields using nems_io. ! the desired field is retrieved from the previously opened file as a ! full 2d horizontal field, then interpolated to the analysis grid ! from the nmmb model grid. finally, the 2d field is scattered from ! processor mype_io to subdomains in output array var. ! a copy of the original field on the nmmb grid is saved internally in array ! work_saved in case this field is to be updated by the analysis ! increment in a call to gsi_nemsio_write immediately after the call to ! gsi_nemsio_read. ! ! program history log: ! 2009-08-04 lueken - added subprogram doc block ! 2010-01-22 parrish - added optional variable good_var to detect read errors in calling program ! and have option to avoid program stop. ! 2013-10-25 todling - reposition ltosi and others to commvars ! ! input argument list: ! varname,vartype,gridtype - descriptors for variable to be retrieved from nmmb file ! lev - vertical level number ! mype - mpi task id ! mype_io - mpi task where field is read from disk ! good_var - optional, on input, set to .false. if present(good_var) then error stop is ! bypassed and good_var is returned .true. for successful read, .false. otherwise. ! ! output argument list: ! var - for successful read, contains desired variable on subdomains. ! good_var - see above ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use mpimod, only: mpi_rtype,mpi_comm_world,ierror,mpi_integer4 use gridmod, only: lat2,lon2,nlon,nlat use gridmod, only: ijn_s,displs_s,itotsub use general_commvars_mod, only: ltosi_s,ltosj_s use nemsio_module, only: nemsio_readrecv use mod_nmmb_to_a, only: nmmb_h_to_a,nmmb_v_to_a implicit none character(*) ,intent(in ) :: varname,vartype,gridtype ! gridtype='H' or 'V' integer(i_kind),intent(in ) :: lev ! vertical level of desired variable real(r_kind) ,intent( out) :: var(lat2*lon2) integer(i_kind),intent(in ) :: mype,mype_io logical,optional,intent(inout):: good_var integer(i_kind) i,iret,j,mm1,n real(r_kind) work(itotsub) real(r_kind) work_a(nlat,nlon) real(r_single) work_b(nlon_regional*nlat_regional) logical good_var_loc mm1=mype+1 if(mype==mype_io) then ! read field from file with nemsio call nemsio_readrecv(gfile,trim(varname),trim(vartype),lev,work_b,iret=iret) if(iret==0) then work_saved=work_b ! interpolate to analysis grid if(trim(gridtype)=='H') call nmmb_h_to_a(work_b,work_a) if(trim(gridtype)=='V') call nmmb_v_to_a(work_b,work_a) ! scatter to subdomains do n=1,itotsub i=ltosi_s(n) j=ltosj_s(n) work(n)=work_a(i,j) end do end if end if call mpi_bcast(iret,1,mpi_integer4,mype_io,mpi_comm_world,ierror) good_var_loc=.true. if(iret/=0) then good_var_loc=.false. if(mype==0) then write(6,*)' problem reading varname=',trim(varname),', vartype=',trim(vartype),', Status = ',iret if(.not.present(good_var)) call stop2(74) end if end if if(present(good_var)) good_var=good_var_loc if(good_var_loc) & call mpi_scatterv(work,ijn_s,displs_s,mpi_rtype, & var,ijn_s(mm1),mpi_rtype,mype_io,mpi_comm_world,ierror) end subroutine gsi_nemsio_read subroutine gsi_nemsio_read_fraction(varname_frain,varname_fice,varname_clwmr,varname_t, & vartype,lev,var_qi,var_qs,var_qr,var_qw,mype,mype_io,good_var) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_read_fraction ! pgrmmr: Shun Liu ! ! abstract: copy from gsi_nemsio_read. To read in NMMB f_rain, f_ice, f_rime and ! T together and then convert to rain water mixing ratio and snow ! mixing ratio ! ! program history log: ! 2015-06-5 S.Liu - read in f_rain, f_ice, f_rimef and T ! 2016-02-10 S.Liu - remove gridtype if-test since all variables are in mass point ! ! input argument list: ! varname,vartype,gridtype - descriptors for variable to be retrieved from ! nmmb file ! lev - vertical level number ! mype - mpi task id ! mype_io - mpi task where field is read from disk ! good_var - optional, on input, set to .false. if present(good_var) then ! error stop is ! bypassed and good_var is returned .true. for successful read, ! .false. otherwise. ! ! output argument list: ! var - for successful read, contains desired variable on subdomains. ! good_var - see above ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use mpimod, only: mpi_rtype,mpi_comm_world,ierror,mpi_integer4 use gridmod, only: lat2,lon2,nlon,nlat use gridmod, only: ijn_s,displs_s,itotsub use general_commvars_mod, only: ltosi_s,ltosj_s use nemsio_module, only: nemsio_readrecv use mod_nmmb_to_a, only: nmmb_h_to_a,nmmb_v_to_a implicit none character(*) ,intent(in ) :: vartype ! gridtype='H' or 'V' character(*) ,intent(in ) :: varname_frain, varname_fice, varname_clwmr, varname_t ! gridtype='H' or 'V' integer(i_kind),intent(in ) :: lev ! vertical level of desired variable real(r_kind) ,intent( out) :: var_qi(lat2*lon2) real(r_kind) ,intent( out) :: var_qs(lat2*lon2) real(r_kind) ,intent( out) :: var_qr(lat2*lon2) real(r_kind) ,intent( out) :: var_qw(lat2*lon2) integer(i_kind),intent(in ) :: mype,mype_io logical,optional,intent(inout):: good_var integer(i_kind) i,iret,j,mm1,n real(r_kind) work_qi(itotsub) real(r_kind) work_qs(itotsub) real(r_kind) work_qr(itotsub) real(r_kind) work_qw(itotsub) real(r_kind) work_a_qi(nlat,nlon) real(r_kind) work_a_qs(nlat,nlon) real(r_kind) work_a_qr(nlat,nlon) real(r_kind) work_a_qw(nlat,nlon) real(r_single) work_b_frain(nlon_regional*nlat_regional) real(r_single) work_b_fice(nlon_regional*nlat_regional) real(r_single) work_b_clwmr(nlon_regional*nlat_regional) real(r_single) work_b_t(nlon_regional*nlat_regional) real(r_single) work_b_qi(nlon_regional*nlat_regional) real(r_single) work_b_qs(nlon_regional*nlat_regional) real(r_single) work_b_qr(nlon_regional*nlat_regional) real(r_single) work_b_qw(nlon_regional*nlat_regional) real(r_single) :: t, f_ice, f_rain, wc, qi, qs, qr, qw logical good_var_loc mm1=mype+1 if(mype==mype_io) then ! read field from file with nemsio call nemsio_readrecv(gfile,trim(varname_frain),trim(vartype),lev,work_b_frain,iret=iret) call nemsio_readrecv(gfile,trim(varname_fice),trim(vartype),lev,work_b_fice,iret=iret) call nemsio_readrecv(gfile,trim(varname_clwmr),trim(vartype),lev,work_b_clwmr,iret=iret) call nemsio_readrecv(gfile,trim(varname_t),trim(vartype),lev,work_b_t,iret=iret) do n=1,nlon_regional*nlat_regional t=work_b_t(n) f_rain=work_b_frain(n) f_ice=work_b_fice(n) wc=work_b_clwmr(n) call fraction2variable(t,f_ice,f_rain,wc,qi,qs,qr,qw) work_b_qi(n)=qi work_b_qs(n)=qs work_b_qr(n)=qr work_b_qw(n)=qw end do if(iret==0) then ! work_saved=work_b ! interpolate to analysis grid call nmmb_h_to_a(work_b_qi,work_a_qi) call nmmb_h_to_a(work_b_qs,work_a_qs) call nmmb_h_to_a(work_b_qr,work_a_qr) call nmmb_h_to_a(work_b_qw,work_a_qw) ! scatter to subdomains do n=1,itotsub i=ltosi_s(n) j=ltosj_s(n) work_qi(n)=work_a_qi(i,j) work_qs(n)=work_a_qs(i,j) work_qr(n)=work_a_qr(i,j) work_qw(n)=work_a_qw(i,j) end do end if end if call mpi_bcast(iret,1,mpi_integer4,mype_io,mpi_comm_world,ierror) good_var_loc=.true. if(iret/=0) then good_var_loc=.false. if(mype==0) then write(6,*)' problem reading varname=',trim(varname_frain),', vartype=',trim(vartype),', Status = ',iret if(.not.present(good_var)) call stop2(74) end if end if if(present(good_var)) good_var=good_var_loc if(good_var_loc) then call mpi_scatterv(work_qi,ijn_s,displs_s,mpi_rtype, & var_qi,ijn_s(mm1),mpi_rtype,mype_io,mpi_comm_world,ierror) call mpi_scatterv(work_qs,ijn_s,displs_s,mpi_rtype, & var_qs,ijn_s(mm1),mpi_rtype,mype_io,mpi_comm_world,ierror) call mpi_scatterv(work_qr,ijn_s,displs_s,mpi_rtype, & var_qr,ijn_s(mm1),mpi_rtype,mype_io,mpi_comm_world,ierror) call mpi_scatterv(work_qw,ijn_s,displs_s,mpi_rtype, & var_qw,ijn_s(mm1),mpi_rtype,mype_io,mpi_comm_world,ierror) end if end subroutine gsi_nemsio_read_fraction subroutine gsi_nemsio_write(varname,vartype,gridtype,lev,var,mype,mype_io,add_saved) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_write ! pgrmmr: ! ! abstract: ! ! program history log: ! 2009-08-04 lueken - added subprogram doc block ! 2013-10-25 todling - reposition ltosi and others to commvars ! ! input argument list: ! varname,vartype,gridtype ! lev ! add_saved ! mype - mpi task id ! mype_io ! ! output argument list: ! var ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use mpimod, only: mpi_rtype,mpi_comm_world,ierror use gridmod, only: lat2,lon2,nlon,nlat,lat1,lon1 use gridmod, only: ijn,displs_g,itotsub,iglobal use general_commvars_mod, only: ltosi,ltosj use nemsio_module, only: nemsio_writerecv use mod_nmmb_to_a, only: nmmb_a_to_h,nmmb_a_to_v implicit none character(*) ,intent(in ) :: varname,vartype,gridtype ! gridtype='H' or 'V' integer(i_kind),intent(in ) :: lev ! vertical level of desired variable real(r_kind) ,intent(in ) :: var(lat2,lon2) integer(i_kind),intent(in ) :: mype,mype_io logical ,intent(in ) :: add_saved integer(i_kind) i,iret,j,mm1,n real(r_kind) work(itotsub),work_sub(lat1,lon1) real(r_kind) work_a(nlat,nlon) real(r_single) work_b(nlon_regional*nlat_regional) mm1=mype+1 do i=1,lon1 do j=1,lat1 work_sub(j,i)=var(j+1,i+1) end do end do call mpi_gatherv(work_sub,ijn(mm1),mpi_rtype, & work,ijn,displs_g,mpi_rtype,mype_io,mpi_comm_world,ierror) if(mype==mype_io) then do n=1,iglobal i=ltosi(n) j=ltosj(n) work_a(i,j)=work(n) end do if(trim(gridtype)=='H') call nmmb_a_to_h(work_a,work_b) if(trim(gridtype)=='V') call nmmb_a_to_v(work_a,work_b) if(add_saved) work_b=work_b+work_saved call nemsio_writerecv(gfile,trim(varname),trim(vartype),lev,work_b,iret=iret) if(iret/=0) then write(6,*)' problem writing varname=',trim(varname),', vartype=',trim(vartype),', Status = ',iret call stop2(74) end if end if end subroutine gsi_nemsio_write subroutine gsi_nemsio_write_fraction(varname_frain,varname_fice,vartype,lev,var_t,var_i,var_r,var_l,mype,mype_io) !$$$ subprogram documentation block ! . . . . ! subprogram: gsi_nemsio_write_fraction ! pgrmmr: Shun Liu ! ! abstract: ! ! program history log: ! 2015-05-12 S.Liu - copy from gsi_nemsio_write and modify to handle NMMB hydrometor fraction variable ! ! input argument list: ! varname,vartype,gridtype ! lev ! add_saved ! mype - mpi task id ! mype_io ! ! output argument list: ! var_frain, var_fice ! ! attributes: ! language: f90 ! machine: ! !$$$ end documentation block use mpimod, only: mpi_rtype,mpi_comm_world,ierror use gridmod, only: lat2,lon2,nlon,nlat,lat1,lon1 use gridmod, only: ijn,displs_g,itotsub,iglobal use general_commvars_mod, only: ltosi,ltosj use nemsio_module, only: nemsio_writerecv use mod_nmmb_to_a, only: nmmb_a_to_h,nmmb_a_to_v implicit none character(*) ,intent(in ) :: varname_frain,varname_fice,vartype ! gridtype='H' or 'V' integer(i_kind),intent(in ) :: lev ! vertical level of desired variable real(r_kind) ,intent(in ) :: var_i(lat2,lon2), var_r(lat2,lon2), var_l(lat2,lon2), var_t(lat2,lon2) integer(i_kind),intent(in ) :: mype,mype_io ! logical ,intent(in ) :: add_saved integer(i_kind) i,iret,j,mm1,n real(r_kind) work_t(itotsub),work_sub_t(lat1,lon1) real(r_kind) work_a_t(nlat,nlon) real(r_single) work_b_t(nlon_regional*nlat_regional) real(r_kind) work_i(itotsub),work_sub_i(lat1,lon1) real(r_kind) work_a_i(nlat,nlon) real(r_single) work_b_i(nlon_regional*nlat_regional) real(r_kind) work_r(itotsub),work_sub_r(lat1,lon1) real(r_kind) work_a_r(nlat,nlon) real(r_single) work_b_r(nlon_regional*nlat_regional) real(r_kind) work_l(itotsub),work_sub_l(lat1,lon1) real(r_kind) work_a_l(nlat,nlon) real(r_single) work_b_l(nlon_regional*nlat_regional) real(r_single) work_b_frain(nlon_regional*nlat_regional) real(r_single) work_b_fice(nlon_regional*nlat_regional) real(r_single) t,qfi,qfr,qfw,f_rain,f_ice mm1=mype+1 do i=1,lon1 do j=1,lat1 work_sub_t(j,i)=var_t(j+1,i+1) work_sub_i(j,i)=var_i(j+1,i+1) work_sub_r(j,i)=var_r(j+1,i+1) work_sub_l(j,i)=var_l(j+1,i+1) end do end do ! write(6,*)'writeout1', maxval(work_sub_t),maxval(work_sub_i),maxval(work_sub_r) call mpi_gatherv(work_sub_t,ijn(mm1),mpi_rtype, & work_t,ijn,displs_g,mpi_rtype,mype_io,mpi_comm_world,ierror) call mpi_gatherv(work_sub_i,ijn(mm1),mpi_rtype, & work_i,ijn,displs_g,mpi_rtype,mype_io,mpi_comm_world,ierror) call mpi_gatherv(work_sub_r,ijn(mm1),mpi_rtype, & work_r,ijn,displs_g,mpi_rtype,mype_io,mpi_comm_world,ierror) call mpi_gatherv(work_sub_l,ijn(mm1),mpi_rtype, & work_l,ijn,displs_g,mpi_rtype,mype_io,mpi_comm_world,ierror) ! write(6,*)'writeout2', maxval(work_t),maxval(work_i),maxval(work_r) if(mype==mype_io) then do n=1,iglobal i=ltosi(n) j=ltosj(n) work_a_t(i,j)=work_t(n) work_a_i(i,j)=work_i(n) work_a_r(i,j)=work_r(n) work_a_l(i,j)=work_l(n) end do ! write(6,*)'writeout3', maxval(work_a_r),maxval(work_a_l) call nmmb_a_to_h(work_a_t,work_b_t) call nmmb_a_to_h(work_a_i,work_b_i) call nmmb_a_to_h(work_a_r,work_b_r) call nmmb_a_to_h(work_a_l,work_b_l) ! if(trim(gridtype)=='V') call nmmb_a_to_v(work_a_t,work_b_i) ! if(trim(gridtype)=='V') call nmmb_a_to_v(work_a_i,work_b_i) ! if(trim(gridtype)=='V') call nmmb_a_to_v(work_a_r,work_b_r) ! if(trim(gridtype)=='V') call nmmb_a_to_v(work_a_l,work_b_l) ! if(add_saved) work_b_t=work_b_t+work_saved_t ! if(add_saved) work_b_i=work_b_i+work_saved_i ! if(add_saved) work_b_r=work_b_r+work_saved_r ! if(add_saved) work_b_l=work_b_l+work_saved_l ! write(6,*)'writeout4', maxval(work_b_r),maxval(work_b_l) ! write(6,*)'writeout44',nlon_regional,nlat_regional,nlon,nlat do n=1,nlon_regional*nlat_regional t=work_b_t(n) qfi=work_b_i(n) qfr=work_b_r(n) qfw=work_b_l(n) call variable2fraction(t, qfi, qfr, qfw, f_ice, f_rain) work_b_frain(n)=f_rain work_b_fice(n)=f_ice ! work_b_frain(n)=qfr ! work_b_fice(n)=qfw end do call nemsio_writerecv(gfile,trim(varname_frain),trim(vartype),lev,work_b_frain,iret=iret) call nemsio_writerecv(gfile,trim(varname_fice),trim(vartype),lev,work_b_fice,iret=iret) ! write(6,*)'writeout5', maxval(work_b_frain),maxval(work_b_fice) if(iret/=0) then write(6,*)' problem writing varname=',trim(varname_frain),', vartype=',trim(vartype),', Status = ',iret call stop2(74) end if end if end subroutine gsi_nemsio_write_fraction Subroutine fraction2variable(t,f_ice,f_rain, wc, qi,qs,qr,qw) !$$$ subprogram documentation block ! . . . . ! subprogram: gsdcloudanalysis driver for generalized cloud/hydrometeor ! analysis ! ! PRGMMR: Shun Liu ORG: EMC/NCEP DATE: 2015-05-28 ! ! ABSTRACT: ! This subroutine fraction to qi, qs, qr, qw ! ! PROGRAM HISTORY LOG: ! 2015-05-28 Shun Liu Add NCO document block ! 2016-06-21 Shun Liu give number precisio and remove f_rimef ! ! ! input argument list: ! mype - processor ID that does this IO ! ! output argument list: ! ! USAGE: ! INTPUT: ! t - sensible temperature ! f_ice - fraction of condensate in form of ice ! f_rain - fraction of liquid water in form of rain ! f_rimef - ratio of total ice growth to deposition groth ! OUTPUT ! qi - cloud ice mixing ratio ! qs - large ice mixing ratio ! qr - rain mixing ratio ! qw - cloud water mixing ratio ! ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 90 ! MACHINE: WCOSS at NOAA/ESRL - college park, DC ! !$$$ use kinds, only: r_kind,r_single real(r_single) t, qi,qs, qr, qw, wc real(r_single) f_ice, f_rain real(r_single),parameter:: epsq=1.e-12_r_single real(r_single),parameter:: tice=233.15_r_single,ticek=273.15_r_single real(r_single),parameter:: tice_mix=243.15_r_single real(r_single) ::t1,t2, coef1, coef2, coef qi=0.0_r_single; qs=0.0_r_single; qr=0.0_r_single; qw=0.0_r_single if(wc > 0.0_r_single) then if(f_ice>1.0_r_single) f_ice=1.0_r_single if(f_ice<0.0_r_single) f_ice=0.0_r_single if(f_rain>1.0_r_single) f_rain=1.0_r_single if(f_rain<0.0_r_single) f_rain=0.0_r_single qi=0.05_r_single*wc*f_ice qs=0.95_r_single*wc*f_ice if(t<=tice_mix)then t1=tice_mix t2=tice coef1=0.05_r_single coef2=0.10_r_single coef=(t-t2)/(t1-t2)*coef1+(t-t1)/(t2-t1)*coef2 qi=coef*wc*f_ice qs=(1.0_r_single-coef)*wc*f_ice end if !* do not consider frime at the moment qr=wc*(1.0_r_single-f_ice)*f_rain qw=wc*(1.0_r_single-f_ice)*(1.0_r_single-f_rain) end if end subroutine fraction2variable subroutine variable2fraction(t, qi, qr, qw, f_ice, f_rain) !$$$ subprogram documentation block ! . . . . ! subprogram: gsdcloudanalysis driver for generalized cloud/hydrometeor analysis ! ! PRGMMR: Shun Liu ORG: EMC/NCEP DATE: 2012-10-24 ! ! ABSTRACT: ! This subroutine qi qr qw to fraction ! ! PROGRAM HISTORY LOG: ! 2013-10-18 Shun Liu Add NCO document block ! 2015-11-16 Shun Liu move from gsdcldanalysis4nmmb.F90 to this module ! 2016-06-21 Shun Liu give number precisio ! ! ! input argument list: ! mype - processor ID that does this IO ! ! output argument list: ! ! USAGE: ! INPUT ! qi - cloud ice mixing ratio ! qr - rain mixing ratio ! qw - cloud water mixing ratio ! OUTPUT: ! f_ice - fraction of condensate in form of ice ! f_rain - fraction of liquid water in form of rain ! f_rimef - ratio of total ice growth to deposition groth ! ! ! REMARKS: ! ! ATTRIBUTES: ! LANGUAGE: FORTRAN 90 ! MACHINE: WCOSS at NOAA/ESRL - college park, DC ! !$$$ use kinds, only: r_kind,r_single real(r_single) t, qi, qr, qw, wc, dum real(r_single) f_ice, f_rain real(r_single),parameter:: epsq=1.e-12_r_single real(r_single),parameter:: tice=233.15_r_single,ticek=273.15_r_single wc=qi+qr+qw if(wc > 0.0_r_single) then if(qi